Changeset 5148
- Timestamp:
- 08/11/10 10:38:41 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 22 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/modules/Bamgx/Bamgx.cpp
r5146 r5148 184 184 //display info 185 185 if(verbosity>0) { 186 if (Th.nbt-Th.nbtout-Th. NbOfQuad*2){187 printf(" new number of triangles = %i\n",(Th.nbt-Th.nbtout-Th. NbOfQuad*2));188 } 189 if (Th. NbOfQuad){190 printf(" new number of quads = %i\n",Th. NbOfQuad);186 if (Th.nbt-Th.nbtout-Th.nbq*2){ 187 printf(" new number of triangles = %i\n",(Th.nbt-Th.nbtout-Th.nbq*2)); 188 } 189 if (Th.nbq ){ 190 printf(" new number of quads = %i\n",Th.nbq); 191 191 } 192 192 } -
issm/trunk/src/c/objects/Bamg/AdjacentTriangle.cpp
r5143 r5148 16 16 /*FUNCTION AdjacentTriangle::Locked {{{1*/ 17 17 int AdjacentTriangle::Locked() const { 18 return t-> TriaAdjSharedEdge[a] & 4;18 return t->AdjEdgeNumber[a] & 4; 19 19 } 20 20 /*}}}*/ 21 21 /*FUNCTION AdjacentTriangle::MarkUnSwap {{{1*/ 22 22 int AdjacentTriangle::MarkUnSwap() const { 23 return t-> TriaAdjSharedEdge[a] & 8;23 return t->AdjEdgeNumber[a] & 8; 24 24 } 25 25 /*}}}*/ … … 27 27 int AdjacentTriangle::GetAllFlag_UnSwap() const { 28 28 // take all flag except MarkUnSwap 29 return t-> TriaAdjSharedEdge[a] & 1012;29 return t->AdjEdgeNumber[a] & 1012; 30 30 } 31 31 /*}}}*/ … … 42 42 /*FUNCTION AdjacentTriangle::EdgeVertex {{{1*/ 43 43 BamgVertex* AdjacentTriangle::EdgeVertex(const int & i) const { 44 return t-> TriaVertices[VerticesOfTriangularEdge[a][i]];44 return t->vertices[VerticesOfTriangularEdge[a][i]]; 45 45 } 46 46 /*}}}*/ 47 47 /*FUNCTION AdjacentTriangle::OppositeVertex {{{1*/ 48 48 BamgVertex* AdjacentTriangle::OppositeVertex() const { 49 return t-> TriaVertices[bamg::OppositeVertex[a]];49 return t->vertices[bamg::OppositeVertex[a]]; 50 50 } 51 51 /*}}}*/ … … 64 64 //set Adjacent Triangle of a triangle 65 65 if(t) { 66 t-> TriaAdjTriangles[a]=ta.t;67 t-> TriaAdjSharedEdge[a]=ta.a|l;66 t->adj[a]=ta.t; 67 t->AdjEdgeNumber[a]=ta.a|l; 68 68 } 69 69 if(ta.t) { 70 ta.t-> TriaAdjTriangles[ta.a] = t ;71 ta.t-> TriaAdjSharedEdge[ta.a] = a| l ;70 ta.t->adj[ta.a] = t ; 71 ta.t->AdjEdgeNumber[ta.a] = a| l ; 72 72 } 73 73 } -
issm/trunk/src/c/objects/Bamg/BamgVertex.h
r5130 r5148 28 28 29 29 union { 30 Triangle *t; // one triangle which is containing the vertex30 Triangle *t; // one triangle which is containing the vertex 31 31 long color; 32 BamgVertex *to; // used in geometry BamgVertex to know the Mesh BamgVertex associated33 VertexOnGeom * onGeometry;// if vint == 8; // set with Mesh::SetVertexFieldOn()34 BamgVertex * onBackgroundVertex; // if vint == 16 on Background vertex Mesh::SetVertexFieldOnBTh()35 VertexOnEdge * onBackgroundEdge; // if vint == 32 on Background edge32 BamgVertex *to; // used in geometry BamgVertex to know the Mesh Vertex associated 33 VertexOnGeom *GeometricalEdgeHook; // if vint == 8; // set with Mesh::SetVertexFieldOn() 34 BamgVertex *BackgroundVertexHook; // if vint == 16 on Background vertex Mesh::SetVertexFieldOnBTh() 35 VertexOnEdge *BackgroundEdgeHook; // if vint == 32 on Background edge 36 36 }; 37 37 … … 46 46 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); 47 47 void Echo(); 48 int ref() const { return ReferenceNumber;}48 int GetReferenceNumber() const { return ReferenceNumber;} 49 49 long Optim(int =1,int =0); 50 50 -
issm/trunk/src/c/objects/Bamg/Direction.cpp
r5130 r5148 23 23 24 24 /*Methods*/ 25 /*FUNCTION Direction:: sens{{{1*/26 int Direction:: sens(Icoor1 i,Icoor1 j) {25 /*FUNCTION Direction::direction{{{1*/ 26 int Direction::direction(Icoor1 i,Icoor1 j) { 27 27 int r =1; 28 28 if (dir!= MaxICoor) { -
issm/trunk/src/c/objects/Bamg/Direction.h
r5130 r5148 15 15 Direction(); 16 16 Direction(Icoor1 i,Icoor1 j); 17 int sens(Icoor1 i,Icoor1 j);17 int direction(Icoor1 i,Icoor1 j); 18 18 }; 19 19 } -
issm/trunk/src/c/objects/Bamg/Edge.cpp
r5130 r5148 18 18 v[0] = ThNew.vertices + Th.Number(v[0]); 19 19 v[1] = ThNew.vertices + Th.Number(v[1]); 20 if ( onGeometry)21 onGeometry = ThNew.Gh.edges+Th.Gh.Number(onGeometry);20 if (GeometricalEdgeHook) 21 GeometricalEdgeHook = ThNew.Gh.edges+Th.Gh.Number(GeometricalEdgeHook); 22 22 if (adj[0]) adj[0] = ThNew.edges + Th.Number(adj[0]); 23 23 if (adj[1]) adj[1] = ThNew.edges + Th.Number(adj[1]); … … 28 28 printf("Edge:\n"); 29 29 printf(" pointers towards two vertices: %p %p\n",v[0],v[1]); 30 printf(" ref = %i\n",ref);31 printf(" onGeometry = %p\n",onGeometry);30 printf(" ReferenceNumber = %i\n",ReferenceNumber); 31 printf(" GeometricalEdgeHook = %p\n",GeometricalEdgeHook); 32 32 printf(" two adjacent edges on the same curve: %p %p\n",adj[0],adj[1]); 33 33 } -
issm/trunk/src/c/objects/Bamg/Edge.h
r5143 r5148 17 17 public: 18 18 BamgVertex *v[2]; 19 long ref;20 GeometricalEdge * onGeometry;19 long ReferenceNumber; 20 GeometricalEdge *GeometricalEdgeHook; 21 21 Edge *adj[2]; // the 2 adj edges if on the same curve 22 22 -
issm/trunk/src/c/objects/Bamg/GeometricalEdge.h
r5146 r5148 14 14 public: 15 15 GeometricalVertex *v[2]; 16 long ref;16 long ReferenceNumber; 17 17 long CurveNumber; 18 18 R2 tg[2]; // the 2 tangentes (tg[0] =0 => no continuity) … … 29 29 R2 F(double theta) const ; // parametrization of the curve edge 30 30 double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente 31 int Tg(int i) const {return i==0 ? TgA() : TgB();}32 int Cracked() const {return type & 1;}33 int TgA() const{return type & 4;}34 int TgB() const {return type & 8;}35 int Mark() const {return type & 16;}36 int Required() {return type & 64;}31 int Tg(int i) const {return i==0 ? TgA() : TgB(); } 32 int Cracked() const {return type &1; } 33 int TgA() const {return type &4; } 34 int TgB() const {return type &8; } 35 int Mark() const {return type &16; } 36 int Required() {return type &64; } 37 37 void SetCracked() { type |= 1;} 38 38 void SetTgA() { type |=4; } -
issm/trunk/src/c/objects/Bamg/GeometricalSubDomain.h
r3913 r5148 13 13 public: 14 14 GeometricalEdge *edge; 15 int sens;// -1 or 116 long ref;15 int direction; // -1 or 1 16 long ReferenceNumber; 17 17 18 18 //Methods -
issm/trunk/src/c/objects/Bamg/Geometry.cpp
r5146 r5148 35 35 vertices = nbv ? new GeometricalVertex[nbv] : NULL; 36 36 edges = nbe ? new GeometricalEdge[nbe]:NULL; 37 curves= NbOfCurves ? new Curve[NbOfCurves]:NULL;38 subdomains = NbSubDomains ? new GeometricalSubDomain[NbSubDomains]:NULL;37 curves= nbcurves ? new Curve[nbcurves]:NULL; 38 subdomains = nbsubdomains ? new GeometricalSubDomain[nbsubdomains]:NULL; 39 39 for (i=0;i<nbe;i++) 40 40 edges[i].Set(Gh.edges[i],Gh,*this); 41 for (i=0;i< NbOfCurves;i++)41 for (i=0;i<nbcurves;i++) 42 42 curves[i].Set(Gh.curves[i],Gh,*this); 43 for (i=0;i< NbSubDomains;i++)43 for (i=0;i<nbsubdomains;i++) 44 44 subdomains[i].Set(Gh.subdomains[i],Gh,*this); 45 45 } … … 51 51 if(vertices) delete [] vertices; vertices=0; 52 52 if(edges) delete [] edges; edges=0; 53 if(triangles) delete [] triangles; triangles=0;54 53 if(quadtree) delete quadtree; quadtree=0; 55 if(curves) delete []curves; curves=0; NbOfCurves=0;54 if(curves) delete []curves; curves=0;nbcurves=0; 56 55 if(subdomains) delete [] subdomains;subdomains=0; 57 56 Init(); … … 66 65 nbv=0; 67 66 nbe=0; 68 NbOfCurves=0;67 nbcurves=0; 69 68 70 69 double Hmin = HUGE_VAL;// the infinie value … … 140 139 edges[i].v[0]= vertices + i1; //pointer toward vertex i1 (=&vertices[i1]) 141 140 edges[i].v[1]= vertices + i2; //pointer toward vertex i2 142 edges[i]. ref=(long)bamggeom->Edges[i*3+2];141 edges[i].ReferenceNumber=(long)bamggeom->Edges[i*3+2]; 143 142 144 143 //get length of edge … … 155 154 156 155 //Cracked? 157 if (edges[i]. ref!=1) edges[i].SetCracked();156 if (edges[i].ReferenceNumber!=1) edges[i].SetCracked(); 158 157 159 158 //prepare metric … … 263 262 if(verbose>5) printf(" processing SubDomains\n"); 264 263 if (bamggeom->SubDomainsSize[1]!=4) ISSMERROR("SubDomains should have 4 columns"); 265 NbSubDomains=bamggeom->SubDomainsSize[0];266 subdomains = new GeometricalSubDomain[ NbSubDomains];267 for (i=0;i< NbSubDomains;i++){264 nbsubdomains=bamggeom->SubDomainsSize[0]; 265 subdomains = new GeometricalSubDomain[nbsubdomains]; 266 for (i=0;i<nbsubdomains;i++){ 268 267 i0=(int)bamggeom->SubDomains[i*4+0]; 269 268 i1=(int)bamggeom->SubDomains[i*4+1]; … … 273 272 if (i1>nbe || i1<=0) ISSMERROR("Bad Subdomain definition: second number should in [1 %i] (edge number)",nbe); 274 273 subdomains[i].edge=edges + (i1-1); 275 subdomains[i]. sens= (int) i2;276 subdomains[i]. ref= i3;274 subdomains[i].direction = (int) i2; 275 subdomains[i].ReferenceNumber = i3; 277 276 } 278 277 } … … 304 303 bamggeom->Vertices[i*3+0]=vertices[i].r.x; 305 304 bamggeom->Vertices[i*3+1]=vertices[i].r.y; 306 bamggeom->Vertices[i*3+2]=vertices[i]. ref();305 bamggeom->Vertices[i*3+2]=vertices[i].GetReferenceNumber(); 307 306 308 307 //update counters … … 320 319 bamggeom->Edges[i*3+0]=Number(edges[i][0])+1; //back to Matlab indexing 321 320 bamggeom->Edges[i*3+1]=Number(edges[i][1])+1; //back to Matlab indexing 322 bamggeom->Edges[i*3+2]=(double)edges[i]. ref;321 bamggeom->Edges[i*3+2]=(double)edges[i].ReferenceNumber; 323 322 324 323 //update counters … … 363 362 /*SubDomains*/ 364 363 if(verbose>5) printf(" writing SubDomains\n"); 365 bamggeom->SubDomainsSize[0]= NbSubDomains;364 bamggeom->SubDomainsSize[0]=nbsubdomains; 366 365 bamggeom->SubDomainsSize[1]=4; 367 if ( NbSubDomains){368 bamggeom->SubDomains=(double*)xmalloc(4* NbSubDomains*sizeof(double));369 for (i=0;i< NbSubDomains;i++){366 if (nbsubdomains){ 367 bamggeom->SubDomains=(double*)xmalloc(4*nbsubdomains*sizeof(double)); 368 for (i=0;i<nbsubdomains;i++){ 370 369 bamggeom->SubDomains[4*i+0]=2; 371 370 bamggeom->SubDomains[4*i+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing 372 bamggeom->SubDomains[4*i+2]=subdomains[i]. sens;373 bamggeom->SubDomains[4*i+3]=subdomains[i]. ref;371 bamggeom->SubDomains[4*i+2]=subdomains[i].direction; 372 bamggeom->SubDomains[4*i+3]=subdomains[i].ReferenceNumber; 374 373 } 375 374 } … … 585 584 vertices[i].SetCorner() ; 586 585 } 587 // if the refa changing then is SetRequired();586 // if the ReferenceNumber a changing then is SetRequired(); 588 587 if (edges[i1].type != edges[i2].type || edges[i1].Required()){ 589 588 vertices[i].SetRequired(); 590 589 } 591 if (edges[i1]. ref != edges[i2].ref) {590 if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) { 592 591 vertices[i].SetRequired(); 593 592 } … … 661 660 for (int step=0;step<2;step++){ 662 661 for (i=0;i<nbe;i++) edges[i].SetUnMark(); 663 NbOfCurves = 0;662 nbcurves = 0; 664 663 long nbgem=0; 665 664 for (int level=0;level < 2 && nbgem != nbe;level++) … … 673 672 GeometricalVertex *a=(*e)(k0); // begin 674 673 if(curves){ 675 curves[ NbOfCurves].be=e;676 curves[ NbOfCurves].kb=k0;674 curves[nbcurves].be=e; 675 curves[nbcurves].kb=k0; 677 676 } 678 677 int nee=0; … … 682 681 e->SetMark(); 683 682 nbgem++; 684 e->CurveNumber= NbOfCurves;683 e->CurveNumber=nbcurves; 685 684 if(curves) { 686 curves[ NbOfCurves].ee=e;687 curves[ NbOfCurves].ke=k1;685 curves[nbcurves].ee=e; 686 curves[nbcurves].ke=k1; 688 687 } 689 688 … … 693 692 e = e->Adj[k1]; // next edge 694 693 } 695 NbOfCurves++;694 nbcurves++; 696 695 if(level) a->SetRequired(); 697 696 } … … 700 699 ISSMASSERT(nbgem && nbe); 701 700 if(step==0){ 702 curves = new Curve[ NbOfCurves];701 curves = new Curve[nbcurves]; 703 702 } 704 703 } 705 for(int i=0;i< NbOfCurves ;i++){704 for(int i=0;i<nbcurves ;i++){ 706 705 GeometricalEdge * be=curves[i].be, *eqbe=be; 707 706 //GeometricalEdge * ee=curves[i].ee, *eqee=be; … … 749 748 printf(" nbv (number of vertices) : %i\n",nbv); 750 749 printf(" nbe (number of edges) : %i\n",nbe); 751 printf(" NbSubDomains: %i\n",NbSubDomains);752 printf(" NbOfCurves: %i\n",NbOfCurves);750 printf(" nbsubdomains: %i\n",nbsubdomains); 751 printf(" nbcurves: %i\n",nbcurves); 753 752 printf(" vertices: %p\n",vertices); 754 printf(" triangles: %p\n",triangles);755 753 printf(" edges: %p\n",edges); 756 754 printf(" quadtree: %p\n",quadtree); … … 774 772 quadtree=0; 775 773 curves=0; 776 triangles=0;777 774 edges=0; 778 775 vertices=0; 779 NbSubDomains=0;780 NbOfCurves=0;776 nbsubdomains=0; 777 nbcurves=0; 781 778 subdomains=0; 782 779 MaxCornerAngle = 10*Pi/180; //default is 10 degres … … 822 819 823 820 s=save_s; 824 GeometricalEdge* on=e. onGeometry;821 GeometricalEdge* on=e.GeometricalEdgeHook; 825 822 if (!on){ 826 823 ISSMERROR("ProjectOnCurve error message: edge provided should be on geometry"); 827 824 } 828 if (!e[0]. onGeometry || !e[1].onGeometry){825 if (!e[0].GeometricalEdgeHook || !e[1].GeometricalEdgeHook){ 829 826 ISSMERROR("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry"); 830 827 } … … 838 835 839 836 //Get geometrical vertices corresponding to v0 and v1 840 VertexOnGeom vg0=*v0. onGeometry, vg1=*v1.onGeometry;837 VertexOnGeom vg0=*v0.GeometricalEdgeHook, vg1=*v1.GeometricalEdgeHook; 841 838 842 839 //build two pointers towrad current geometrical edge … … 846 843 R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; 847 844 int OppositeSens = (V01,AB)<0; 848 int sens0=0,sens1=1;845 int direction0=0,direction1=1; 849 846 if (OppositeSens) s=1-s,Exchange(vg0,vg1),Exchange(V0,V1); 850 847 … … 854 851 const int mxe=100; 855 852 GeometricalEdge* ge[mxe+1]; 856 int sensge[mxe+1];853 int directionge[mxe+1]; 857 854 double lge[mxe+1]; 858 855 int bge=mxe/2,tge=bge; 859 ge[bge] = e. onGeometry;860 sensge[bge]=1;861 862 while (eg0!=(GeometricalEdge*)vg0 && (*eg0)( sens0)!=(GeometricalVertex*)vg0){856 ge[bge] = e.GeometricalEdgeHook; 857 directionge[bge]=1; 858 859 while (eg0!=(GeometricalEdge*)vg0 && (*eg0)(direction0)!=(GeometricalVertex*)vg0){ 863 860 if (bge<=0) { 864 861 if(NbTry) { 865 printf("Fatal Error: on the class trianglesbefore call Geometry::ProjectOnCurve\n");862 printf("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve\n"); 866 863 printf("That bug might come from:\n"); 867 864 printf(" 1) a mesh edge containing more than %i geometrical edges\n",mxe/2); … … 874 871 } 875 872 GeometricalEdge* tmpge = eg0; 876 ge[--bge] =eg0 = eg0->Adj[ sens0];873 ge[--bge] =eg0 = eg0->Adj[direction0]; 877 874 ISSMASSERT(bge>=0 && bge<=mxe); 878 sens0 = 1-( sensge[bge] = tmpge->DirAdj[sens0]);879 } 880 while (eg1 != (GeometricalEdge*) vg1 && (*eg1)( sens1) != (GeometricalVertex*) vg1) {875 direction0 = 1-( directionge[bge] = tmpge->DirAdj[direction0]); 876 } 877 while (eg1 != (GeometricalEdge*) vg1 && (*eg1)(direction1) != (GeometricalVertex*) vg1) { 881 878 if(tge>=mxe ) { 882 printf("WARNING: on the class trianglesbefore call Geometry::ProjectOnCurve is having issues (isn't it Eric?)\n");879 printf("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)\n"); 883 880 NbTry++; 884 881 if (NbTry<2) goto retry; 885 printf("Fatal Error: on the class trianglesbefore call Geometry::ProjectOnCurve\n");882 printf("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve\n"); 886 883 printf("That bug might come from:\n"); 887 884 printf(" 1) a mesh edge contening more than %i geometrical edges\n",mxe/2); … … 891 888 } 892 889 GeometricalEdge* tmpge = eg1; 893 ge[++tge] =eg1 = eg1->Adj[ sens1];894 sensge[tge]= sens1 = 1-tmpge->DirAdj[sens1];890 ge[++tge] =eg1 = eg1->Adj[direction1]; 891 directionge[tge]= direction1 = 1-tmpge->DirAdj[direction1]; 895 892 ISSMASSERT(tge>=0 && tge<=mxe); 896 893 } 897 894 898 895 899 if ((*eg0)( sens0)==(GeometricalVertex*)vg0)900 vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0, sens0); //vg0 = absisce901 902 if ((*eg1)( sens1)==(GeometricalVertex*)vg1)903 vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1, sens1);896 if ((*eg0)(direction0)==(GeometricalVertex*)vg0) 897 vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,direction0); //vg0 = absisce 898 899 if ((*eg1)(direction1)==(GeometricalVertex*)vg1) 900 vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,direction1); 904 901 905 902 double sg; … … 916 913 for(i=bge;i<tge;i++){ 917 914 ISSMASSERT(i>=0 && i<=mxe); 918 BB = (*ge[i])[ sensge[i]];915 BB = (*ge[i])[directionge[i]]; 919 916 lge[i]=ll += Norme2(AA-BB); 920 917 AA=BB ;} … … 925 922 on =0; 926 923 s0 = vg0; 927 s1= sensge[bge];924 s1= directionge[bge]; 928 925 double l0=0,l1; 929 926 i=bge; 930 927 while ( (l1=lge[i]) < ls ) { 931 928 ISSMASSERT(i>=0 && i<=mxe); 932 i++,s0=1-(s1= sensge[i]),l0=l1;929 i++,s0=1-(s1=directionge[i]),l0=l1; 933 930 } 934 931 on=ge[i]; -
issm/trunk/src/c/objects/Bamg/Geometry.h
r5130 r5148 23 23 long nbv; // number of vertices 24 24 long nbe; // number of edges 25 long NbSubDomains;26 long NbOfCurves;25 long nbsubdomains; 26 long nbcurves; 27 27 GeometricalVertex *vertices; 28 Triangle *triangles;29 28 GeometricalEdge *edges; 30 29 QuadTree *quadtree; -
issm/trunk/src/c/objects/Bamg/Mesh.cpp
r5146 r5148 106 106 { 107 107 vertices[nbv] = Tho.vertices[i]; 108 if (!vertices[nbv]. ref())108 if (!vertices[nbv].GetReferenceNumber()) 109 109 vertices[nbv].ReferenceNumber = refv[i]; 110 110 nbv++; … … 127 127 } 128 128 triangles[nbt] = Triangle(this,kk[i0],kk[i1],kk[i2]); 129 triangles[nbt].color = Tho.subdomains[reft[i]]. ref;129 triangles[nbt].color = Tho.subdomains[reft[i]].ReferenceNumber; 130 130 nbt++; 131 131 } … … 145 145 ReconstructExistingMesh(); 146 146 147 if (! NbSubDomains){148 ISSMERROR(" NbSubDomains==0");147 if (!nbsubdomains){ 148 ISSMERROR("nbsubdomains==0"); 149 149 } 150 150 if (!subdomains[0].head || !subdomains[0].head->link){ … … 168 168 nbt = Th.nbt; 169 169 nbe = Th.nbe; 170 NbSubDomains = Th.NbSubDomains;170 nbsubdomains = Th.nbsubdomains; 171 171 nbtout = Th.nbtout; 172 NbOfQuad = Th.NbOfQuad ; 173 NbOfSwapTriangle =0; 172 nbq = Th.nbq ; 174 173 NbVerticesOnGeomVertex = Th.NbVerticesOnGeomVertex; 175 174 if(NbVerticesOnGeomVertex) … … 198 197 if(nbe) 199 198 edges = new Edge[nbe]; 200 if( NbSubDomains)201 subdomains = new SubDomain[ NbSubDomains];199 if(nbsubdomains) 200 subdomains = new SubDomain[nbsubdomains]; 202 201 pmin = Th.pmin; 203 202 pmax = Th.pmax; … … 209 208 for(i=0;i<nbv;i++) 210 209 vertices[i].Set(Th.vertices[i],Th,*this); 211 for(i=0;i< NbSubDomains;i++)210 for(i=0;i<nbsubdomains;i++) 212 211 subdomains[i].Set(Th,i,*this); 213 212 for (i=0;i<NbVerticesOnGeomVertex;i++) … … 424 423 i1=(int)bamgmesh->Edges[i*3+0]-1; //-1 for C indexing 425 424 i2=(int)bamgmesh->Edges[i*3+1]-1; //-1 for C indexing 426 edges[i]. ref=(long)bamgmesh->Edges[i*3+2];425 edges[i].ReferenceNumber=(long)bamgmesh->Edges[i*3+2]; 427 426 edges[i].v[0]= vertices +i1; 428 427 edges[i].v[1]= vertices +i2; … … 484 483 ISSMERROR("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+1,j+1,nbe,Gh.nbe); 485 484 } 486 edges[i]. onGeometry=Gh.edges+j;485 edges[i].GeometricalEdgeHook=Gh.edges+j; 487 486 } 488 487 } … … 500 499 //SubDomain 501 500 if(bamgmesh->SubDomains){ 502 long i3,head, sens;501 long i3,head,direction; 503 502 if(verbose>5) printf(" processing SubDomains\n"); 504 NbSubDomains=bamgmesh->SubDomainsSize[0];505 subdomains = new SubDomain [ NbSubDomains ];506 for (i=0;i< NbSubDomains;i++) {503 nbsubdomains=bamgmesh->SubDomainsSize[0]; 504 subdomains = new SubDomain [ nbsubdomains ]; 505 for (i=0;i<nbsubdomains;i++) { 507 506 i3 =(int)bamgmesh->SubDomains[i*3+0]; 508 507 head=(int)bamgmesh->SubDomains[i*3+1]-1;//C indexing 509 sens=(int)bamgmesh->SubDomains[i*3+2];508 direction=(int)bamgmesh->SubDomains[i*3+2]; 510 509 if (i3!=23) ISSMERROR("Bad Subdomain definition: first number should be 3"); 511 510 if (head<0 || head>=nbt) ISSMERROR("Bad Subdomain definition: head should in [1 %i] (triangle number)",nbt); … … 581 580 bamgmesh->Vertices[i*3+0]=vertices[i].r.x; 582 581 bamgmesh->Vertices[i*3+1]=vertices[i].r.y; 583 bamgmesh->Vertices[i*3+2]=vertices[i]. ref();582 bamgmesh->Vertices[i*3+2]=vertices[i].GetReferenceNumber(); 584 583 } 585 584 } … … 595 594 bamgmesh->Edges[i*3+0]=Number(edges[i][0])+1; //back to M indexing 596 595 bamgmesh->Edges[i*3+1]=Number(edges[i][1])+1; //back to M indexing 597 bamgmesh->Edges[i*3+2]=edges[i]. ref;598 if(edges[i]. onGeometry){596 bamgmesh->Edges[i*3+2]=edges[i].ReferenceNumber; 597 if(edges[i].GeometricalEdgeHook){ 599 598 NumIssmSegments++; 600 599 } … … 665 664 num=0; 666 665 for (i=0;i<nbe;i++){ 667 if(edges[i]. onGeometry){666 if(edges[i].GeometricalEdgeHook){ 668 667 //build segment 669 668 int i1=Number(edges[i][0]); … … 677 676 bamgmesh->IssmSegments[num*4+1]=Number(edges[i][1])+1; //back to M indexing 678 677 bamgmesh->IssmSegments[num*4+2]=(int)j/3+1; //back to M indexing 679 bamgmesh->IssmSegments[num*4+3]=edges[i]. ref;678 bamgmesh->IssmSegments[num*4+3]=edges[i].ReferenceNumber; 680 679 num+=1; 681 680 stop=true; … … 686 685 bamgmesh->IssmSegments[num*4+1]=Number(edges[i][0])+1; //back to M indexing 687 686 bamgmesh->IssmSegments[num*4+2]=(int)j/3+1; //back to M indexing 688 bamgmesh->IssmSegments[num*4+3]=edges[i]. ref;687 bamgmesh->IssmSegments[num*4+3]=edges[i].ReferenceNumber; 689 688 num+=1; 690 689 stop=true; … … 703 702 /*Triangles*/ 704 703 if(verbose>5) printf(" writing Triangles\n"); 705 k=nbInT- NbOfQuad*2;704 k=nbInT-nbq*2; 706 705 num=0; 707 706 bamgmesh->TrianglesSize[0]=k; … … 716 715 bamgmesh->Triangles[num*4+1]=Number(t[1])+1; //back to M indexing 717 716 bamgmesh->Triangles[num*4+2]=Number(t[2])+1; //back to M indexing 718 bamgmesh->Triangles[num*4+3]=subdomains[reft[i]]. ref;717 bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ReferenceNumber; 719 718 num=num+1; 720 719 } … … 724 723 /*Quadrilaterals*/ 725 724 if(verbose>5) printf(" writing Quadrilaterals\n"); 726 bamgmesh->QuadrilateralsSize[0]= NbOfQuad;725 bamgmesh->QuadrilateralsSize[0]=nbq; 727 726 bamgmesh->QuadrilateralsSize[1]=5; 728 if ( NbOfQuad){729 bamgmesh->Quadrilaterals=(double*)xmalloc(5* NbOfQuad*sizeof(double));727 if (nbq){ 728 bamgmesh->Quadrilaterals=(double*)xmalloc(5*nbq*sizeof(double)); 730 729 for (i=0;i<nbt;i++){ 731 730 Triangle &t =triangles[i]; … … 738 737 bamgmesh->Quadrilaterals[i*5+2]=Number(v2)+1; //back to M indexing 739 738 bamgmesh->Quadrilaterals[i*5+3]=Number(v3)+1; //back to M indexing 740 bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]]. ref;739 bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]].ReferenceNumber; 741 740 } 742 741 } … … 745 744 /*SubDomains*/ 746 745 if(verbose>5) printf(" writing SubDomains\n"); 747 bamgmesh->SubDomainsSize[0]= NbSubDomains;746 bamgmesh->SubDomainsSize[0]=nbsubdomains; 748 747 bamgmesh->SubDomainsSize[1]=4; 749 if ( NbSubDomains){750 bamgmesh->SubDomains=(double*)xmalloc(4* NbSubDomains*sizeof(double));751 for (i=0;i< NbSubDomains;i++){748 if (nbsubdomains){ 749 bamgmesh->SubDomains=(double*)xmalloc(4*nbsubdomains*sizeof(double)); 750 for (i=0;i<nbsubdomains;i++){ 752 751 bamgmesh->SubDomains[i*4+0]=3; 753 752 bamgmesh->SubDomains[i*4+1]=reft[Number(subdomains[i].head)]; 754 753 bamgmesh->SubDomains[i*4+2]=1; 755 bamgmesh->SubDomains[i*4+3]=subdomains[i]. ref;754 bamgmesh->SubDomains[i*4+3]=subdomains[i].ReferenceNumber; 756 755 } 757 756 } … … 759 758 /*SubDomainsFromGeom*/ 760 759 if(verbose>5) printf(" writing SubDomainsFromGeom\n"); 761 bamgmesh->SubDomainsFromGeomSize[0]=Gh. NbSubDomains;760 bamgmesh->SubDomainsFromGeomSize[0]=Gh.nbsubdomains; 762 761 bamgmesh->SubDomainsFromGeomSize[1]=4; 763 if (Gh. NbSubDomains){764 bamgmesh->SubDomainsFromGeom=(double*)xmalloc(4*Gh. NbSubDomains*sizeof(double));765 for (i=0;i<Gh. NbSubDomains;i++){762 if (Gh.nbsubdomains){ 763 bamgmesh->SubDomainsFromGeom=(double*)xmalloc(4*Gh.nbsubdomains*sizeof(double)); 764 for (i=0;i<Gh.nbsubdomains;i++){ 766 765 bamgmesh->SubDomainsFromGeom[i*4+0]=2; 767 766 bamgmesh->SubDomainsFromGeom[i*4+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing 768 bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i]. sens;769 bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i]. ref;767 bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].direction; 768 bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ReferenceNumber; 770 769 } 771 770 } … … 806 805 k=0; 807 806 for (i=0;i<nbe;i++){ 808 if (edges[i]. onGeometry) k=k+1;807 if (edges[i].GeometricalEdgeHook) k=k+1; 809 808 } 810 809 bamgmesh->EdgesOnGeometricEdgeSize[0]=k; … … 814 813 int count=0; 815 814 for (i=0;i<nbe;i++){ 816 if (edges[i]. onGeometry){815 if (edges[i].GeometricalEdgeHook){ 817 816 bamgmesh->EdgesOnGeometricEdge[count*2+0]=(double)i+1; //back to Matlab indexing 818 bamgmesh->EdgesOnGeometricEdge[count*2+1]=(double)Gh.Number(edges[i]. onGeometry)+1; //back to Matlab indexing817 bamgmesh->EdgesOnGeometricEdge[count*2+1]=(double)Gh.Number(edges[i].GeometricalEdgeHook)+1; //back to Matlab indexing 819 818 count=count+1; 820 819 } … … 1400 1399 edges[add].v[0] = &triangles[it][VerticesOfTriangularEdge[j][0]]; 1401 1400 edges[add].v[1] = &triangles[it][VerticesOfTriangularEdge[j][1]]; 1402 edges[add]. onGeometry=NULL;1401 edges[add].GeometricalEdgeHook=NULL; 1403 1402 //if already existed 1404 1403 if (i<nbeold){ 1405 edges[add]. ref=edgessave[i].ref;1406 edges[add]. onGeometry=edgessave[i].onGeometry; // HACK to get required edges1404 edges[add].ReferenceNumber=edgessave[i].ReferenceNumber; 1405 edges[add].GeometricalEdgeHook=edgessave[i].GeometricalEdgeHook; // HACK to get required edges 1407 1406 printf("oh no...\n"); 1408 1407 } 1409 1408 else 1410 edges[add]. ref=Min(edges[add].v[0]->ref(),edges[add].v[1]->ref());1409 edges[add].ReferenceNumber=Min(edges[add].v[0]->GetReferenceNumber(),edges[add].v[1]->GetReferenceNumber()); 1411 1410 } 1412 1411 } … … 1474 1473 /*Reconstruct subdomains info*/ 1475 1474 1476 //check that NbSubDomains is empty1477 if ( NbSubDomains){1478 ISSMERROR(" NbSubDomains should be 0");1479 } 1480 NbSubDomains=0;1475 //check that nbsubdomains is empty 1476 if (nbsubdomains){ 1477 ISSMERROR("nbsubdomains should be 0"); 1478 } 1479 nbsubdomains=0; 1481 1480 1482 1481 //color the subdomains … … 1494 1493 1495 1494 //color = number of subdomains 1496 colorT[it]= NbSubDomains;1495 colorT[it]=nbsubdomains; 1497 1496 1498 1497 //color all the adjacent triangles of T that share a non marked edge … … 1509 1508 //color the adjacent triangle 1510 1509 if ( ! t->Locked(j) && tt && (colorT[jt = Number(tt)] == -1) && ( tt->color==kolor)) { 1511 colorT[jt]= NbSubDomains;1510 colorT[jt]=nbsubdomains; 1512 1511 st[++level]=jt; 1513 1512 st[++level]=0; … … 1517 1516 else level-=2; 1518 1517 } 1519 NbSubDomains++;1520 } 1521 } 1522 if (verbose> 3) printf(" The Number of sub domain = %i\n", NbSubDomains);1518 nbsubdomains++; 1519 } 1520 } 1521 if (verbose> 3) printf(" The Number of sub domain = %i\n",nbsubdomains); 1523 1522 1524 1523 //build subdomains 1525 1524 long isd; 1526 subdomains = new SubDomain[ NbSubDomains];1525 subdomains = new SubDomain[nbsubdomains]; 1527 1526 1528 1527 //initialize subdomains[isd].head as 0 1529 for (isd=0;isd< NbSubDomains;isd++) subdomains[isd].head =0;1528 for (isd=0;isd<nbsubdomains;isd++) subdomains[isd].head =0; 1530 1529 1531 1530 k=0; … … 1535 1534 if ((!tt || tt->color != triangles[it].color) && !subdomains[isd=colorT[it]].head){ 1536 1535 subdomains[isd].head = triangles+it; 1537 subdomains[isd]. ref= triangles[it].color;1538 subdomains[isd]. sens= j; // hack1536 subdomains[isd].ReferenceNumber = triangles[it].color; 1537 subdomains[isd].direction = j; // hack 1539 1538 subdomains[isd].edge = 0; 1540 1539 k++; … … 1543 1542 } 1544 1543 //check that we have been through all subdomains 1545 if (k!= NbSubDomains){1546 ISSMERROR("k!= NbSubDomains");1544 if (k!= nbsubdomains){ 1545 ISSMERROR("k!= nbsubdomains"); 1547 1546 } 1548 1547 //delete colorT and st … … 1569 1568 Gh.vertices = new GeometricalVertex[k]; 1570 1569 Gh.edges = new GeometricalEdge[nbe]; 1571 Gh. NbSubDomains = NbSubDomains;1572 Gh.subdomains = new GeometricalSubDomain[ NbSubDomains];1570 Gh.nbsubdomains = nbsubdomains; 1571 Gh.subdomains = new GeometricalSubDomain[nbsubdomains]; 1573 1572 if (verbose>3) printf(" number of vertices = %i\n number of edges = %i\n",Gh.nbv,Gh.nbe); 1574 1573 NbVerticesOnGeomVertex = Gh.nbv; … … 1632 1631 Gh.edges[i].tg[1]=R2(); 1633 1632 1634 bool required= edges[i]. onGeometry;1633 bool required= edges[i].GeometricalEdgeHook; 1635 1634 if(required) kreq++; 1636 edges[i]. onGeometry= Gh.edges + i;1635 edges[i].GeometricalEdgeHook = Gh.edges + i; 1637 1636 if(required){ 1638 1637 Gh.edges[i].v[0]->SetRequired(); … … 1651 1650 len[j1] += l12; 1652 1651 hmin = Min(hmin,l12); 1653 Gh.edges[i]. ref = edges[i].ref;1652 Gh.edges[i].ReferenceNumber = edges[i].ReferenceNumber; 1654 1653 1655 1654 k = edge4->SortAndAdd(i0,i1); … … 1670 1669 1671 1670 //Build Gh.subdomains 1672 for (i=0;i< NbSubDomains;i++){1671 for (i=0;i<nbsubdomains;i++){ 1673 1672 it = Number(subdomains[i].head); 1674 j = subdomains[i]. sens;1673 j = subdomains[i].direction; 1675 1674 long i0 = Number(triangles[it][VerticesOfTriangularEdge[j][0]]); 1676 1675 long i1 = Number(triangles[it][VerticesOfTriangularEdge[j][1]]); 1677 1676 k = edge4->SortAndFind(i0,i1); 1678 1677 if(k>=0){ 1679 subdomains[i]. sens= (vertices + i0 == edges[k].v[0]) ? 1 : -1;1678 subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1; 1680 1679 subdomains[i].edge = edges+k; 1681 1680 Gh.subdomains[i].edge = Gh.edges + k; 1682 Gh.subdomains[i]. sens = subdomains[i].sens;1683 Gh.subdomains[i]. ref = subdomains[i].ref;1681 Gh.subdomains[i].direction = subdomains[i].direction; 1682 Gh.subdomains[i].ReferenceNumber = subdomains[i].ReferenceNumber; 1684 1683 } 1685 1684 else … … 2203 2202 // computed the number of cracked edge 2204 2203 for (k=i=0;i<nbe;i++){ 2205 if(edges[i]. onGeometry->Cracked()) k++;2204 if(edges[i].GeometricalEdgeHook->Cracked()) k++; 2206 2205 } 2207 2206 … … 2222 2221 2223 2222 for (i=0;i<nbe;i++){ 2224 if(edges[i]. onGeometry->Cracked()){2223 if(edges[i].GeometricalEdgeHook->Cracked()){ 2225 2224 2226 2225 //Fill edges fields of CrackedEdges 2227 CrackedEdges[k ].E =edges[i]. onGeometry;2226 CrackedEdges[k ].E =edges[i].GeometricalEdgeHook; 2228 2227 CrackedEdges[k++].e1=&edges[i]; 2229 2228 … … 2459 2458 2460 2459 2461 if (OutSide|| !Gh.subdomains || !Gh. NbSubDomains )2460 if (OutSide|| !Gh.subdomains || !Gh.nbsubdomains ) 2462 2461 { // No geom sub domain 2463 2462 long i; 2464 2463 if (subdomains) delete [] subdomains; 2465 2464 subdomains = new SubDomain[ NbSubDomTot]; 2466 NbSubDomains= NbSubDomTot;2467 for ( i=0;i< NbSubDomains;i++) {2465 nbsubdomains= NbSubDomTot; 2466 for ( i=0;i<nbsubdomains;i++) { 2468 2467 subdomains[i].head=NULL; 2469 subdomains[i]. ref=i+1;2468 subdomains[i].ReferenceNumber=i+1; 2470 2469 } 2471 2470 long * mark = new long[nbt]; … … 2488 2487 // else if(mark[it] == -2 ) triangles[it].Draw(999); 2489 2488 it++;} // end white (it<nbt) 2490 if (k!= NbSubDomains){2491 ISSMERROR("k!= NbSubDomains");2489 if (k!=nbsubdomains){ 2490 ISSMERROR("k!=nbsubdomains"); 2492 2491 } 2493 2492 if(OutSide) … … 2496 2495 // because in this case we have only the true boundary edge 2497 2496 // so teh boundary is manifold 2498 long nbk = NbSubDomains;2497 long nbk = nbsubdomains; 2499 2498 while (nbk) 2500 2499 for (it=0;it<nbt && nbk ;it++) … … 2505 2504 long kr = mark[it]; 2506 2505 if(kr !=kl) { 2507 if (kl >=0 && subdomains[kl]. ref <0 && kr >=0 && subdomains[kr].ref>=0)2508 nbk--,subdomains[kr]. ref=subdomains[kl].ref-1;2509 if (kr >=0 && subdomains[kr]. ref <0 && kl >=0 && subdomains[kl].ref>=0)2510 nbk--,subdomains[kl]. ref=subdomains[kr].ref-1;2511 if(kr<0 && kl >=0 && subdomains[kl]. ref>=0)2512 nbk--,subdomains[kl]. ref=-1;2513 if(kl<0 && kr >=0 && subdomains[kr]. ref>=0)2514 nbk--,subdomains[kr]. ref=-1;2506 if (kl >=0 && subdomains[kl].ReferenceNumber <0 && kr >=0 && subdomains[kr].ReferenceNumber>=0) 2507 nbk--,subdomains[kr].ReferenceNumber=subdomains[kl].ReferenceNumber-1; 2508 if (kr >=0 && subdomains[kr].ReferenceNumber <0 && kl >=0 && subdomains[kl].ReferenceNumber>=0) 2509 nbk--,subdomains[kl].ReferenceNumber=subdomains[kr].ReferenceNumber-1; 2510 if(kr<0 && kl >=0 && subdomains[kl].ReferenceNumber>=0) 2511 nbk--,subdomains[kl].ReferenceNumber=-1; 2512 if(kl<0 && kr >=0 && subdomains[kr].ReferenceNumber>=0) 2513 nbk--,subdomains[kr].ReferenceNumber=-1; 2515 2514 } 2516 2515 } 2517 2516 long j=0; 2518 for ( i=0;i< NbSubDomains;i++)2519 if((-subdomains[i]. ref) %2) { // good2517 for ( i=0;i<nbsubdomains;i++) 2518 if((-subdomains[i].ReferenceNumber) %2) { // good 2520 2519 if(i != j) 2521 2520 Exchange(subdomains[i],subdomains[j]); … … 2530 2529 }//while (t) 2531 2530 } 2532 if(verbose>4) printf(" Number of removes subdomains (OutSideMesh) = %i\n", NbSubDomains-j);2533 NbSubDomains=j;2531 if(verbose>4) printf(" Number of removes subdomains (OutSideMesh) = %i\n",nbsubdomains-j); 2532 nbsubdomains=j; 2534 2533 } 2535 2534 … … 2539 2538 else 2540 2539 { // find the head for all sub domaine 2541 if (Gh. NbSubDomains != NbSubDomains && subdomains)2540 if (Gh.nbsubdomains != nbsubdomains && subdomains) 2542 2541 delete [] subdomains, subdomains=0; 2543 2542 if (! subdomains ) 2544 subdomains = new SubDomain[ Gh. NbSubDomains];2545 NbSubDomains =Gh.NbSubDomains;2543 subdomains = new SubDomain[ Gh.nbsubdomains]; 2544 nbsubdomains =Gh.nbsubdomains; 2546 2545 long err=0; 2547 2546 CreateSingleVertexToTriangleConnectivity(); … … 2552 2551 mark[it]=triangles[it].link ? -1 : -2; 2553 2552 long inew =0; 2554 for (int i=0;i< NbSubDomains;i++) {2553 for (int i=0;i<nbsubdomains;i++) { 2555 2554 GeometricalEdge &eg = *Gh.subdomains[i].edge; 2556 subdomains[i]. ref = Gh.subdomains[i].ref;2557 int ssdlab = subdomains[i]. ref;2555 subdomains[i].ReferenceNumber = Gh.subdomains[i].ReferenceNumber; 2556 int ssdlab = subdomains[i].ReferenceNumber; 2558 2557 // by carefull is not easy to find a edge create from a GeometricalEdge 2559 2558 // see routine MakeGeometricalEdgeToEdge … … 2564 2563 BamgVertex * v0 = e(0),*v1 = e(1); 2565 2564 Triangle *t = v0->t; 2566 int sens = Gh.subdomains[i].sens;2567 // test if ge and e is in the same sens2568 if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) sens = -sens;2569 subdomains[i]. sens = sens;2565 int direction = Gh.subdomains[i].direction; 2566 // test if ge and e is in the same direction 2567 if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) direction = -direction ; 2568 subdomains[i].direction = direction; 2570 2569 subdomains[i].edge = &e; 2571 if (!t || ! sens){2572 ISSMERROR("!t || ! sens");2570 if (!t || !direction){ 2571 ISSMERROR("!t || !direction"); 2573 2572 } 2574 2573 … … 2580 2579 } 2581 2580 if (ta.EdgeVertex(0) == v1) { // ok we find the edge 2582 if ( sens>0)2581 if (direction>0) 2583 2582 subdomains[i].head=t=Adj(ta); 2584 2583 else … … 2614 2613 } 2615 2614 2616 if (inew < NbSubDomains) {2617 if (verbose>5) printf("WARNING: %i SubDomains are being removed\n", NbSubDomains-inew);2618 NbSubDomains=inew;}2615 if (inew < nbsubdomains) { 2616 if (verbose>5) printf("WARNING: %i SubDomains are being removed\n",nbsubdomains-inew); 2617 nbsubdomains=inew;} 2619 2618 2620 2619 … … 2639 2638 srand(19999999); 2640 2639 NbRef=0; 2641 NbOfTriangleSearchFind =0;2642 NbOfSwapTriangle =0;2643 2640 nbv=0; 2644 2641 maxnbv=maxnbv_in; 2645 2642 nbt=0; 2646 NbOfQuad= 0;2643 nbq = 0; 2647 2644 maxnbt=2*maxnbv_in-2; 2648 NbSubDomains=0;2645 nbsubdomains=0; 2649 2646 NbVertexOnBThVertex=0; 2650 2647 NbVertexOnBThEdge=0; … … 2680 2677 NbVerticesOnGeomEdge=0; 2681 2678 subdomains=NULL; 2682 NbSubDomains=0;2679 nbsubdomains=0; 2683 2680 } 2684 2681 /*}}}1*/ … … 2927 2924 { 2928 2925 Edge * ei = edges+i; 2929 GeometricalEdge * onGeometry = ei->onGeometry;2930 e[Gh.Number( onGeometry)] = ei;2926 GeometricalEdge *GeometricalEdgeHook = ei->GeometricalEdgeHook; 2927 e[Gh.Number(GeometricalEdgeHook)] = ei; 2931 2928 } 2932 2929 for ( i=0;i<nbe ; i++) 2933 2930 for (int ii=0;ii<2;ii++) { 2934 2931 Edge * ei = edges+i; 2935 GeometricalEdge * onGeometry = ei->onGeometry;2932 GeometricalEdge *GeometricalEdgeHook = ei->GeometricalEdgeHook; 2936 2933 int j= ii; 2937 while (!(* onGeometry)[j].Required()) {2938 Adj( onGeometry,j); // next geom edge2934 while (!(*GeometricalEdgeHook)[j].Required()) { 2935 Adj(GeometricalEdgeHook,j); // next geom edge 2939 2936 j=1-j; 2940 if (e[Gh.Number( onGeometry)]) break; // optimisation2941 e[Gh.Number( onGeometry)] = ei;2937 if (e[Gh.Number(GeometricalEdgeHook)]) break; // optimisation 2938 e[Gh.Number(GeometricalEdgeHook)] = ei; 2942 2939 } 2943 2940 } … … 2989 2986 triangles[i].SetHidden(j),kk++; 2990 2987 } 2991 NbOfQuad= kk;2988 nbq = kk; 2992 2989 if (verbose>2){ 2993 printf(" number of quadrilaterals = %i\n", NbOfQuad);2994 printf(" number of triangles = %i\n",nbt-nbtout- NbOfQuad*2);2990 printf(" number of quadrilaterals = %i\n",nbq); 2991 printf(" number of triangles = %i\n",nbt-nbtout- nbq*2); 2995 2992 printf(" number of outside triangles = %i\n",nbtout); 2996 2993 } … … 3109 3106 for (i=0;i<Bh.nbv;i++){ 3110 3107 BamgVertex &bv=Bh[i]; 3111 if (!bv. onGeometry){3108 if (!bv.GeometricalEdgeHook){ 3112 3109 vertices[nbv].r = bv.r; 3113 3110 vertices[nbv++].m = bv.m; … … 3218 3215 BamgVertex * pvA=&vA, * pvB=&vB; 3219 3216 if (vA.vint == IsVertexOnVertex){ 3220 pA=vA. onBackgroundVertex;3217 pA=vA.BackgroundVertexHook; 3221 3218 } 3222 3219 else if (vA.vint == IsVertexOnEdge){ 3223 pA=vA. onBackgroundEdge->be;3224 tA=vA. onBackgroundEdge->abcisse;3220 pA=vA.BackgroundEdgeHook->be; 3221 tA=vA.BackgroundEdgeHook->abcisse; 3225 3222 } 3226 3223 else { … … 3229 3226 3230 3227 if (vB.vint == IsVertexOnVertex){ 3231 pB=vB. onBackgroundVertex;3228 pB=vB.BackgroundVertexHook; 3232 3229 } 3233 3230 else if(vB.vint == IsVertexOnEdge){ 3234 pB=vB. onBackgroundEdge->be;3235 tB=vB. onBackgroundEdge->abcisse;3231 pB=vB.BackgroundEdgeHook->be; 3232 tB=vB.BackgroundEdgeHook->abcisse; 3236 3233 } 3237 3234 else { … … 3255 3252 if( vA.vint == IsVertexOnEdge) 3256 3253 { // find the start edge 3257 e = vA. onBackgroundEdge->be;3254 e = vA.BackgroundEdgeHook->be; 3258 3255 3259 3256 } … … 3265 3262 Exchange(pvA,pvB); 3266 3263 Exchange(A,B); 3267 e = vB. onBackgroundEdge->be;3264 e = vB.BackgroundEdgeHook->be; 3268 3265 3269 3266 } … … 3272 3269 } 3273 3270 3274 // find the direction of walking with sensof edge and pA,PB;3271 // find the direction of walking with direction of edge and pA,PB; 3275 3272 R2 AB=B-A; 3276 3273 3277 3274 double cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB); 3278 3275 int kkk=0; 3279 int sens= (cosE01AB>0) ? 1 : 0;3276 int direction = (cosE01AB>0) ? 1 : 0; 3280 3277 3281 3278 // double l=0; // length of the edge AB … … 3292 3289 double te0; 3293 3290 // we suppose take the curve's abcisse 3294 for ( eee=e,iii= sens,te0=tA;3291 for ( eee=e,iii=direction,te0=tA; 3295 3292 eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ; 3296 3293 neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) { … … 3362 3359 for (i=0;i<nbv;i++) ordre[i]=0; 3363 3360 3364 //Initialize NbSubDomains3365 NbSubDomains =0;3361 //Initialize nbsubdomains 3362 nbsubdomains =0; 3366 3363 3367 3364 /* generation of triangles adjacency*/ … … 3619 3616 for (i=0;i<nbe;i++){ 3620 3617 /*If the current mesh edge is on Geometry*/ 3621 if(edges[i]. onGeometry){3618 if(edges[i].GeometricalEdgeHook){ 3622 3619 for(int j=0;j<2;j++){ 3623 3620 /*Go through the edges adjacent to current edge (if on the same curve)*/ … … 3625 3622 /*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/ 3626 3623 /*Check that the 2 vertices are on geometry AND required*/ 3627 if(!edges[i][j]. onGeometry->IsRequiredVertex()){3624 if(!edges[i][j].GeometricalEdgeHook->IsRequiredVertex()){ 3628 3625 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); 3629 printf("This edge is on geometrical edge number %i\n",Gh.Number(edges[i]. onGeometry)+1);3630 if (edges[i][j]. onGeometry->OnGeomVertex())3631 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);3632 else if (edges[i][j]. onGeometry->OnGeomEdge())3633 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);3626 printf("This edge is on geometrical edge number %i\n",Gh.Number(edges[i].GeometricalEdgeHook)+1); 3627 if (edges[i][j].GeometricalEdgeHook->OnGeomVertex()) 3628 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].GeometricalEdgeHook->gv)+1); 3629 else if (edges[i][j].GeometricalEdgeHook->OnGeomEdge()) 3630 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].GeometricalEdgeHook->ge)+1); 3634 3631 else 3635 printf("Its pointer is %p\n",edges[i][j]. onGeometry);3632 printf("Its pointer is %p\n",edges[i][j].GeometricalEdgeHook); 3636 3633 3637 3634 printf("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n"); … … 3655 3652 for ( it=0;it<nbt;it++) 3656 3653 renu[it]=-1; // outside triangle 3657 for ( i=0;i< NbSubDomains;i++)3654 for ( i=0;i<nbsubdomains;i++) 3658 3655 { 3659 3656 t=t0=subdomains[i].head; … … 3690 3687 triangles[it].Renumbering(triangles,te,renu); 3691 3688 3692 for ( i=0;i< NbSubDomains;i++)3689 for ( i=0;i<nbsubdomains;i++) 3693 3690 subdomains[i].head=triangles+renu[Number(subdomains[i].head)]; 3694 3691 … … 3989 3986 if (tstart[i] != &vide) // not a boundary vertex 3990 3987 delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega)); 3991 if (! NbOfQuad)3988 if (!nbq) 3992 3989 for ( i=0;i<nbv;i++) 3993 3990 if (tstart[i] != &vide) // not a boundary vertex … … 4100 4097 long newnbt=0,newnbv=0; 4101 4098 long * kedge = 0; 4102 long new NbOfQuad=NbOfQuad;4099 long newnbq=nbq; 4103 4100 long * ksplit= 0, * ksplitarray=0; 4104 4101 long kkk=0; … … 4117 4114 Triangle * lastT = triangles + nbt; 4118 4115 for (i=0;i<nbe;i++) 4119 if(edges[i]. onGeometry) NbEdgeOnGeom++;4116 if(edges[i].GeometricalEdgeHook) NbEdgeOnGeom++; 4120 4117 long newnbe=nbe+nbe; 4121 4118 // long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex; … … 4144 4141 long ferr=0; 4145 4142 for (i=0;i<nbe;i++) 4146 newedges[ie]. onGeometry=0;4143 newedges[ie].GeometricalEdgeHook=0; 4147 4144 4148 4145 for (i=0;i<nbe;i++) 4149 4146 { 4150 GeometricalEdge *ong = edges[i]. onGeometry;4147 GeometricalEdge *ong = edges[i].GeometricalEdgeHook; 4151 4148 4152 4149 newedges[ie]=edges[i]; … … 4167 4164 ISSMERROR("!edgesGtoB"); 4168 4165 } 4169 ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i]. onGeometry)],4166 ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].GeometricalEdgeHook)], 4170 4167 edges[i][0],edges[i][1],0.5,vertices[k], 4171 4168 newVertexOnBThEdge[kvb], 4172 4169 newVerticesOnGeomEdge[kvg++]); 4173 vertices[k].ReferenceNumber= edges[i]. ref;4170 vertices[k].ReferenceNumber= edges[i].ReferenceNumber; 4174 4171 vertices[k].DirOfSearch = NoDirOfSearch; 4175 4172 ; … … 4187 4184 0.5,vertices[k],newVerticesOnGeomEdge[kvg++]); 4188 4185 // vertices[k].i = toI2( vertices[k].r); 4189 vertices[k].ReferenceNumber = edges[i]. ref;4186 vertices[k].ReferenceNumber = edges[i].ReferenceNumber; 4190 4187 vertices[k].DirOfSearch = NoDirOfSearch; 4191 4188 vertices[k].m = Metric(0.5,edges[i][0],0.5,edges[i][1]); … … 4196 4193 vertices[k].r = ((R2) edges[i][0] + (R2) edges[i][1] )*0.5; 4197 4194 vertices[k].m = Metric(0.5,edges[i][0],0.5,edges[i][1]); 4198 vertices[k]. onGeometry= 0;4195 vertices[k].GeometricalEdgeHook = 0; 4199 4196 } 4200 4197 //vertices[k].i = toI2( vertices[k].r); … … 4202 4199 R2 AA = (A+AB)*0.5; 4203 4200 R2 BB = (AB+B)*0.5; 4204 vertices[k].ReferenceNumber = edges[i]. ref;4201 vertices[k].ReferenceNumber = edges[i].ReferenceNumber; 4205 4202 vertices[k].DirOfSearch = NoDirOfSearch; 4206 4203 4207 newedges[ie]. onGeometry= Gh.Containing(AA,ong);4204 newedges[ie].GeometricalEdgeHook = Gh.Containing(AA,ong); 4208 4205 newedges[ie++].v[1]=vertices+k; 4209 4206 … … 4211 4208 newedges[ie].adj[0]=newedges + ie -1; 4212 4209 newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ; 4213 newedges[ie]. onGeometry= Gh.Containing(BB,ong);4210 newedges[ie].GeometricalEdgeHook = Gh.Containing(BB,ong); 4214 4211 newedges[ie++].v[0]=vertices+k; 4215 4212 k++; … … 4349 4346 } 4350 4347 // now do the element split 4351 new NbOfQuad = 4*NbOfQuad;4348 newnbq = 4*nbq; 4352 4349 nbv = k; 4353 4350 kkk = nbt; … … 4517 4514 } 4518 4515 4519 if (kk==6) new NbOfQuad+=3;4516 if (kk==6) newnbq+=3; 4520 4517 for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk; 4521 4518 ksplit[i]= nbt; // save last adresse of the new triangles … … 4544 4541 edges = newedges; 4545 4542 nbe = newnbe; 4546 NbOfQuad = newNbOfQuad;4547 4548 for (i=0;i< NbSubDomains;i++)4543 nbq = newnbq; 4544 4545 for (i=0;i<nbsubdomains;i++) 4549 4546 { 4550 4547 long k = subdomains[i].edge- edges; … … 4566 4563 4567 4564 if (verbose>2){ 4568 printf(" number of quadrilaterals = %i\n", NbOfQuad);4569 printf(" number of triangles = %i\n",nbt-nbtout- NbOfQuad*2);4565 printf(" number of quadrilaterals = %i\n",nbq); 4566 printf(" number of triangles = %i\n",nbt-nbtout- nbq*2); 4570 4567 printf(" number of outside triangles = %i\n",nbtout); 4571 4568 } … … 4607 4604 BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]]; 4608 4605 BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]]; 4609 if (v0. onGeometry && v1.onGeometry){4606 if (v0.GeometricalEdgeHook && v1.GeometricalEdgeHook){ 4610 4607 R2 P= ((R2) v0 + (R2) v1)*0.5; 4611 4608 if ( nbv<maxnbv) { … … 4750 4747 if (t->det<0) // outside triangle 4751 4748 dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop; 4752 // NbOfTriangleSearchFind += counter;4753 4749 return t; 4754 4750 } … … 4777 4773 4778 4774 //loop over all subdomains 4779 for (int i=0;i< NbSubDomains;i++){4775 for (int i=0;i<nbsubdomains;i++){ 4780 4776 4781 4777 //first triangle of the subdomain i … … 4814 4810 4815 4811 int i,j,k; 4816 int NbOfCurves=0,NbNewPoints,NbEdgeCurve;4812 int nbcurves=0,NbNewPoints,NbEdgeCurve; 4817 4813 double lcurve,lstep,s; 4818 4814 const int MaxSubEdge = 10; … … 4890 4886 long NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge; 4891 4887 Gh.UnMarkEdges(); 4892 NbOfCurves=0;4888 nbcurves=0; 4893 4889 4894 4890 //go through the edges of the geometry … … 4922 4918 edges[nbe].v[0]=a->to; 4923 4919 edges[nbe].v[1]=b->to;; 4924 edges[nbe]. ref = e->ref;4925 edges[nbe]. onGeometry= e;4920 edges[nbe].ReferenceNumber = e->ReferenceNumber; 4921 edges[nbe].GeometricalEdgeHook = e; 4926 4922 edges[nbe].adj[0] = 0; 4927 4923 edges[nbe].adj[1] = 0; … … 5021 5017 vb = &vertices[nbv++]; 5022 5018 vb->m = Metric(aa,a->m,bb,b->m); 5023 vb->ReferenceNumber = e-> ref;5019 vb->ReferenceNumber = e->ReferenceNumber; 5024 5020 vb->DirOfSearch =NoDirOfSearch; 5025 5021 double abcisse = k ? bb : aa; … … 5031 5027 edges[nbe].v[0]=va; 5032 5028 edges[nbe].v[1]=vb; 5033 edges[nbe]. ref =e->ref;5034 edges[nbe]. onGeometry= e;5029 edges[nbe].ReferenceNumber =e->ReferenceNumber; 5030 edges[nbe].GeometricalEdgeHook = e; 5035 5031 edges[nbe].adj[0] = PreviousNewEdge; 5036 5032 if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe]; … … 5053 5049 if(!kstep){ 5054 5050 NbVerticesOnGeomEdge0 += NbNewPoints; 5055 NbOfCurves++;5051 nbcurves++; 5056 5052 } 5057 5053 nbvend=nbv+NbNewPoints; … … 5061 5057 edges[nbe].v[0]=va; 5062 5058 edges[nbe].v[1]=vb; 5063 edges[nbe]. ref = e->ref;5064 edges[nbe]. onGeometry= e;5059 edges[nbe].ReferenceNumber = e->ReferenceNumber; 5060 edges[nbe].GeometricalEdgeHook = e; 5065 5061 edges[nbe].adj[0] = PreviousNewEdge; 5066 5062 edges[nbe].adj[1] = 0; … … 5146 5142 this->Init(imaxnbv); 5147 5143 BTh.SetVertexFieldOn(); 5148 int* bcurve = new int[Gh. NbOfCurves]; //5144 int* bcurve = new int[Gh.nbcurves]; // 5149 5145 5150 5146 /* There are 2 ways to make the loop … … 5196 5192 Gh.UnMarkEdges(); 5197 5193 int bfind=0; 5198 for (int i=0;i<Gh. NbOfCurves;i++) bcurve[i]=-1;5194 for (int i=0;i<Gh.nbcurves;i++) bcurve[i]=-1; 5199 5195 5200 5196 /*Loop over the backgrounf mesh BTh edges*/ … … 5206 5202 5207 5203 /* If one of the vertex is required we are in a new curve*/ 5208 if (ei[je]. onGeometry->IsRequiredVertex()){5204 if (ei[je].GeometricalEdgeHook->IsRequiredVertex()){ 5209 5205 5210 5206 /*Get curve number*/ 5211 int nc=ei. onGeometry->CurveNumber;5207 int nc=ei.GeometricalEdgeHook->CurveNumber; 5212 5208 5213 5209 //printf("Dealing with curve number %i\n",nc); 5214 //printf("edge on geometry is same as GhCurve? %s\n",(ei. onGeometry==Gh.curves[nc].be || ei.onGeometry==Gh.curves[nc].ee)?"yes":"no");5215 //if(ei. onGeometry==Gh.curves[nc].be || ei.onGeometry==Gh.curves[nc].ee){5216 // printf("Do we have the right extremity? curve first vertex -> %s\n",((GeometricalVertex *)*ei[je]. onGeometry==&(*Gh.curves[nc].be)[Gh.curves[nc].kb])?"yes":"no");5217 // printf("Do we have the right extremity? curve last vertex -> %s\n",((GeometricalVertex *)*ei[je]. onGeometry==&(*Gh.curves[nc].ee)[Gh.curves[nc].ke])?"yes":"no");5210 //printf("edge on geometry is same as GhCurve? %s\n",(ei.GeometricalEdgeHook==Gh.curves[nc].be || ei.GeometricalEdgeHook==Gh.curves[nc].ee)?"yes":"no"); 5211 //if(ei.GeometricalEdgeHook==Gh.curves[nc].be || ei.GeometricalEdgeHook==Gh.curves[nc].ee){ 5212 // printf("Do we have the right extremity? curve first vertex -> %s\n",((GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].be)[Gh.curves[nc].kb])?"yes":"no"); 5213 // printf("Do we have the right extremity? curve last vertex -> %s\n",((GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].ee)[Gh.curves[nc].ke])?"yes":"no"); 5218 5214 //} 5219 5215 //BUG FIX from original bamg 5220 5216 /*Check that we are on the same edge and right vertex (0 or 1) */ 5221 if(ei. onGeometry==Gh.curves[nc].be && (GeometricalVertex *)*ei[je].onGeometry==&(*Gh.curves[nc].be)[Gh.curves[nc].kb]){5217 if(ei.GeometricalEdgeHook==Gh.curves[nc].be && (GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].be)[Gh.curves[nc].kb]){ 5222 5218 bcurve[nc]=iedge*2+je; 5223 5219 bfind++; 5224 5220 } 5225 else if ((ei. onGeometry==Gh.curves[nc].ee && (GeometricalVertex *)*ei[je].onGeometry==&(*Gh.curves[nc].ee)[Gh.curves[nc].ke]) && bcurve[nc]==-1){5221 else if ((ei.GeometricalEdgeHook==Gh.curves[nc].ee && (GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].ee)[Gh.curves[nc].ke]) && bcurve[nc]==-1){ 5226 5222 bcurve[nc]=iedge*2+je; 5227 5223 bfind++; … … 5230 5226 } 5231 5227 } 5232 if (bfind!=Gh. NbOfCurves) ISSMERROR("problem generating number of curves (Gh.NbOfCurves=%i bfind=%i)",Gh.NbOfCurves,bfind);5228 if (bfind!=Gh.nbcurves) ISSMERROR("problem generating number of curves (Gh.nbcurves=%i bfind=%i)",Gh.nbcurves,bfind); 5233 5229 5234 5230 // method in 2 + 1 step … … 5247 5243 5248 5244 /*Go through all geometrical curve*/ 5249 for (int icurve=0;icurve<Gh. NbOfCurves;icurve++){5245 for (int icurve=0;icurve<Gh.nbcurves;icurve++){ 5250 5246 5251 5247 /*Get edge and vertex (index) of background mesh on this curve*/ … … 5282 5278 int k0equi=jedgeequi,k1equi; 5283 5279 Edge * peequi= BTh.edges+iedgeequi; 5284 GeometricalEdge *ongequi = peequi-> onGeometry;5280 GeometricalEdge *ongequi = peequi->GeometricalEdgeHook; 5285 5281 5286 5282 double sNew=Lstep;// abscisse of the new points (phase==1) 5287 5283 L=0;// length of the curve 5288 5284 long i=0;// index of new points on the curve 5289 register GeometricalVertex * GA0 = *(*peequi)[k0equi]. onGeometry;5285 register GeometricalVertex * GA0 = *(*peequi)[k0equi].GeometricalEdgeHook; 5290 5286 BamgVertex *A0; 5291 5287 A0 = GA0->to; // the vertex in new mesh … … 5297 5293 ISSMASSERT(A0-vertices>=0 && A0-vertices<nbv); 5298 5294 if(ongequi->Required()){ 5299 GeometricalVertex *GA1 = *(*peequi)[1-k0equi]. onGeometry;5295 GeometricalVertex *GA1 = *(*peequi)[1-k0equi].GeometricalEdgeHook; 5300 5296 A1 = GA1->to; // 5301 5297 } … … 5306 5302 k1 = 1-k0; // next vertex of the edge 5307 5303 k1equi= 1 - k0equi; 5308 ISSMASSERT(pe && ee. onGeometry);5309 ee. onGeometry->SetMark();5304 ISSMASSERT(pe && ee.GeometricalEdgeHook); 5305 ee.GeometricalEdgeHook->SetMark(); 5310 5306 BamgVertex & v0=ee[0], & v1=ee[1]; 5311 5307 R2 AB=(R2)v1-(R2)v0; … … 5341 5337 VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save 5342 5338 ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); 5343 A1->ReferenceNumber = eeequi. ref;5339 A1->ReferenceNumber = eeequi.ReferenceNumber; 5344 5340 A1->DirOfSearch =NoDirOfSearch; 5345 e-> onGeometry= ongequi;5341 e->GeometricalEdgeHook = ongequi; 5346 5342 e->v[0]=A0; 5347 5343 e->v[1]=A1; 5348 e-> ref = eeequi.ref;5344 e->ReferenceNumber = eeequi.ReferenceNumber; 5349 5345 e->adj[0]=PreviousNewEdge; 5350 5346 … … 5358 5354 5359 5355 //some checks 5360 ISSMASSERT(ee. onGeometry->CurveNumber==ei.onGeometry->CurveNumber);5361 if (ee[k1]. onGeometry->IsRequiredVertex()) {5362 ISSMASSERT(eeequi[k1equi]. onGeometry->IsRequiredVertex());5363 register GeometricalVertex * GA1 = *eeequi[k1equi]. onGeometry;5356 ISSMASSERT(ee.GeometricalEdgeHook->CurveNumber==ei.GeometricalEdgeHook->CurveNumber); 5357 if (ee[k1].GeometricalEdgeHook->IsRequiredVertex()) { 5358 ISSMASSERT(eeequi[k1equi].GeometricalEdgeHook->IsRequiredVertex()); 5359 register GeometricalVertex * GA1 = *eeequi[k1equi].GeometricalEdgeHook; 5364 5360 A1=GA1->to;// the vertex in new mesh 5365 5361 ISSMASSERT(A1-vertices>=0 && A1-vertices<nbv); … … 5379 5375 if (phase){ // construction of the last edge 5380 5376 Edge* e=edges + nbe++; 5381 e-> onGeometry= ongequi;5377 e->GeometricalEdgeHook = ongequi; 5382 5378 e->v[0]=A0; 5383 5379 e->v[1]=A1; 5384 e-> ref = peequi->ref;5380 e->ReferenceNumber = peequi->ReferenceNumber; 5385 5381 e->adj[0]=PreviousNewEdge; 5386 5382 e->adj[1]=0; … … 5649 5645 AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.vint][0]); 5650 5646 BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 5651 // we turn around a in the direct sens5647 // we turn around a in the direct direction 5652 5648 5653 5649 Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1; … … 5816 5812 if(!ToSwap) tt1 = Next(tt2); 5817 5813 } 5818 else { // changement de sens5814 else { // changement de direction 5819 5815 ret = -1; 5820 5816 Exchange(pva,pvb); -
issm/trunk/src/c/objects/Bamg/Mesh.h
r5146 r5148 34 34 SubDomain *subdomains; 35 35 long NbRef; // counter of ref on the this class if 0 we can delete 36 long maxnbv,maxnbt; // nombre max de sommets , de triangles 37 long nbv,nbt,nbe; // nb of vertices, of triangles, of edges 38 long NbOfQuad; // nb of quadrangle 39 long NbSubDomains; 36 long maxnbv,maxnbt; // nombre max de sommets , de triangles 37 long nbv,nbt,nbe,nbq; // nb of vertices, of triangles, of edges and quadrilaterals 38 long nbsubdomains; 40 39 long nbtout; // Nb of oudeside triangle 41 long NbOfTriangleSearchFind; 42 long NbOfSwapTriangle; 40 41 R2 pmin,pmax; // extrema 42 double coefIcoor; // coef to integer Icoor1; 43 ListofIntersectionTriangles lIntTria; 44 43 45 long NbVerticesOnGeomVertex; 44 46 VertexOnGeom *VerticesOnGeomVertex; … … 53 55 long NbCrackedEdges; 54 56 CrackedEdge *CrackedEdges; 55 R2 pmin,pmax; // extrema56 double coefIcoor; // coef to integer Icoor1;57 ListofIntersectionTriangles lIntTria;58 57 59 58 //Constructors/Destructors … … 139 138 } 140 139 inline void SetVertexFieldOn(){ 141 for (int i=0;i<nbv;i++) vertices[i]. onGeometry=NULL;140 for (int i=0;i<nbv;i++) vertices[i].GeometricalEdgeHook=NULL; 142 141 for (int j=0;j<NbVerticesOnGeomVertex;j++) VerticesOnGeomVertex[j].SetOn(); 143 142 for (int k=0;k<NbVerticesOnGeomEdge;k++ ) VerticesOnGeomEdge[k].SetOn(); 144 143 } 145 144 inline void SetVertexFieldOnBTh(){ 146 for (int i=0;i<nbv;i++) vertices[i]. onGeometry=NULL;145 for (int i=0;i<nbv;i++) vertices[i].GeometricalEdgeHook=NULL; 147 146 for (int j=0;j<NbVertexOnBThVertex;j++) VertexOnBThVertex[j].SetOnBTh(); 148 147 for (int k=0;k<NbVertexOnBThEdge;k++ ) VertexOnBThEdge[k].SetOnBTh(); -
issm/trunk/src/c/objects/Bamg/Metric.cpp
r5016 r5148 156 156 { 157 157 k=k/2; 158 // we begin by the end to walk in the correct sensfrom a to b158 // we begin by the end to walk in the correct direction from a to b 159 159 // due to the stack 160 160 Ms1[level]=Mi; -
issm/trunk/src/c/objects/Bamg/QuadTree.cpp
r5120 r5148 388 388 { 389 389 I2 i2 = b->v[k]->i; 390 // try if is in the right sens--390 // try if is in the right direction -- 391 391 h0 = NORM(iplus,i2.x,jplus,i2.y); 392 392 if (h0 <h) { … … 415 415 NbVerticesSearch++; 416 416 I2 i2 = b->v[k]->i; 417 // if good senswhen try --417 // if good direction when try -- 418 418 419 419 h0 = NORM(iplus,i2.x,jplus,i2.y); -
issm/trunk/src/c/objects/Bamg/R2.h
r5124 r5148 91 91 } 92 92 template <class R,class RR> 93 inline R NormeInfini (const P2<R,RR> x) {94 return Max(Abs(x.x),Abs(x.y)) ;95 }96 template <class R,class RR>97 93 inline RR Norme2_2 (const P2<R,RR> x) { 98 94 return (RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y ; -
issm/trunk/src/c/objects/Bamg/SubDomain.h
r5095 r5148 12 12 13 13 class SubDomain { 14 14 15 public: 15 Triangle * head; 16 long ref; 17 int sens; // -1 or 1 18 Edge* edge; // to geometrical 16 17 Triangle *head; 18 long ReferenceNumber; 19 int direction; // -1 or 1 20 Edge *edge; // to geometrical 19 21 20 22 //Methods -
issm/trunk/src/c/objects/Bamg/Triangle.cpp
r5143 r5148 19 19 ISSMERROR("i>=nbv || j>=nbv || k>=nbv"); 20 20 } 21 TriaVertices[0]=v+i;22 TriaVertices[1]=v+j;23 TriaVertices[2]=v+k;24 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;25 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;21 vertices[0]=v+i; 22 vertices[1]=v+j; 23 vertices[2]=v+k; 24 adj[0]=adj[1]=adj[2]=0; 25 AdjEdgeNumber[0]=AdjEdgeNumber[1]=AdjEdgeNumber[2]=0; 26 26 det=0; 27 27 } … … 29 29 /*FUNCTION Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2) {{{1*/ 30 30 Triangle::Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2){ 31 TriaVertices[0]=v0;32 TriaVertices[1]=v1;33 TriaVertices[2]=v2;34 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;35 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;31 vertices[0]=v0; 32 vertices[1]=v1; 33 vertices[2]=v2; 34 adj[0]=adj[1]=adj[2]=0; 35 AdjEdgeNumber[0]=AdjEdgeNumber[1]=AdjEdgeNumber[2]=0; 36 36 if (v0) det=0; 37 37 else { … … 48 48 49 49 printf("Triangle:\n"); 50 printf(" TriaVertices pointer towards three vertices\n");51 printf(" TriaVertices[0] TriaVertices[1] TriaVertices[2] = %p %p %p\n",TriaVertices[0],TriaVertices[1],TriaVertices[2]);52 printf(" TriaAdjTrianglespointer towards three adjacent triangles\n");53 printf(" TriaAdjTriangles[0] TriaAdjTriangles[1] TriaAdjTriangles[2] = %p %p %p\n",TriaAdjTriangles[0],TriaAdjTriangles[1],TriaAdjTriangles[2]);50 printf(" vertices pointer towards three vertices\n"); 51 printf(" vertices[0] vertices[1] vertices[2] = %p %p %p\n",vertices[0],vertices[1],vertices[2]); 52 printf(" adj pointer towards three adjacent triangles\n"); 53 printf(" adj[0] adj[1] adj[2] = %p %p %p\n",adj[0],adj[1],adj[2]); 54 54 printf(" det (integer triangle determinant) = %i\n",det); 55 55 if (link){ … … 62 62 printf("\nThree vertices:\n"); 63 63 for(i=0;i<3;i++){ 64 if ( TriaVertices[i]){65 TriaVertices[i]->Echo();64 if (vertices[i]){ 65 vertices[i]->Echo(); 66 66 } 67 67 else{ … … 99 99 k++; 100 100 //Get ttc, adjacent triangle of t with respect to vertex j 101 ttc = t-> TriaAdjTriangles[j];101 ttc = t->adj[j]; 102 102 //is the current triangle inside or outside? 103 103 outside = !ttc->link; … … 108 108 t = ttc; 109 109 //NextEdge[3] = {1,2,0}; 110 jc = NextEdge[t-> TriaAdjSharedEdge[j]&3];110 jc = NextEdge[t->AdjEdgeNumber[j]&3]; 111 111 j = NextEdge[jc]; 112 112 … … 132 132 133 133 // initialize tp, jp the previous triangle & edge 134 Triangle *tp= TriaAdjTriangles[jp];135 jp = TriaAdjSharedEdge[jp]&3;134 Triangle *tp=adj[jp]; 135 jp = AdjEdgeNumber[jp]&3; 136 136 do { 137 137 while (t->swap(j,koption)){ … … 139 139 NbSwap++; 140 140 k++; 141 t= tp-> TriaAdjTriangles[jp]; // set unchange t qnd j for previous triangles142 j= NextEdge[tp-> TriaAdjSharedEdge[jp]&3];141 t= tp->adj[jp]; // set unchange t qnd j for previous triangles 142 j= NextEdge[tp->AdjEdgeNumber[jp]&3]; 143 143 } 144 144 // end on this Triangle … … 146 146 jp = NextEdge[j]; 147 147 148 t= tp-> TriaAdjTriangles[jp]; // set unchange t qnd j for previous triangles149 j= NextEdge[tp-> TriaAdjSharedEdge[jp]&3];148 t= tp->adj[jp]; // set unchange t qnd j for previous triangles 149 j= NextEdge[tp->AdjEdgeNumber[jp]&3]; 150 150 151 151 } while( t != this); … … 159 159 if (link) { 160 160 int a=-1; 161 if ( TriaAdjSharedEdge[0] & 16 ) a=0;162 if ( TriaAdjSharedEdge[1] & 16 ) a=1;163 if ( TriaAdjSharedEdge[2] & 16 ) a=2;161 if (AdjEdgeNumber[0] & 16 ) a=0; 162 if (AdjEdgeNumber[1] & 16 ) a=1; 163 if (AdjEdgeNumber[2] & 16 ) a=2; 164 164 if (a>=0) { 165 t = TriaAdjTriangles[a];165 t = adj[a]; 166 166 // if (t-this<0) return 0; 167 v2 = TriaVertices[VerticesOfTriangularEdge[a][0]];168 v0 = TriaVertices[VerticesOfTriangularEdge[a][1]];169 v1 = TriaVertices[OppositeEdge[a]];170 v3 = t-> TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];167 v2 = vertices[VerticesOfTriangularEdge[a][0]]; 168 v0 = vertices[VerticesOfTriangularEdge[a][1]]; 169 v1 = vertices[OppositeEdge[a]]; 170 v3 = t->vertices[OppositeEdge[AdjEdgeNumber[a]&3]]; 171 171 } 172 172 } … … 177 177 double Triangle::QualityQuad(int a,int option) const{ 178 178 double q; 179 if (!link || TriaAdjSharedEdge[a] &4)179 if (!link || AdjEdgeNumber[a] &4) 180 180 q= -1; 181 181 else { 182 Triangle * t = TriaAdjTriangles[a];182 Triangle * t = adj[a]; 183 183 if (t-this<0) q= -1;// because we do 2 times 184 184 else if (!t->link ) q= -1; 185 else if ( TriaAdjSharedEdge[0] & 16 || TriaAdjSharedEdge[1] & 16 || TriaAdjSharedEdge[2] & 16 || t->TriaAdjSharedEdge[0] & 16 || t->TriaAdjSharedEdge[1] & 16 || t->TriaAdjSharedEdge[2] & 16 )185 else if (AdjEdgeNumber[0] & 16 || AdjEdgeNumber[1] & 16 || AdjEdgeNumber[2] & 16 || t->AdjEdgeNumber[0] & 16 || t->AdjEdgeNumber[1] & 16 || t->AdjEdgeNumber[2] & 16 ) 186 186 q= -1; 187 187 else if(option){ 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]];188 const BamgVertex & v2 = *vertices[VerticesOfTriangularEdge[a][0]]; 189 const BamgVertex & v0 = *vertices[VerticesOfTriangularEdge[a][1]]; 190 const BamgVertex & v1 = *vertices[OppositeEdge[a]]; 191 const BamgVertex & v3 = * t->vertices[OppositeEdge[AdjEdgeNumber[a]&3]]; 192 192 q = QuadQuality(v0,v1,v2,v3); // do the float part 193 193 } … … 200 200 void Triangle::Set(const Triangle & rec,const Mesh & Th ,Mesh & ThNew){ 201 201 *this = rec; 202 if ( TriaVertices[0] ) TriaVertices[0] = ThNew.vertices + Th.Number(TriaVertices[0]);203 if ( TriaVertices[1] ) TriaVertices[1] = ThNew.vertices + Th.Number(TriaVertices[1]);204 if ( TriaVertices[2] ) TriaVertices[2] = ThNew.vertices + Th.Number(TriaVertices[2]);205 if( TriaAdjTriangles[0]) TriaAdjTriangles[0] = ThNew.triangles + Th.Number(TriaAdjTriangles[0]);206 if( TriaAdjTriangles[1]) TriaAdjTriangles[1] = ThNew.triangles + Th.Number(TriaAdjTriangles[1]);207 if( TriaAdjTriangles[2]) TriaAdjTriangles[2] = ThNew.triangles + Th.Number(TriaAdjTriangles[2]);202 if ( vertices[0] ) vertices[0] = ThNew.vertices + Th.Number(vertices[0]); 203 if ( vertices[1] ) vertices[1] = ThNew.vertices + Th.Number(vertices[1]); 204 if ( vertices[2] ) vertices[2] = ThNew.vertices + Th.Number(vertices[2]); 205 if(adj[0]) adj[0] = ThNew.triangles + Th.Number(adj[0]); 206 if(adj[1]) adj[1] = ThNew.triangles + Th.Number(adj[1]); 207 if(adj[2]) adj[2] = ThNew.triangles + Th.Number(adj[2]); 208 208 if (link >= Th.triangles && link < Th.triangles + Th.nbt) 209 209 link = ThNew.triangles + Th.Number(link); … … 216 216 if(a/4 !=0) return 0;// arete lock or MarkUnSwap 217 217 218 register Triangle *t1=this,*t2= TriaAdjTriangles[a];// les 2 triangles adjacent219 register short a1=a,a2= TriaAdjSharedEdge[a];// les 2 numero de l arete dans les 2 triangles218 register Triangle *t1=this,*t2=adj[a];// les 2 triangles adjacent 219 register short a1=a,a2=AdjEdgeNumber[a];// les 2 numero de l arete dans les 2 triangles 220 220 if(a2/4 !=0) return 0; // arete lock or MarkUnSwap 221 221 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]];222 register BamgVertex *sa=t1->vertices[VerticesOfTriangularEdge[a1][0]]; 223 register BamgVertex *sb=t1->vertices[VerticesOfTriangularEdge[a1][1]]; 224 register BamgVertex *s1=t1->vertices[OppositeVertex[a1]]; 225 register BamgVertex *s2=t2->vertices[OppositeVertex[a2]]; 226 226 227 227 Icoor2 det1=t1->det , det2=t2->det ; -
issm/trunk/src/c/objects/Bamg/Triangle.h
r5143 r5148 17 17 18 18 private: 19 BamgVertex * TriaVertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer20 Triangle * TriaAdjTriangles[3]; // 3 pointers toward the adjacent triangles21 short TriaAdjSharedEdge[3]; // edge id in the adjacent triangles. The edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 119 BamgVertex *vertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer 20 Triangle *adj[3]; // 3 pointers toward the adjacent triangles 21 short AdjEdgeNumber[3]; // edge id in the adjacent triangles. The edge number 1 is the edge number AdjEdgeNumber[1] in the Adjacent triangle 1 22 22 23 23 public: … … 34 34 35 35 //Operators 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];};36 const BamgVertex & operator[](int i) const {return *vertices[i];}; 37 BamgVertex & operator[](int i) {return *vertices[i];}; 38 const BamgVertex * operator()(int i) const {return vertices[i];}; 39 BamgVertex * & operator()(int i) {return vertices[i];}; 40 40 41 41 //Methods … … 43 43 int swap(short a1,int=0); 44 44 long Optim(short a,int =0); 45 int Locked(int a)const { return TriaAdjSharedEdge[a]&4;}46 int Hidden(int a)const { return TriaAdjSharedEdge[a]&16;}47 int GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;}48 void SetAllFlag(int a,int f){ TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);}45 int Locked(int a)const { return AdjEdgeNumber[a]&4;} 46 int Hidden(int a)const { return AdjEdgeNumber[a]&16;} 47 int GetAllflag(int a){return AdjEdgeNumber[a] & 1020;} 48 void SetAllFlag(int a,int f){AdjEdgeNumber[a] = (AdjEdgeNumber[a] &3) + (1020 & f);} 49 49 double QualityQuad(int a,int option=1) const; 50 short NuEdgeTriangleAdj(int i) const {return TriaAdjSharedEdge[i&3]&3;} // Number of the adjacent edge in adj tria50 short NuEdgeTriangleAdj(int i) const {return AdjEdgeNumber[i&3]&3;} // Number of the adjacent edge in adj tria 51 51 AdjacentTriangle FindBoundaryEdge(int i) const; 52 AdjacentTriangle Adj(int i) const {return AdjacentTriangle( TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);};53 Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];}52 AdjacentTriangle Adj(int i) const {return AdjacentTriangle(adj[i],AdjEdgeNumber[i]&3);}; 53 Triangle* TriangleAdj(int i) const {return adj[i&3];} 54 54 Triangle* Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const ; 55 55 void Renumbering(Triangle *tb,Triangle *te, long *renu){ 56 56 if (link >=tb && link <te) link = tb + renu[link -tb]; 57 if ( TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb];58 if ( TriaAdjTriangles[1] >=tb && TriaAdjTriangles[1] <te) TriaAdjTriangles[1] = tb + renu[TriaAdjTriangles[1]-tb];59 if ( TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb];57 if (adj[0] >=tb && adj[0] <te) adj[0] = tb + renu[adj[0]-tb]; 58 if (adj[1] >=tb && adj[1] <te) adj[1] = tb + renu[adj[1]-tb]; 59 if (adj[2] >=tb && adj[2] <te) adj[2] = tb + renu[adj[2]-tb]; 60 60 } 61 61 void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){ 62 if ( TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb];63 if ( TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb];64 if ( TriaVertices[2] >=vb && TriaVertices[2] <ve) TriaVertices[2] = vb + renu[TriaVertices[2]-vb];62 if (vertices[0] >=vb && vertices[0] <ve) vertices[0] = vb + renu[vertices[0]-vb]; 63 if (vertices[1] >=vb && vertices[1] <ve) vertices[1] = vb + renu[vertices[1]-vb]; 64 if (vertices[2] >=vb && vertices[2] <ve) vertices[2] = vb + renu[vertices[2]-vb]; 65 65 } 66 66 void SetAdjAdj(short a){ 67 67 a &= 3; 68 register Triangle *tt= TriaAdjTriangles[a];69 TriaAdjSharedEdge[a] &= 55; // remove MarkUnSwap70 register short aatt = TriaAdjSharedEdge[a] & 3;68 register Triangle *tt=adj[a]; 69 AdjEdgeNumber [a] &= 55; // remove MarkUnSwap 70 register short aatt = AdjEdgeNumber[a] & 3; 71 71 if(tt){ 72 tt-> TriaAdjTriangles[aatt]=this;73 tt-> TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark72 tt->adj[aatt]=this; 73 tt->AdjEdgeNumber[aatt]=a + (AdjEdgeNumber[a] & 60 ) ;}// Copy all the mark 74 74 } 75 75 void SetAdj2(short a,Triangle *t,short aat){ 76 TriaAdjTriangles[a]=t; //the adjacent triangle to the edge a is t77 TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle76 adj[a]=t; //the adjacent triangle to the edge a is t 77 AdjEdgeNumber[a]=aat; //position of the edge in the adjacent triangle 78 78 if(t) { //if t!=NULL add adjacent triangle to t (this) 79 t-> TriaAdjTriangles[aat]=this;80 t-> TriaAdjSharedEdge[aat]=a;79 t->adj[aat]=this; 80 t->AdjEdgeNumber[aat]=a; 81 81 } 82 82 } 83 83 void SetSingleVertexToTriangleConnectivity() { 84 if ( TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0);85 if ( TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1);86 if ( TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2);84 if (vertices[0]) (vertices[0]->t=this,vertices[0]->vint=0); 85 if (vertices[1]) (vertices[1]->t=this,vertices[1]->vint=1); 86 if (vertices[2]) (vertices[2]->t=this,vertices[2]->vint=2); 87 87 } 88 88 void SetHidden(int a){ 89 89 //Get Adjacent Triangle number a 90 register Triangle* t = TriaAdjTriangles[a];90 register Triangle* t = adj[a]; 91 91 //if it exist 92 92 //C|=D -> C=(C|D) bitwise inclusive OR 93 if(t) t-> TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16;94 TriaAdjSharedEdge[a] |= 16;93 if(t) t->AdjEdgeNumber[AdjEdgeNumber[a] & 3] |=16; 94 AdjEdgeNumber[a] |= 16; 95 95 } 96 96 97 97 void SetLocked(int a){ 98 98 //mark the edge as on Boundary 99 register Triangle * t = TriaAdjTriangles[a];100 t-> TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4;101 TriaAdjSharedEdge[a] |= 4;99 register Triangle * t = adj[a]; 100 t->AdjEdgeNumber[AdjEdgeNumber[a] & 3] |=4; 101 AdjEdgeNumber[a] |= 4; 102 102 } 103 103 void SetMarkUnSwap(int a){ 104 register Triangle * t = TriaAdjTriangles[a];105 t-> TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8;106 TriaAdjSharedEdge[a] |=8 ;104 register Triangle * t = adj[a]; 105 t->AdjEdgeNumber[AdjEdgeNumber[a] & 3] |=8; 106 AdjEdgeNumber[a] |=8 ; 107 107 } 108 108 void SetUnMarkUnSwap(int a){ 109 register Triangle * t = TriaAdjTriangles[a];110 t-> TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32111 TriaAdjSharedEdge[a] &=55 ;109 register Triangle * t = adj[a]; 110 t->AdjEdgeNumber[AdjEdgeNumber[a] & 3] &=55; // 23 + 32 111 AdjEdgeNumber[a] &=55 ; 112 112 } 113 113 void SetDet() { 114 if( TriaVertices[0] && TriaVertices[1] && TriaVertices[2]) det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]);114 if(vertices[0] && vertices[1] && vertices[2]) det = bamg::det(*vertices[0],*vertices[1],*vertices[2]); 115 115 else det = -1; } 116 116 … … 118 118 double qualite() ; 119 119 void Set(const Triangle &,const Mesh &,Mesh &); 120 int In(BamgVertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;}120 int In(BamgVertex *v) const { return vertices[0]==v || vertices[1]==v || vertices[2]==v ;} 121 121 122 122 }; -
issm/trunk/src/c/objects/Bamg/VertexOnEdge.h
r5120 r5148 29 29 30 30 //Methods 31 void SetOnBTh(){v-> onBackgroundEdge=this;v->vint=IsVertexOnEdge;}31 void SetOnBTh(){v->BackgroundEdgeHook=this;v->vint=IsVertexOnEdge;} 32 32 void Set(const Mesh &,long,Mesh &); 33 33 }; -
issm/trunk/src/c/objects/Bamg/VertexOnGeom.h
r5120 r5148 38 38 int OnGeomEdge() const {return this? abscisse >=0 :0;} 39 39 int IsRequiredVertex() {return this? ((abscisse<0 ? (gv?gv->Required():0):0 )) : 0;} 40 void SetOn(){mv-> onGeometry=this;mv->vint=IsVertexOnGeom;}40 void SetOn(){mv->GeometricalEdgeHook=this;mv->vint=IsVertexOnGeom;} 41 41 42 42 //Inline methods -
issm/trunk/src/c/objects/Bamg/VertexOnVertex.h
r5120 r5148 21 21 22 22 //Methods 23 void SetOnBTh(){v-> onBackgroundVertex=bv;v->vint=IsVertexOnVertex;}23 void SetOnBTh(){v->BackgroundVertexHook=bv;v->vint=IsVertexOnVertex;} 24 24 void Set(const Mesh &,long,Mesh &); 25 25 };
Note:
See TracChangeset
for help on using the changeset viewer.