Index: /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.h =================================================================== --- /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.h (revision 3271) +++ /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.h (revision 3272) @@ -35,25 +35,25 @@ //Methods - R2 F(double theta) const ; // parametrization of the curve edge + R2 F(double theta) const ; // parametrization of the curve edge double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente - int Cracked() const {return flag & 1;} - int Dup() const { return flag & 32;} - int Equi()const {return flag & 2;} - int ReverseEqui()const {return flag & 128;} - int TgA()const {return flag &4;} - int TgB()const {return flag &8;} - int Tg(int i) const{return i==0 ? TgA() : TgB();} - int Mark()const {return flag &16;} - int Required() { return flag & 64;} - void SetCracked() { flag |= 1;} - void SetDup() { flag |= 32;} // not a real edge - void SetEqui() { flag |= 2;} - void SetTgA() { flag|=4;} - void SetTgB() { flag|=8;} - void SetMark() { flag|=16;} - void SetUnMark() { flag &= 1007 /* 1023-16*/;} - void SetRequired() { flag |= 64;} - void SetReverseEqui() {flag |= 128;} - void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew); + int Tg(int i) const{return i==0 ? TgA() : TgB();} + int Cracked() const {return flag & 1; } + int Dup() const {return flag & 32; } + int Equi()const {return flag & 2; } + int ReverseEqui()const {return flag & 128;} + int TgA()const {return flag & 4; } + int TgB()const {return flag & 8; } + int Mark()const {return flag & 16; } + int Required() {return flag & 64; } + void SetCracked() { flag |= 1; } + void SetDup() { flag |= 32; } // not a real edge + void SetEqui() { flag |= 2; } + void SetTgA() { flag |=4; } + void SetTgB() { flag |=8; } + void SetMark() { flag |=16; } + void SetUnMark() { flag &= 1007 /* 1023-16*/;} + void SetRequired() { flag |= 64; } + void SetReverseEqui() { flag |= 128;} + void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew); }; Index: /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp (revision 3271) +++ /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp (revision 3272) @@ -1,3 +1,4 @@ +#include #include "BamgObjects.h" @@ -18,5 +19,5 @@ //initialize fields nx =nnx; //number of vertices - nbax =mmx; // 3 * number of triangles + nbax =mmx; // 3 * number of triangles NbOfEdges=0; head = new long [nx]; @@ -25,5 +26,5 @@ //initialize head (-1 everywhere) i=nx; - while (i--) head[i]=-1; + while(i--) head[i]=-1; } /*}}}1*/ @@ -66,10 +67,6 @@ int h,n; - //check that head is not empty - if (head==NULL) { - throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add no head is NULL (How come?)")); - } - //get n from h (usually h=ii) + assert(head); n=head[h=Abs(ii)%nx]; Index: /issm/trunk/src/c/Bamgx/objects/SetOfE4.h =================================================================== --- /issm/trunk/src/c/Bamgx/objects/SetOfE4.h (revision 3271) +++ /issm/trunk/src/c/Bamgx/objects/SetOfE4.h (revision 3272) @@ -21,9 +21,16 @@ public: + + // Constructors SetOfEdges4(long ,long);// nb Edges mx , nb de sommet ~SetOfEdges4() {delete [] head; delete [] Edges;} + + //operators + IntEdge & operator[](long k){return Edges[k];} + + //Methods long add (long ii,long jj); long SortAndAdd (long ii,long jj) {return ii <=jj ? add (ii,jj) : add (jj,ii) ;} - long nb(){return NbOfEdges;} + long nb(){return NbOfEdges;} long find (long ii,long jj); long SortAndFind (long ii,long jj) {return ii <=jj ? find (ii,jj) : find (jj,ii) ;} @@ -31,7 +38,4 @@ long j(long k){return Edges[k].j;} long newarete(long k){return NbOfEdges == k+1;} - - //operators - IntEdge & operator[](long k){return Edges[k];} }; } Index: /issm/trunk/src/c/Bamgx/objects/Triangle.h =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangle.h (revision 3271) +++ /issm/trunk/src/c/Bamgx/objects/Triangle.h (revision 3272) @@ -79,5 +79,5 @@ } void SetAdj2(short a,Triangle *t,short aat){ - TriaAdjTriangles[a]=t; //the adjacent triangle to the edge a is t + TriaAdjTriangles[a]=t; //the adjacent triangle to the edge a is t TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle if(t) { //if t!=NULL add adjacent triangle to t (this) Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3271) +++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3272) @@ -23,16 +23,18 @@ PreInit(0); + /*Read Geometry if provided*/ + if(bamggeom->NumEdges>0) { + Gh.ReadGeometry(bamggeom,bamgopts); + Gh.AfterRead(); + } + /*Read background mesh*/ ReadMesh(bamgmesh,bamgopts); - /*Read Geometry*/ + /*Build Geometry if not provided*/ if(bamggeom->NumEdges==0) { /*Recreate geometry if needed*/ printf("WARNING: mesh present but no geometry found. Reconstructing...\n"); BuildGeometryFromMesh(bamgopts); - Gh.AfterRead(); - } - else{ - Gh.ReadGeometry(bamggeom,bamgopts); Gh.AfterRead(); } @@ -511,8 +513,11 @@ i2=bamgmesh->NumEdgesOnGeometricEdge; for (i1=0;i1EdgesOnGeometricEdge[i*2+0]; - j=(int)bamgmesh->EdgesOnGeometricEdge[i*2+1]; + printf("%i %i\n",(int)bamgmesh->EdgesOnGeometricEdge[i1*2+0],(int)bamgmesh->EdgesOnGeometricEdge[i1*2+1]); + } + for (i1=0;i1EdgesOnGeometricEdge[i1*2+0]; + j=(int)bamgmesh->EdgesOnGeometricEdge[i1*2+1]; if(!(i>0 && j >0 && i <= nbe && j <= Gh.nbe)) { - throw ErrorException(__FUNCT__,exprintf("EdgesOnGeometricEdge error: We must have : (i>0 && j >0 && i <= nbe && j <= Gh.nbe)")); + throw ErrorException(__FUNCT__,exprintf("ReadMesh error: EdgesOnGeometricEdge edge provided (line %i: [%i %i]) is incorrect (must be positive, [02) printf(" Completing Triangles structure for a mesh of %i vertices\n",nbv); + if(verbosity>2) printf(" Reconstruct mesh of %i vertices\n",nbv); //initialize ordre @@ -2640,43 +2646,47 @@ for (i=0;iSortAndAdd(Number(edges[i][0]),Number(edges[i][1]))); + } + if (kk != nbe){ + throw ErrorException(__FUNCT__,exprintf("There are %i double edges in the mesh",kk-nbe)); + } + + //Add edges of all triangles in existing mesh long* st = new long[nbt*3]; for (i=0;iSortAndAdd(Number(edges[i][0]),Number(edges[i][1]))); - } - if (kk != nbe) { - throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i",kk-nbe)); - } - - // for (i=0;iSortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), - Number(triangles[i][VerticesOfTriangularEdge[j][1]])); - long invisible = triangles[i].Hidden(j); - if(st[k]==-1){ - st[k]=3*i+j; - } + for (int j=0;j<3;j++){ + + //Add current triangle edge to edge4 + long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])); + + long invisible=triangles[i].Hidden(j); + + //If the edge has not been added to st, add it + if(st[k]==-1) st[k]=3*i+j; + + //If the edge already exists, add adjacency else if(st[k]>=0) { assert(!triangles[i].TriangleAdj(j)); assert(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3))); - triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); - if (invisible) triangles[i].SetHidden(j); - if (k5) { - printf(" info on Mesh:\n"); + printf(" info of Mesh:\n"); printf(" - number of vertices = %i \n",nbv); printf(" - number of triangles = %i \n",nbt); @@ -2694,8 +2706,8 @@ } - // check the consistant of edge[].adj and the geometrical required vertex + //check the consistant of edge[].adj and the geometrical required vertex long k=0; for (i=0;inb();i++){ - if (st[i] >=0){ // edge alone + if (st[i]>=0){ // edge alone if (ii(i); @@ -2721,11 +2733,9 @@ /* mesh generation with boundary points*/ - long nbvb = 0; + long nbvb=0; for (i=0;ii, ordre[1]->i, ordre[i]->i ) == 0;) + //Find a vertex that is not aligned with vertices 0 and 1 + for (i=2;det(ordre[0]->i,ordre[1]->i,ordre[i]->i)==0;) if (++i>=nbvb) { throw ErrorException(__FUNCT__,exprintf("ReconstructExistingMesh: All the vertices are aligned")); } + //Move this vertex (i) to the 2d position in ordre Exchange(ordre[2], ordre[i]); + /*Reconstruct mesh beginning with 2 triangles*/ Vertex * v0=ordre[0], *v1=ordre[1]; @@ -2774,12 +2787,10 @@ triangles[1].link=&triangles[0]; - if (!quadtree ) - delete quadtree; // ->ReInitialise(); - + if (!quadtree) delete quadtree; //ReInitialise; quadtree = new QuadTree(this,0); quadtree->Add(*v0); quadtree->Add(*v1); - // We add the vertices one by one + // vertices are added one by one long NbSwap=0; for (int icount=2; icount=0){ // edge alone => on border ... FH oct 2009 - Vertex & a=edges[i][0], & b = edges[i][1]; - if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1. - { + if (st[i] >=0){ //edge alone => on border + Vertex &a=edges[i][0], &b=edges[i][1]; + if (a.t && b.t){ knbe++; - if (ForceEdge(a,b,ta)<0) - nbloss++; - } + if (ForceEdge(a,b,ta)<0) nbloss++; + } } } if(nbloss) { - throw ErrorException(__FUNCT__,exprintf("we lost(?) %i edges other %i",nbloss,knbe)); + throw ErrorException(__FUNCT__,exprintf("we lost %i existing edges other %i",nbloss,knbe)); } @@ -2817,22 +2826,21 @@ if (triangles[i].link){ // remove triangles krm++; - for (int j=0;j<3;j++) - { + for (int j=0;j<3;j++){ TriangleAdjacent ta = triangles[i].Adj(j); - Triangle & tta = * (Triangle *) ta; - if(! tta.link) // edge between remove and not remove - { // change the link of ta; + Triangle &tta = *(Triangle*)ta; + //if edge between remove and not remove + if(! tta.link){ + // change the link of ta; int ja = ta; Vertex *v0= ta.EdgeVertex(0); Vertex *v1= ta.EdgeVertex(1); long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv); - if (st[k]<0){ - throw ErrorException(__FUNCT__,exprintf("st[k]<0")); - } + + assert(st[k]>=0); tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3)); ta.SetLock(); st[k]=-2-st[k]; - } - } + } + } } } @@ -2843,12 +2851,10 @@ triangles[i].color=-1; } - else - { + else{ triangles[i].color= savenbt+ NbTfillHoll++; - } - } - if (savenbt+NbTfillHoll>savenbtx ){ - throw ErrorException(__FUNCT__,exprintf("savenbt+NbTfillHoll>savenbtx")); - } + } + } + assert(savenbt+NbTfillHoll<=savenbtx); + // copy of the outside triangles in saveTriangles for (i=0;i= triangles && ta < tmax) - { + else if ( ta >= triangles && ta < tmax){ ta= savetriangles + ta->color; ti.SetAdj2(j,ta,aa); if(lck) ti.SetLocked(j); - } - } - } - // OutSidesTriangles = triangles; - // long NbOutSidesTriangles = nbt; + } + } + } // restore triangles; @@ -2896,6 +2897,5 @@ delete edge4; delete [] st; - for (i=0;iAdd(vertices[i]); + for (i=0;iAdd(vertices[i]); SetVertexFieldOn(); @@ -2906,11 +2906,15 @@ if(edges[i].onGeometry){ for(int j=0;j<2;j++){ + /*Go through the edges adjacent to current edge (if on the same curve)*/ if (!edges[i].adj[j]){ - if(!edges[i][j].onGeometry->IsRequiredVertex()) { - printf("Error: adj and vertex requires edges [%i %i]=%i on = %i\n",i,j,Number(edges[i][j]),Gh.Number(edges[i].onGeometry)); + /*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/ + /*Check that the 2 vertices are on geometry AND required*/ + if(!edges[i][j].onGeometry->IsRequiredVertex()){ + printf("ReconstructExistingMesh error message: problem with the edge %i: [%i][%i]\n",Number(edges[i][j])+1,i,j); + printf(" This edge is on geometry (Edge of geometry %i)",Gh.Number(edges[i].onGeometry)); if (edges[i][j].onGeometry->OnGeomVertex()) printf(" Vertex %i\n",Gh.Number(edges[i][j].onGeometry->gv)); else if (edges[i][j].onGeometry->OnGeomEdge()) - printf(" Edges %i \n",Gh.Number(edges[i][j].onGeometry->ge)); + printf(" Edges %i\n",Gh.Number(edges[i][j].onGeometry->ge)); else printf(" = %p\n",edges[i][j].onGeometry); Index: /issm/trunk/src/c/Bamgx/objects/Vertex.h =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Vertex.h (revision 3271) +++ /issm/trunk/src/c/Bamgx/objects/Vertex.h (revision 3272) @@ -30,23 +30,23 @@ union { Triangle* t; // one triangle which is containing the vertex - long color; - Vertex* to; // use in geometry Vertex to now the Mesh Vertex associed - VertexOnGeom* onGeometry; // if vint == 8; // set with Triangles::SetVertexFieldOn() - Vertex* onBackgroundVertex; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh() - VertexOnEdge* onBackgroundEdge; // if vint == 32 on Background edge + long color; + Vertex* to; // use in geometry Vertex to now the Mesh Vertex associed + VertexOnGeom* onGeometry; // if vint == 8; // set with Triangles::SetVertexFieldOn() + Vertex* onBackgroundVertex;// if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh() + VertexOnEdge* onBackgroundEdge; // if vint == 32 on Background edge }; //Operators - operator I2() const {return i;} // Cast operator - operator const R2 & () const {return r;} // Cast operator - operator Metric () const {return m;} // Cast operator + operator I2() const {return i;} // Cast operator + operator const R2 & () const {return r;} // Cast operator + operator Metric () const {return m;} // Cast operator double operator()(R2 x) const { return m(x);} // Get x in the metric m //methods (No constructor and no destructors...) double Smoothing(Triangles & ,const Triangles & ,Triangle * & ,double =1); - 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); - void Echo(); - int ref() const { return ReferenceNumber;} - long Optim(int =1,int =0); + 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); + void Echo(); + int ref() const { return ReferenceNumber;} + long Optim(int =1,int =0); //inline functions @@ -54,5 +54,5 @@ }; - //FOR NOW (WARNING) + //Intermediary double QuadQuality(const Vertex &,const Vertex &,const Vertex &,const Vertex &); } Index: /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h =================================================================== --- /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h (revision 3271) +++ /issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h (revision 3272) @@ -40,7 +40,7 @@ //Methods - int OnGeomVertex()const {return this? abscisse <0 :0;} + int OnGeomVertex()const{return this? abscisse <0 :0;} int OnGeomEdge() const {return this? abscisse >=0 :0;} - int IsRequiredVertex(){ return this? (( abscisse<0 ? (gv?gv->Required():0):0 )) : 0;} + int IsRequiredVertex() {return this? ((abscisse<0 ? (gv?gv->Required():0):0 )) : 0;} void SetOn(){mv->onGeometry=this;mv->vint=IsVertexOnGeom;} Index: /issm/trunk/src/m/classes/public/bamg.m =================================================================== --- /issm/trunk/src/m/classes/public/bamg.m (revision 3271) +++ /issm/trunk/src/m/classes/public/bamg.m (revision 3272) @@ -139,5 +139,5 @@ %Get position of the two nodes of the edge in domain i1=j; - i2=mod(j+1,domain(1).nods-1); + i2=mod(j+1,domain(1).nods); %rift is crossing edge [i1 i2] of the domain @@ -168,5 +168,5 @@ count=count+nods; - disp(' done'); break; + break; end end @@ -187,9 +187,17 @@ end - %update other fields of bamg_geometry - bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1); - bamg_geometry.NumEdges=size(bamg_geometry.Edges,1); - bamg_geometry.NumSubDomains=size(bamg_geometry.SubDomains,1); +else + + %Use segments as the geometry + %bamg_geometry.Vertices=[md.x md.y ones(md.numberofgrids,1)]; + %bamg_geometry.Edges=[md.segments(:,1:2) md.segmentmarkers]; + %bamg_geometry.Edges + end + +%update other fields of bamg_geometry +bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1); +bamg_geometry.NumEdges=size(bamg_geometry.Edges,1); +bamg_geometry.NumSubDomains=size(bamg_geometry.SubDomains,1); %}}} @@ -201,4 +209,8 @@ bamg_mesh.NumEdges=0; bamg_mesh.Edges=zeros(0,3); +bamg_mesh.NumEdgesOnGeometricEdge=0; +bamg_mesh.EdgesOnGeometricEdge=zeros(0,2); +bamg_mesh.NumVerticesOnGeometricEdge=0; +bamg_mesh.VerticesOnGeometricEdge=zeros(0,2); bamg_mesh.hVertices=[]; bamg_mesh.NumCrackedEdges=0; @@ -215,11 +227,7 @@ end end + if isstruct(md.rifts) - for i=1:length(md.rifts) - bamg_mesh.Edges=[bamg_mesh.Edges;md.rifts(i).segments]; - end - bamg_mesh.NumEdges=size(bamg_mesh.Edges,1); - bamg_mesh.NumCrackedEdges=bamg_mesh.NumEdges; - bamg_mesh.CrackedEdges=[1 2;3 4]; + error('bamg error message: rifts not supported yet. Do meshprocessrift AFTER bamg'); end end @@ -265,13 +273,12 @@ %call Bamg -[triangles vertices segments segmentsmarkers metric]=Bamg(bamg_mesh,bamg_geometry,bamg_options); +[bamgmesh_out bamggeom_out]=Bamg(bamg_mesh,bamg_geometry,bamg_options); % plug results onto model -md.x=vertices(:,1); -md.y=vertices(:,2); -md.elements=triangles(:,1:3); -md.segments=segments; -md.segmentmarkers=segmentsmarkers; -md.dummy=metric; +md.x=bamgmesh_out.Vertices(:,1); +md.y=bamgmesh_out.Vertices(:,2); +md.elements=bamgmesh_out.Triangles(:,1:3); +md.segments=bamgmesh_out.Segments; +md.segmentmarkers=bamgmesh_out.SegmentsMarkers; %Fill in rest of fields: Index: /issm/trunk/src/mex/Bamg/Bamg.cpp =================================================================== --- /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 3271) +++ /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 3272) @@ -42,4 +42,10 @@ FetchData(&bamgmesh_in.NumCrackedEdges,mxGetField(BAMGMESH,0,"NumCrackedEdges")); FetchData(&bamgmesh_in.CrackedEdges,NULL,NULL,mxGetField(BAMGMESH,0,"CrackedEdges")); + FetchData(&bamgmesh_in.NumEdges,mxGetField(BAMGMESH,0,"NumEdges")); + FetchData(&bamgmesh_in.Edges,NULL,NULL,mxGetField(BAMGMESH,0,"Edges")); + FetchData(&bamgmesh_in.NumEdgesOnGeometricEdge,mxGetField(BAMGMESH,0,"NumEdgesOnGeometricEdge")); + FetchData(&bamgmesh_in.EdgesOnGeometricEdge,NULL,NULL,mxGetField(BAMGMESH,0,"EdgesOnGeometricEdge")); + FetchData(&bamgmesh_in.NumVerticesOnGeometricEdge,mxGetField(BAMGMESH,0,"NumVerticesOnGeometricEdge")); + FetchData(&bamgmesh_in.VerticesOnGeometricEdge,NULL,NULL,mxGetField(BAMGMESH,0,"VerticesOnGeometricEdge")); 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));} @@ -77,6 +83,4 @@ Bamgx(&bamgmesh_out,&bamggeom_out,&bamgmesh_in,&bamggeom_in,&bamgopts); - printf("ok0\n"); - /*Variables*/ mxArray* bamgmesh_mat=NULL; @@ -114,6 +118,4 @@ bamgmesh_mat=mxCreateStructArray(ndim,dimensions,numfields,fnames); - printf("ok2\n"); - mxSetField(bamgmesh_mat,0,"NumTriangles",mxCreateDoubleScalar(bamgmesh_out.NumTriangles)); @@ -222,11 +224,6 @@ mxSetField(bamgmesh_mat,0,"SubDomainsFromGeom",pfield2); - printf("ok3\n"); - - /*assign output datasets: */ - printf("ok4\n"); *BAMGMESHOUT=bamgmesh_mat; - printf("ok5\n"); /*Variables*/ @@ -270,5 +267,4 @@ mxSetField(bamggeom_mat,0,"Edges",pfield2); - mxSetField(bamggeom_mat,0,"NumTangentAtEdges",mxCreateDoubleScalar(bamggeom_out.NumTangentAtEdges)); @@ -317,7 +313,5 @@ /*assign output datasets: */ - printf("ok4\n"); *BAMGGEOMOUT=bamggeom_mat; - printf("ok5\n"); /*Free ressources: */