Changeset 3272


Ignore:
Timestamp:
03/12/10 11:37:01 (15 years ago)
Author:
Mathieu Morlighem
Message:

minor

Location:
issm/trunk/src
Files:
9 edited

Legend:

Unmodified
Added
Removed

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

    r3263 r3272  
    3535
    3636                        //Methods
    37                         R2 F(double theta) const ; // parametrization of the curve edge
     37                        R2     F(double theta) const ; // parametrization of the curve edge
    3838                        double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente
    39                         int  Cracked() const {return flag & 1;}
    40                         int  Dup() const { return flag & 32;}
    41                         int  Equi()const {return flag & 2;}
    42                         int  ReverseEqui()const {return flag & 128;}
    43                         int  TgA()const {return flag &4;}
    44                         int  TgB()const {return flag &8;}
    45                         int  Tg(int i) const{return i==0 ? TgA() : TgB();}
    46                         int  Mark()const {return flag &16;}
    47                         int  Required() { return flag & 64;}
    48                         void SetCracked() { flag |= 1;}
    49                         void SetDup()     { flag |= 32;} // not a real edge
    50                         void SetEqui()    { flag |= 2;}
    51                         void SetTgA()     { flag|=4;}
    52                         void SetTgB()     { flag|=8;}
    53                         void SetMark()    { flag|=16;}
    54                         void SetUnMark()  { flag &= 1007 /* 1023-16*/;}
    55                         void SetRequired() { flag |= 64;}
    56                         void SetReverseEqui() {flag |= 128;}
    57                         void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew);
     39                        int    Tg(int i) const{return i==0 ? TgA() : TgB();}
     40                        int    Cracked() const    {return flag & 1;  }
     41                        int    Dup() const        {return flag & 32; }
     42                        int    Equi()const        {return flag & 2;  }
     43                        int    ReverseEqui()const {return flag & 128;}
     44                        int    TgA()const         {return flag & 4;  }
     45                        int    TgB()const         {return flag & 8;  }
     46                        int    Mark()const        {return flag & 16; }
     47                        int    Required()         {return flag & 64; }
     48                        void   SetCracked()     { flag |= 1;  }
     49                        void   SetDup()         { flag |= 32; } // not a real edge
     50                        void   SetEqui()        { flag |= 2;  }
     51                        void   SetTgA()         { flag |=4;   }
     52                        void   SetTgB()         { flag |=8;   }
     53                        void   SetMark()        { flag |=16;  }
     54                        void   SetUnMark()      { flag &= 1007 /* 1023-16*/;}
     55                        void   SetRequired()    { flag |= 64; }
     56                        void   SetReverseEqui() { flag |= 128;}
     57                        void   Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew);
    5858        };
    5959

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

    r3263 r3272  
    11
     2#include <assert.h>
    23#include "BamgObjects.h"
    34
    … …  
    1819                //initialize fields
    1920                nx   =nnx;   //number of vertices
    20                 nbax =mmx; // 3 * number of triangles
     21                nbax =mmx;   // 3 * number of triangles
    2122                NbOfEdges=0;
    2223                head = new long [nx];
    … …  
    2526                //initialize head (-1 everywhere)
    2627                i=nx;
    27                 while (i--) head[i]=-1;
     28                while(i--) head[i]=-1;
    2829        }
    2930        /*}}}1*/
    … …  
    6667                int h,n;
    6768
    68                 //check that head is not empty
    69                 if (head==NULL) {
    70                         throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add no head is NULL (How come?)"));
    71                 }
    72 
    7369                //get n from h (usually h=ii)
     70                assert(head);
    7471                n=head[h=Abs(ii)%nx];
    7572

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

    r3245 r3272  
    2121
    2222                public:
     23
     24                        // Constructors
    2325                        SetOfEdges4(long ,long);// nb Edges mx , nb de sommet
    2426                        ~SetOfEdges4() {delete [] head; delete [] Edges;}
     27
     28                        //operators
     29                        IntEdge & operator[](long k){return  Edges[k];}
     30
     31                        //Methods
    2532                        long add (long ii,long jj);
    2633                        long SortAndAdd (long ii,long jj) {return ii <=jj ? add (ii,jj)  : add (jj,ii) ;}
    27                         long  nb(){return NbOfEdges;}
     34                        long nb(){return NbOfEdges;}
    2835                        long find (long ii,long jj);
    2936                        long SortAndFind (long ii,long jj) {return ii <=jj ? find (ii,jj)  : find (jj,ii) ;}
    … …  
    3138                        long j(long k){return Edges[k].j;}
    3239                        long newarete(long k){return NbOfEdges == k+1;}
    33 
    34                         //operators
    35                         IntEdge & operator[](long k){return  Edges[k];}
    3640        };
    3741}

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

    r3263 r3272  
    7979                          }
    8080                        void SetAdj2(short a,Triangle *t,short aat){
    81                                 TriaAdjTriangles[a]=t;   //the adjacent triangle to the edge a is t
     81                                TriaAdjTriangles[a]=t;    //the adjacent triangle to the edge a is t
    8282                                TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle
    8383                                if(t) { //if t!=NULL add adjacent triangle to t (this)

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

    r3267 r3272  
    2323                PreInit(0);
    2424
     25                /*Read Geometry if provided*/
     26                if(bamggeom->NumEdges>0) {
     27                        Gh.ReadGeometry(bamggeom,bamgopts);
     28                        Gh.AfterRead();
     29                }
     30
    2531                /*Read background mesh*/
    2632                ReadMesh(bamgmesh,bamgopts);
    2733
    28                 /*Read Geometry*/
     34                /*Build Geometry if not provided*/
    2935                if(bamggeom->NumEdges==0) {
    3036                        /*Recreate geometry if needed*/
    3137                        printf("WARNING: mesh present but no geometry found. Reconstructing...\n");
    3238                        BuildGeometryFromMesh(bamgopts);
    33                         Gh.AfterRead();
    34                 }
    35                 else{
    36                         Gh.ReadGeometry(bamggeom,bamgopts);
    3739                        Gh.AfterRead();
    3840                }
    … …  
    511513                        i2=bamgmesh->NumEdgesOnGeometricEdge;
    512514                        for (i1=0;i1<i2;i1++) {
    513                                 i=(int)bamgmesh->EdgesOnGeometricEdge[i*2+0];
    514                                 j=(int)bamgmesh->EdgesOnGeometricEdge[i*2+1];
     515                                printf("%i %i\n",(int)bamgmesh->EdgesOnGeometricEdge[i1*2+0],(int)bamgmesh->EdgesOnGeometricEdge[i1*2+1]);
     516                        }
     517                        for (i1=0;i1<i2;i1++) {
     518                                i=(int)bamgmesh->EdgesOnGeometricEdge[i1*2+0];
     519                                j=(int)bamgmesh->EdgesOnGeometricEdge[i1*2+1];
    515520                                if(!(i>0 && j >0 && i <= nbe && j <= Gh.nbe)) {
    516                                         throw ErrorException(__FUNCT__,exprintf("EdgesOnGeometricEdge error: We must have : (i>0 && j >0 && i <= nbe && j <= Gh.nbe)"));
     521                                        throw ErrorException(__FUNCT__,exprintf("ReadMesh error: EdgesOnGeometricEdge edge provided (line %i: [%i %i]) is incorrect (must be positive, [0<i<nbe=%i 0<j<Gh.nbe=%i]",i1+1,i,j,nbe,Gh.nbe));
    517522                                }
    518523                                edges[i-1].onGeometry = Gh.edges + j-1;
    … …  
    26282633                 */
    26292634
     2635                /*Intermediary*/
    26302636                int verbosity=0;
    26312637
    … …  
    26342640                // find extrema coordinates of vertices pmin,pmax
    26352641                long i;
    2636                 if(verbosity>2) printf("      Completing Triangles structure for a mesh of %i vertices\n",nbv);
     2642                if(verbosity>2) printf("      Reconstruct mesh of %i vertices\n",nbv);
    26372643
    26382644                //initialize ordre
    … …  
    26402646                for (i=0;i<nbv;i++) ordre[i]=0;
    26412647
     2648                //Initialize NbSubDomains
    26422649                NbSubDomains =0;
    26432650
    2644                 /* generation of the adjacence of the triangles*/
    2645 
     2651                /* generation of triangles adjacency*/
     2652
     2653                //First add existing edges
     2654                long kk =0;
    26462655                SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
    2647 
    2648                 //initialize st
     2656                for (i=0;i<nbe;i++){
     2657                        kk=kk+(i==edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
     2658                }
     2659                if (kk != nbe){
     2660                        throw ErrorException(__FUNCT__,exprintf("There are %i double edges in the mesh",kk-nbe));
     2661                }
     2662
     2663                //Add edges of all triangles in existing mesh
    26492664                long* st = new long[nbt*3];
    26502665                for (i=0;i<nbt*3;i++) st[i]=-1;
    2651 
    2652                 //check number of edges
    2653                 long kk =0;
    2654                 for (i=0;i<nbe;i++){
    2655                         kk=kk+(i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));
    2656                 }
    2657                 if (kk != nbe) {
    2658                         throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i",kk-nbe));
    2659                 }
    2660 
    2661                 //
    26622666                for (i=0;i<nbt;i++){
    2663                         for (int j=0;j<3;j++) {
    2664                                 long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),
    2665                                                         Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
    2666                                 long invisible = triangles[i].Hidden(j);
    2667                                 if(st[k]==-1){
    2668                                         st[k]=3*i+j;
    2669                                 }
     2667                        for (int j=0;j<3;j++){
     2668
     2669                                //Add current triangle edge to edge4
     2670                                long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]));
     2671
     2672                                long invisible=triangles[i].Hidden(j);
     2673
     2674                                //If the edge has not been added to st, add it
     2675                                if(st[k]==-1) st[k]=3*i+j;
     2676                               
     2677                                //If the edge already exists, add adjacency
    26702678                                else if(st[k]>=0) {
    26712679                                        assert(!triangles[i].TriangleAdj(j));
    26722680                                        assert(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3)));
    26732681
    2674                                         triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
    2675                                         if (invisible)  triangles[i].SetHidden(j);
    2676                                         if (k<nbe) {
    2677                                                 triangles[i].SetLocked(j);
    2678                                         }
     2682                                        triangles[i].SetAdj2(j,triangles+st[k]/3,(int)(st[k]%3));
     2683                                        if (invisible) triangles[i].SetHidden(j);
     2684                                        if (k<nbe)     triangles[i].SetLocked(j);
     2685
     2686                                        //Make st[k] negative so that it will throw an error message if it is found again
    26792687                                        st[k]=-2-st[k];
    26802688                                }
     2689
     2690                                //An edge belongs to 2 triangles
    26812691                                else {
    26822692                                        throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles",
    … …  
    26852695                        }
    26862696                }
     2697
     2698                //Display info if required
    26872699                if(verbosity>5) {
    2688                         printf("         info on Mesh:\n");
     2700                        printf("         info of Mesh:\n");
    26892701                        printf("            - number of vertices    = %i \n",nbv);
    26902702                        printf("            - number of triangles   = %i \n",nbt);
    … …  
    26942706                }
    26952707
    2696                 // check the consistant of edge[].adj and the geometrical required  vertex
     2708                //check the consistant of edge[].adj and the geometrical required  vertex
    26972709                long k=0;
    26982710                for (i=0;i<edge4->nb();i++){
    2699                         if (st[i] >=0){ // edge alone
     2711                        if (st[i]>=0){ // edge alone
    27002712                                if (i<nbe){
    27012713                                        long i0=edge4->i(i);
    … …  
    27212733
    27222734                /* mesh generation with boundary points*/
    2723                 long nbvb = 0;
     2735                long nbvb=0;
    27242736                for (i=0;i<nbv;i++){
    27252737                        vertices[i].t=0;
    27262738                        vertices[i].vint=0;
    2727                         if (ordre[i]){
    2728                                 ordre[nbvb++]=ordre[i];
    2729                         }
     2739                        if (ordre[i]) ordre[nbvb++]=ordre[i];
    27302740                }
    27312741
    … …  
    27362746                subdomains=0;
    27372747
    2738                 long  Nbtriafillhole = 2*nbvb;
    2739                 Triangle* triafillhole =new Triangle[Nbtriafillhole];
    2740                 triangles =  triafillhole;
     2748                long  Nbtriafillhole=2*nbvb;
     2749                Triangle* triafillhole=new Triangle[Nbtriafillhole];
     2750                triangles = triafillhole;
    27412751
    27422752                nbt=2;
    27432753                nbtx= Nbtriafillhole;
    27442754
    2745                 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;)
     2755                //Find a vertex that is not aligned with vertices 0 and 1
     2756                for (i=2;det(ordre[0]->i,ordre[1]->i,ordre[i]->i)==0;)
    27462757                 if  (++i>=nbvb) {
    27472758                         throw ErrorException(__FUNCT__,exprintf("ReconstructExistingMesh: All the vertices are aligned"));
    27482759                 }
     2760                //Move this vertex (i) to the 2d position in ordre
    27492761                Exchange(ordre[2], ordre[i]);
    27502762
     2763                /*Reconstruct mesh beginning with 2 triangles*/
    27512764                Vertex *  v0=ordre[0], *v1=ordre[1];
    27522765
    … …  
    27742787                triangles[1].link=&triangles[0];
    27752788
    2776                 if (!quadtree )
    2777                  delete  quadtree; // ->ReInitialise();
    2778 
     2789                if (!quadtree) delete quadtree; //ReInitialise;
    27792790                quadtree = new QuadTree(this,0);
    27802791                quadtree->Add(*v0);
    27812792                quadtree->Add(*v1);
    27822793
    2783                 // We add the vertices one by one
     2794                // vertices are added one by one
    27842795                long NbSwap=0;
    27852796                for (int icount=2; icount<nbvb; icount++) {
    … …  
    27922803                }
    27932804
    2794                 // inforce the boundary
     2805                //enforce the boundary
    27952806                TriangleAdjacent ta(0,0);
    27962807                long nbloss = 0,knbe=0;
    27972808                for ( i = 0; i < nbe; i++){
    2798                         if (st[i] >=0){  // edge alone => on border ...  FH oct 2009
    2799                                 Vertex & a=edges[i][0], & b =    edges[i][1];
    2800                                 if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1.
    2801                                   {
     2809                        if (st[i] >=0){ //edge alone => on border
     2810                                Vertex &a=edges[i][0], &b=edges[i][1];
     2811                                if (a.t && b.t){
    28022812                                        knbe++;
    2803                                         if (ForceEdge(a,b,ta)<0)
    2804                                          nbloss++;
    2805                                   }
     2813                                        if (ForceEdge(a,b,ta)<0) nbloss++;
     2814                                }
    28062815                        }
    28072816                }
    28082817                if(nbloss) {
    2809                         throw ErrorException(__FUNCT__,exprintf("we lost(?) %i edges other %i",nbloss,knbe));
     2818                        throw ErrorException(__FUNCT__,exprintf("we lost %i existing edges other %i",nbloss,knbe));
    28102819                }
    28112820
    … …  
    28172826                        if (triangles[i].link){ // remove triangles
    28182827                                krm++;
    2819                                 for (int j=0;j<3;j++)
    2820                                   {
     2828                                for (int j=0;j<3;j++){
    28212829                                        TriangleAdjacent ta =  triangles[i].Adj(j);
    2822                                         Triangle & tta = * (Triangle *) ta;
    2823                                         if(! tta.link) // edge between remove and not remove
    2824                                           { // change the link of ta;
     2830                                        Triangle &tta = *(Triangle*)ta;
     2831                                        //if edge between remove and not remove
     2832                                        if(! tta.link){
     2833                                                // change the link of ta;
    28252834                                                int ja = ta;
    28262835                                                Vertex *v0= ta.EdgeVertex(0);
    28272836                                                Vertex *v1= ta.EdgeVertex(1);
    28282837                                                long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);
    2829                                                 if (st[k]<0){
    2830                                                         throw ErrorException(__FUNCT__,exprintf("st[k]<0"));
    2831                                                 }
     2838
     2839                                                assert(st[k]>=0);
    28322840                                                tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
    28332841                                                ta.SetLock();
    28342842                                                st[k]=-2-st[k];
    2835                                           }
    2836                                   }
     2843                                        }
     2844                                }
    28372845                        }
    28382846                }
    … …  
    28432851                                triangles[i].color=-1;
    28442852                        }
    2845                         else
    2846                           {
     2853                        else{
    28472854                                triangles[i].color= savenbt+ NbTfillHoll++;
    2848                           }
    2849                 }
    2850                 if (savenbt+NbTfillHoll>savenbtx ){
    2851                         throw ErrorException(__FUNCT__,exprintf("savenbt+NbTfillHoll>savenbtx"));
    2852                 }
     2855                        }
     2856                }
     2857                assert(savenbt+NbTfillHoll<=savenbtx);
     2858
    28532859                // copy of the outside triangles in saveTriangles
    28542860                for (i=0;i<nbt;i++){
    … …  
    28622868                k =0;
    28632869                Triangle * tmax = triangles + nbt;
    2864                 for (i=0;i<savenbt;i++) 
    2865                   {
     2870                for (i=0;i<savenbt;i++) {
    28662871                        Triangle & ti = savetriangles[i];
    2867                         for (int j=0;j<3;j++)
    2868                           {
     2872                        for (int j=0;j<3;j++){
    28692873                                Triangle * ta = ti.TriangleAdj(j);
    28702874                                int aa = ti.NuEdgeTriangleAdj(j);
    28712875                                int lck = ti.Locked(j);
    28722876                                if (!ta) k++; // bug
    2873                                 else if ( ta >= triangles && ta < tmax)
    2874                                   {
     2877                                else if ( ta >= triangles && ta < tmax){
    28752878                                        ta= savetriangles + ta->color;
    28762879                                        ti.SetAdj2(j,ta,aa);
    28772880                                        if(lck) ti.SetLocked(j);
    2878                                   }
    2879                           }
    2880                   }
    2881                 //       OutSidesTriangles = triangles;
    2882                 //      long NbOutSidesTriangles = nbt;
     2881                                }
     2882                        }
     2883                }
    28832884
    28842885                // restore triangles;
    … …  
    28962897                delete edge4;
    28972898                delete [] st;
    2898                 for (i=0;i<nbv;i++)
    2899                  quadtree->Add(vertices[i]);
     2899                for (i=0;i<nbv;i++) quadtree->Add(vertices[i]);
    29002900
    29012901                SetVertexFieldOn();
    … …  
    29062906                        if(edges[i].onGeometry){
    29072907                                for(int j=0;j<2;j++){
     2908                                        /*Go through the edges adjacent to current edge (if on the same curve)*/
    29082909                                        if (!edges[i].adj[j]){
    2909                                                 if(!edges[i][j].onGeometry->IsRequiredVertex()) {
    2910                                                         printf("Error: adj and vertex requires edges [%i %i]=%i on = %i\n",i,j,Number(edges[i][j]),Gh.Number(edges[i].onGeometry));
     2910                                                /*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/
     2911                                                /*Check that the 2 vertices are on geometry AND required*/
     2912                                                if(!edges[i][j].onGeometry->IsRequiredVertex()){
     2913                                                        printf("ReconstructExistingMesh error message: problem with the edge %i: [%i][%i]\n",Number(edges[i][j])+1,i,j);
     2914                                                        printf("   This edge is on geometry (Edge of geometry %i)",Gh.Number(edges[i].onGeometry));
    29112915                                                        if (edges[i][j].onGeometry->OnGeomVertex())
    29122916                                                         printf(" Vertex %i\n",Gh.Number(edges[i][j].onGeometry->gv));
    29132917                                                        else if (edges[i][j].onGeometry->OnGeomEdge())
    2914                                                          printf(" Edges %i \n",Gh.Number(edges[i][j].onGeometry->ge));
     2918                                                         printf(" Edges %i\n",Gh.Number(edges[i][j].onGeometry->ge));
    29152919                                                        else
    29162920                                                         printf(" = %p\n",edges[i][j].onGeometry);

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

    r3263 r3272  
    3030                        union {
    3131                                Triangle* t; // one triangle which is containing the vertex
    32                                 long color;
    33                                 Vertex* to;  // use in geometry Vertex to now the Mesh Vertex associed
    34                                 VertexOnGeom* onGeometry;       // if vint == 8; // set with Triangles::SetVertexFieldOn()
    35                                 Vertex* onBackgroundVertex;     // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh()
    36                                 VertexOnEdge* onBackgroundEdge; // if vint == 32 on Background edge
     32                                long      color;
     33                                Vertex*   to;  // use in geometry Vertex to now the Mesh Vertex associed
     34                                VertexOnGeom* onGeometry;        // if vint == 8; // set with Triangles::SetVertexFieldOn()
     35                                Vertex*       onBackgroundVertex;// if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh()
     36                                VertexOnEdge* onBackgroundEdge;  // if vint == 32 on Background edge
    3737                        };
    3838
    3939                        //Operators
    40                         operator  I2() const {return i;}             // Cast operator
    41                         operator  const R2 & () const {return r;}    // Cast operator
    42                         operator Metric () const {return m;}         // Cast operator
     40                        operator I2() const {return i;}             // Cast operator
     41                        operator const R2 & () const {return r;}    // Cast operator
     42                        operator Metric () const {return m;}        // Cast operator
    4343                        double operator()(R2 x) const { return m(x);} // Get x in the metric m
    4444
    4545                        //methods (No constructor and no destructors...)
    4646                        double Smoothing(Triangles & ,const Triangles & ,Triangle  * & ,double =1);
    47                         void  MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts);
    48                         void  Echo();
    49                         int   ref() const { return ReferenceNumber;}
    50                         long  Optim(int =1,int =0);
     47                        void   MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts);
     48                        void   Echo();
     49                        int    ref() const { return ReferenceNumber;}
     50                        long   Optim(int =1,int =0);
    5151
    5252                        //inline functions
    … …  
    5454        };
    5555
    56         //FOR NOW (WARNING)
     56        //Intermediary
    5757        double QuadQuality(const Vertex &,const Vertex &,const Vertex &,const Vertex &);
    5858}

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

    r3267 r3272  
    4040
    4141                        //Methods
    42                         int  OnGeomVertex()const {return this? abscisse <0 :0;}
     42                        int  OnGeomVertex()const{return this? abscisse <0 :0;}
    4343                        int  OnGeomEdge() const {return this? abscisse >=0 :0;}
    44                         int  IsRequiredVertex(){ return this? (( abscisse<0 ? (gv?gv->Required():0):0 )) : 0;}
     44                        int  IsRequiredVertex() {return this? ((abscisse<0 ? (gv?gv->Required():0):0 )) : 0;}
    4545                        void SetOn(){mv->onGeometry=this;mv->vint=IsVertexOnGeom;}
    4646

  • issm/trunk/src/m/classes/public/bamg.m

    r3268 r3272  
    139139                                                %Get position of the two nodes of the edge in domain
    140140                                                i1=j;
    141                                                 i2=mod(j+1,domain(1).nods-1);
     141                                                i2=mod(j+1,domain(1).nods);
    142142
    143143                                                %rift is crossing edge [i1 i2] of the domain
    … …  
    168168                                                count=count+nods;
    169169
    170                                                 disp(' done'); break;
     170                                                break;
    171171                                        end
    172172                                end
    … …  
    187187        end
    188188
    189         %update other fields of bamg_geometry
    190         bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1);
    191         bamg_geometry.NumEdges=size(bamg_geometry.Edges,1);
    192         bamg_geometry.NumSubDomains=size(bamg_geometry.SubDomains,1);
     189else
     190
     191        %Use segments as the geometry
     192        %bamg_geometry.Vertices=[md.x md.y ones(md.numberofgrids,1)];
     193        %bamg_geometry.Edges=[md.segments(:,1:2) md.segmentmarkers];
     194        %bamg_geometry.Edges
     195
    193196end
     197
     198%update other fields of bamg_geometry
     199bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1);
     200bamg_geometry.NumEdges=size(bamg_geometry.Edges,1);
     201bamg_geometry.NumSubDomains=size(bamg_geometry.SubDomains,1);
    194202%}}}
    195203
    … …  
    201209bamg_mesh.NumEdges=0;
    202210bamg_mesh.Edges=zeros(0,3);
     211bamg_mesh.NumEdgesOnGeometricEdge=0;
     212bamg_mesh.EdgesOnGeometricEdge=zeros(0,2);
     213bamg_mesh.NumVerticesOnGeometricEdge=0;
     214bamg_mesh.VerticesOnGeometricEdge=zeros(0,2);
    203215bamg_mesh.hVertices=[];
    204216bamg_mesh.NumCrackedEdges=0;
    … …  
    215227                end
    216228        end
     229
    217230        if isstruct(md.rifts)
    218                 for i=1:length(md.rifts)
    219                         bamg_mesh.Edges=[bamg_mesh.Edges;md.rifts(i).segments];
    220                 end
    221                 bamg_mesh.NumEdges=size(bamg_mesh.Edges,1);
    222                 bamg_mesh.NumCrackedEdges=bamg_mesh.NumEdges;
    223                 bamg_mesh.CrackedEdges=[1 2;3 4];
     231                error('bamg error message: rifts not supported yet. Do meshprocessrift AFTER bamg');
    224232        end
    225233end
    … …  
    265273
    266274%call Bamg
    267 [triangles vertices segments segmentsmarkers metric]=Bamg(bamg_mesh,bamg_geometry,bamg_options);
     275[bamgmesh_out bamggeom_out]=Bamg(bamg_mesh,bamg_geometry,bamg_options);
    268276
    269277% plug results onto model
    270 md.x=vertices(:,1);
    271 md.y=vertices(:,2);
    272 md.elements=triangles(:,1:3);
    273 md.segments=segments;
    274 md.segmentmarkers=segmentsmarkers;
    275 md.dummy=metric;
     278md.x=bamgmesh_out.Vertices(:,1);
     279md.y=bamgmesh_out.Vertices(:,2);
     280md.elements=bamgmesh_out.Triangles(:,1:3);
     281md.segments=bamgmesh_out.Segments;
     282md.segmentmarkers=bamgmesh_out.SegmentsMarkers;
    276283
    277284%Fill in rest of fields:

  • issm/trunk/src/mex/Bamg/Bamg.cpp

    r3271 r3272  
    4242        FetchData(&bamgmesh_in.NumCrackedEdges,mxGetField(BAMGMESH,0,"NumCrackedEdges"));
    4343        FetchData(&bamgmesh_in.CrackedEdges,NULL,NULL,mxGetField(BAMGMESH,0,"CrackedEdges"));
     44        FetchData(&bamgmesh_in.NumEdges,mxGetField(BAMGMESH,0,"NumEdges"));
     45        FetchData(&bamgmesh_in.Edges,NULL,NULL,mxGetField(BAMGMESH,0,"Edges"));
     46        FetchData(&bamgmesh_in.NumEdgesOnGeometricEdge,mxGetField(BAMGMESH,0,"NumEdgesOnGeometricEdge"));
     47        FetchData(&bamgmesh_in.EdgesOnGeometricEdge,NULL,NULL,mxGetField(BAMGMESH,0,"EdgesOnGeometricEdge"));
     48        FetchData(&bamgmesh_in.NumVerticesOnGeometricEdge,mxGetField(BAMGMESH,0,"NumVerticesOnGeometricEdge"));
     49        FetchData(&bamgmesh_in.VerticesOnGeometricEdge,NULL,NULL,mxGetField(BAMGMESH,0,"VerticesOnGeometricEdge"));
    4450        if (bamgmesh_in.hVertices && (cols!=1 || lines!=bamgmesh_in.NumVertices)){throw ErrorException(__FUNCT__,exprintf("the size of 'hVertices' should be [%i %i]",bamgmesh_in.NumVertices,1));}
    4551
    … …  
    7783        Bamgx(&bamgmesh_out,&bamggeom_out,&bamgmesh_in,&bamggeom_in,&bamgopts);
    7884
    79         printf("ok0\n");
    80 
    8185        /*Variables*/
    8286        mxArray*    bamgmesh_mat=NULL;
    … …  
    114118        bamgmesh_mat=mxCreateStructArray(ndim,dimensions,numfields,fnames);
    115119
    116         printf("ok2\n");
    117          
    118120        mxSetField(bamgmesh_mat,0,"NumTriangles",mxCreateDoubleScalar(bamgmesh_out.NumTriangles));
    119121
    … …  
    222224        mxSetField(bamgmesh_mat,0,"SubDomainsFromGeom",pfield2);
    223225
    224         printf("ok3\n");
    225 
    226 
    227226        /*assign output datasets: */
    228         printf("ok4\n");
    229227        *BAMGMESHOUT=bamgmesh_mat;
    230         printf("ok5\n");
    231228
    232229        /*Variables*/
    … …  
    270267        mxSetField(bamggeom_mat,0,"Edges",pfield2);
    271268
    272 
    273269        mxSetField(bamggeom_mat,0,"NumTangentAtEdges",mxCreateDoubleScalar(bamggeom_out.NumTangentAtEdges));
    274270
    … …  
    317313
    318314        /*assign output datasets: */
    319         printf("ok4\n");
    320315        *BAMGGEOMOUT=bamggeom_mat;
    321         printf("ok5\n");
    322316
    323317        /*Free ressources: */
Note: See TracChangeset for help on using the changeset viewer.