![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2958
- Timestamp:
- 02/04/10 11:37:37 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 4 edited
-
Mesh2.h (modified) (7 diffs)
-
SetOfE4.h (modified) (1 diff)
-
objects/SetOfE4.cpp (modified) (2 diffs)
-
objects/Triangles.cpp (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.h
r2957 r2958 144 144 145 145 public: 146 Triangle* t; // le triangle147 int a; // le numero de l arete146 Triangle* t; //pointer toward the triangle 147 int a; //Edge number 148 148 149 149 TriangleAdjacent(Triangle * tt,int aa): t(tt),a(aa &3) {}; … … 153 153 operator Triangle & () const {return *t;} 154 154 operator int() const {return a;} 155 TriangleAdjacent & operator++() 156 { 155 TriangleAdjacent & operator++(){ 157 156 a= NextEdge[a]; 158 return *this;} 159 TriangleAdjacent operator--() 160 { 161 a= PreviousEdge[a]; 162 return *this;} 163 inline TriangleAdjacent Adj() const ; 164 int swap(); 165 inline void SetAdj2(const TriangleAdjacent& , int =0); 166 inline Vertex * EdgeVertex(const int &) const ; 167 inline Vertex * OppositeVertex() const ; 168 inline Icoor2 & det() const; 169 inline int Locked() const ; 170 inline int GetAllFlag_UnSwap() const ; 171 inline void SetLock(); 172 inline int MarkUnSwap() const; 173 inline void SetMarkUnSwap(); 174 inline void SetCracked(); 175 inline int Cracked() const ; 157 return *this; 158 } 159 TriangleAdjacent operator--(){ 160 a= PreviousEdge[a]; 161 return *this; 162 } 163 inline TriangleAdjacent Adj() const ; 164 int swap(); 165 inline void SetAdj2(const TriangleAdjacent& , int =0); 166 inline Vertex * EdgeVertex(const int &) const ; 167 inline Vertex * OppositeVertex() const ; 168 inline Icoor2 & det() const; 169 inline int Locked() const ; 170 inline int GetAllFlag_UnSwap() const ; 171 inline void SetLock(); 172 inline int MarkUnSwap() const; 173 inline void SetMarkUnSwap(); 174 inline void SetCracked(); 175 inline int Cracked() const ; 176 176 }; 177 177 /*}}}1*/ … … 298 298 Vertex* ns [3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer 299 299 Triangle* at [3]; // 3 adjacent triangles (nu) 300 Int1 aa[3]; // les nu des arete dans le triangles (mod 4)300 Int1 aa[3]; // les nu des arete dans le triangles (mod 3) 301 301 public: 302 302 Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres) … … 354 354 } 355 355 356 void SetAdj2(Int1 a,Triangle *t,Int1 aat) 357 { at[a]=t;aa[a]=aat; 358 if(t) {t->at[aat]=this;t->aa[aat]=a;} 359 } 356 void SetAdj2(Int1 a,Triangle *t,Int1 aat){ 357 //at = pointer toward the adjacent triangles of this 358 //aa = position of the edge in the triangle (mod 3) 359 at[a]=t; //the adjacent triangle to the edge a is t 360 aa[a]=aat; //position of the edge in the adjacent triangle 361 //if t!=NULL add adjacent triangle to t (this) 362 if(t) { 363 t->at[aat]=this; 364 t->aa[aat]=a; 365 } 366 } 360 367 361 368 void SetTriangleContainingTheVertex() … … 377 384 378 385 void SetHidden(int a){ 379 register Triangle * t = at[a]; 386 387 //Get Adjacent Triangle number a 388 register Triangle* t = at[a]; 389 390 //if it exist 391 //C|=D -> C=(C|D) bitwise inclusive OR 380 392 if(t) t->aa[aa[a] & 3] |=16; 381 393 aa[a] |= 16; … … 391 403 392 404 void SetLocked(int a){ 405 //mark the edge as on Boundary 393 406 register Triangle * t = at[a]; 394 407 t->aa[aa[a] & 3] |=4; … … 1071 1084 } 1072 1085 1073 inline void TriangleAdjacent::SetAdj2(const TriangleAdjacent & ta, int l ) 1074 { // set du triangle adjacent1086 inline void TriangleAdjacent::SetAdj2(const TriangleAdjacent & ta, int l ){ 1087 //set Adjacent Triangle of a triangle 1075 1088 if(t) { 1076 1089 t->at[a]=ta.t; 1077 t->aa[a]=ta.a|l;} 1078 if(ta.t) { 1079 ta.t->at[ta.a] = t ; 1080 ta.t->aa[ta.a] = a| l ; 1081 } 1082 } 1090 t->aa[a]=ta.a|l; 1091 } 1092 if(ta.t) { 1093 ta.t->at[ta.a] = t ; 1094 ta.t->aa[ta.a] = a| l ; 1095 } 1096 } 1083 1097 1084 1098 -
issm/trunk/src/c/Bamgx/SetOfE4.h
r2862 r2958 14 14 class SetOfEdges4 { 15 15 Int4 nx,nbax,NbOfEdges; 16 Int4 * tete;17 Int4Edge 16 Int4* head; 17 Int4Edge* Edges; 18 18 19 19 public: 20 20 SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet 21 ~SetOfEdges4() {delete [] tete; delete [] Edges;}21 ~SetOfEdges4() {delete [] head; delete [] Edges;} 22 22 Int4 add (Int4 ii,Int4 jj); 23 Int4 addtrie(Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj) : add (jj,ii) ;}23 Int4 SortAndAdd (Int4 ii,Int4 jj) {return ii <=jj ? add (ii,jj) : add (jj,ii) ;} 24 24 Int4 nb(){return NbOfEdges;} 25 25 Int4 find (Int4 ii,Int4 jj); 26 Int4 findtrie(Int4 ii,Int4 jj) {return ii <=jj ? find (ii,jj) : find (jj,ii) ;}26 Int4 SortAndFind (Int4 ii,Int4 jj) {return ii <=jj ? find (ii,jj) : find (jj,ii) ;} 27 27 Int4 i(Int4 k){return Edges[k].i;} 28 28 Int4 j(Int4 k){return Edges[k].j;} -
issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp
r2859 r2958 14 14 /*FUNCTION SetOfEdges4::SetOfEdges4(Int4 mmx,Int4 nnx){{{1*/ 15 15 SetOfEdges4::SetOfEdges4(Int4 mmx,Int4 nnx){ 16 nx=nnx; 17 nbax=mmx; 18 NbOfEdges = 0; 19 tete= new Int4 [nx]; 20 Int4 i=nx; 21 while (i--) 22 tete[i] = -1;// vide 23 Edges =new Int4Edge[nbax]; 16 17 /*Intermediary*/ 18 int i; 19 20 //initialize fields 21 nx =nnx; //number of vertices 22 nbax =mmx; // 3 * number of triangles 23 NbOfEdges=0; 24 head = new Int4 [nx]; 25 Edges= new Int4Edge[nbax]; 26 27 //initialize head (-1 everywhere) 28 i=nx; 29 while (i--) head[i]=-1; 24 30 } 25 31 /*}}}1*/ … … 28 34 /*FUNCTION SetOfEdges4::find {{{1*/ 29 35 Int4 SetOfEdges4::find(Int4 ii,Int4 jj) { 30 if (tete == 0 ) { 31 throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::find no more tete de liste (?)")); 36 37 /*Intermediary*/ 38 int n; 39 40 //check that head is not empty 41 if (head==NULL) { 42 throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add no head is NULL (How come?)")); 32 43 } 33 Int4 n = tete[ Abs( ii ) % nx ];34 44 35 while (n >= 0) 36 if (ii == Edges[n].i && jj == Edges[n].j) return n; 37 else n = Edges[n].next; 38 return -1; //do not exist, return -1 45 //get n from h (usually h=ii) 46 n=head[Abs(ii)%nx]; 47 48 //go through the existing edges that holds h (=ii) and return position in Edge 49 while (n >= 0){ 50 51 //if the edge ii jj is already in Edges, return n 52 if (ii == Edges[n].i && jj == Edges[n].j) return n; 53 54 //else go to next edge that holds ii 55 else n = Edges[n].next; 56 } 57 58 //if we reach this point, the edge does not exist return -1 59 return -1; 39 60 } 40 61 /*}}}1*/ 41 62 /*FUNCTION SetOfEdges4::add{{{1*/ 42 63 Int4 SetOfEdges4::add(Int4 ii,Int4 jj) { 43 if (tete == 0 ) { 44 throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add no more tete de liste (?)")); 64 65 /*Intermediary*/ 66 int h,n; 67 68 //check that head is not empty 69 if (head==NULL) { 70 throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add no head is NULL (How come?)")); 45 71 } 46 Int4 h; 47 Int4 n = tete[ h = Abs( ii ) % nx ]; 48 while (n >= 0) 49 if (ii == Edges[n].i && jj == Edges[n].j) 50 return n; 51 else n = Edges[n].next; 72 73 //get n from h (usually h=ii) 74 n=head[h=Abs(ii)%nx]; 75 76 //go through the existing edges that holds h (=ii) and check that 77 //the edge ii jj is not already in Edge 78 while (n >= 0){ 79 80 //if the edge ii jj is already in Edges, return n 81 if (ii == Edges[n].i && jj == Edges[n].j) return n; 82 83 //else go to next edge that holds ii 84 else n = Edges[n].next; 85 } 86 87 //check that nbax <=NbOfEdges 52 88 if (nbax <=NbOfEdges ) { 53 89 throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add overflow: NbOfEdges=%i > nbax=%i",NbOfEdges,nbax)); 54 90 } 55 91 92 //update chain 56 93 Edges[NbOfEdges].i=ii; 57 94 Edges[NbOfEdges].j=jj; 58 Edges[NbOfEdges].next= tete[h];59 tete[h] = NbOfEdges;95 Edges[NbOfEdges].next= head[h]; 96 head[h] = NbOfEdges; 60 97 return NbOfEdges ++; 61 98 } 62 99 /*}}}1*/ 63 64 100 } -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2957 r2958 881 881 882 882 /*Intermediary*/ 883 int i,j,k ;883 int i,j,k,kk,it,jt; 884 884 885 885 /*Recover options*/ … … 899 899 if (cutoffradian>=0) Gh.MaximalAngleOfCorner = cutoffradian; 900 900 901 //Construction of the edges 901 /*Construction of the edges*/ 902 903 //initialize st and edge4 902 904 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv); 903 905 Int4* st = new Int4[nbt*3]; 904 for (i=0;i<nbt*3;i++) 905 st[i]=-1;906 Int4 kk =0;907 908 Int4 nbeold = nbe;906 907 //initialize st as -1 (chaining algorithm) 908 for (i=0;i<nbt*3;i++) st[i]=-1; 909 910 //build edge4 (chain) 909 911 for (i=0;i<nbe;i++){ 910 edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])); 911 } 912 edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])); 913 } 914 //check that there is no double edge 912 915 if (nbe != edge4->nb()){ 913 916 throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i, nbe4=%i",nbe,edge4->nb())); 914 917 } 918 //keep nbe in nbeold 919 Int4 nbeold = nbe; 920 921 //Go through the triangles and ass the edges in edge4 if they are not there yet 915 922 for (i=0;i<nbt;i++){ 923 //3 edges per triangle 916 924 for (j=0;j<3;j++) { 917 Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 925 //Add Edge to edge4 (k=numberofedges in edge4) 926 Int4 k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 918 927 Int4 invisible = triangles[i].Hidden(j); 928 929 //if st[k] has not been changed yet, add 3*i+j (= vertex position in the index) 919 930 if(st[k]==-1) st[k]=3*i+j; 931 932 //else st[k]>=0 -> the edge already exist, check 920 933 else if(st[k]>=0) { 934 //check that it is not an edge on boundary (should not alrady exist) 921 935 if (triangles[i].TriangleAdj(j) || triangles[st[k]/3].TriangleAdj((int) (st[k]%3))){ 922 throw ErrorException(__FUNCT__,exprintf(" triangles[i].TriangleAdj(j) || triangles[st[k]/3].TriangleAdj((int) (st[k]%3))"));936 throw ErrorException(__FUNCT__,exprintf("problem in Geometry reconstruction: an edge on boundary is duplicated (double element?)")); 923 937 } 938 //OK, the element is not on boundary, is belongs to 2 triangles -> build Adjacent triangles list 924 939 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); 925 940 if (invisible) triangles[i].SetHidden(j); 941 // if k < nbe mark the adge as on Boundary (Locked) 926 942 if (k<nbe) { 927 943 triangles[i].SetLocked(j); 928 944 } 945 //set st[k] as negative so that the edge should not be called again 929 946 st[k]=-2-st[k]; 930 947 } 948 //else (see 3 lines above), the edge has been called more than twice: return error 931 949 else { 932 950 printf("The edge (%i,%i) belongs to more than 2 triangles (%i)\n",Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]),k); … … 938 956 } 939 957 } 958 959 //delete edge4 940 960 Int4 nbedges = edge4->nb(); // the total number of edges 941 delete edge4; 942 edge4 =0; 943 961 delete edge4; edge4=NULL; 962 963 //display info 944 964 if(verbosity>5) { 945 printf(" info on Mesh s:\n");965 printf(" info on Mesh:\n"); 946 966 printf(" - number of vertices = %i \n",nbv); 947 967 printf(" - number of triangles = %i \n",nbt); … … 950 970 printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv); 951 971 } 952 // check the consistant of edge[].adj and the geometrical required vertex 953 k=0; 954 kk=0; 955 Int4 it; 956 957 for (i=0;i<nbedges;i++) 958 if (st[i] <-1) {// edge internal 959 it = (-2-st[i])/3; 960 j = (int) ((-2-st[i])%3); 961 Triangle & tt = * triangles[it].TriangleAdj(j); 962 if (triangles[it].color != tt.color|| i < nbeold) k++; 963 } 964 else if (st[i] >=0) // edge alone 965 kk++; 972 973 // check consistency of edge[].adj and geometrical required vertices 974 k=0; kk=0; 975 for (i=0;i<nbedges;i++){ 976 //internal edge 977 if (st[i] <-1) { 978 //get triangle number back 979 it = (-2-st[i])/3; 980 //get edge position back 981 j = (int) ((-2-st[i])%3); 982 Triangle &tt=*triangles[it].TriangleAdj(j); 983 if (triangles[it].color != tt.color|| i < nbeold) k++; 984 } 985 //boundary edge (alone) 986 else if (st[i] >=0) 987 kk++; 988 } 989 990 /*Constructions of edges*/ 966 991 967 992 k += kk; 968 993 kk=0; 969 994 if (k) { 970 // construction of the edges971 995 nbe = k; 972 Edge * edgessave =edges;996 Edge* edgessave=edges; 973 997 edges = new Edge[nbe]; 974 998 k =0; 975 // construction of the edges 999 1000 //display info 976 1001 if(verbosity>4) printf(" Construction of the edges %i\n",nbe); 977 1002 … … 979 1004 Int4 add= -1; 980 1005 981 if (st[i] <-1) // edge internal982 1006 //internal edge (belongs to two triangles) 1007 if (st[i] <-1){ 983 1008 it = (-2-st[i])/3; 984 1009 j = (int) ((-2-st[i])%3); 985 1010 Triangle & tt = * triangles[it].TriangleAdj(j); 986 if (triangles[it].color != tt.color || i < nbeold) // Modif FH 061220551011 if (triangles[it].color != tt.color || i < nbeold) 987 1012 add=k++; 988 989 else if (st[i] >=0) // edge alone990 1013 } 1014 //boundary edge 1015 else if (st[i] >=0){ 991 1016 it = st[i]/3; 992 1017 j = (int) (st[i]%3); 993 1018 add=k++; 994 } 995 996 if (add>=0 && add < nbe) 997 { 998 1019 } 1020 if (add>=0 && add < nbe){ 999 1021 edges[add].v[0] = &triangles[it][VerticesOfTriangularEdge[j][0]]; 1000 1022 edges[add].v[1] = &triangles[it][VerticesOfTriangularEdge[j][1]]; 1001 1023 edges[add].onGeometry=0; 1002 if (i<nbeold) // in file edge // Modif FH 061220551003 1004 edges[add].ref =edgessave[i].ref;1005 edges[add].onGeometry = edgessave[i].onGeometry; // HACK pour recuperer les aretes requise midf FH avril 2006 ????1006 1024 //if already existed 1025 if (i<nbeold){ 1026 edges[add].ref=edgessave[i].ref; 1027 edges[add].onGeometry=edgessave[i].onGeometry; // HACK to get required edges 1028 } 1007 1029 else 1008 edges[add].ref = Min(edges[add].v[0]->ref(),edges[add].v[1]->ref()); // no a good choice1030 edges[add].ref=Min(edges[add].v[0]->ref(),edges[add].v[1]->ref()); 1009 1031 } 1010 1032 } 1033 1034 //check that we have been through all edges 1011 1035 if (k!=nbe){ 1012 throw ErrorException(__FUNCT__,exprintf("k!=nbe")); 1013 } 1036 throw ErrorException(__FUNCT__,exprintf("problem in edge construction process: k!=nbe (should not happen)")); 1037 } 1038 //delete edgessave 1014 1039 if (edgessave) delete [] edgessave; 1015 1040 } 1016 1041 1017 // construction of edges[].adj 1018 for (i=0;i<nbv;i++) 1019 vertices[i].color =0; 1020 for (i=0;i<nbe;i++) 1021 for (j=0;j<2;j++) 1022 edges[i].v[j]->color++; 1023 1024 for (i=0;i<nbv;i++) 1025 vertices[i].color = (vertices[i].color ==2) ? -1 : -2; 1026 for (i=0;i<nbe;i++) 1027 for (j=0;j<2;j++) 1028 { 1029 Vertex *v=edges[i].v[j]; 1030 Int4 i0=v->color,j0; 1031 if(i0<0) 1032 edges[i ].adj[ j ]=0; // Add FH Jan 2008 1033 if(i0==-1) 1034 v->color=i*2+j; 1035 else if (i0>=0) {// i and i0 edge are adjacent by the vertex v 1036 j0 = i0%2; 1037 i0 = i0/2; 1038 if (v!=edges[i0 ].v[j0]){ 1039 throw ErrorException(__FUNCT__,exprintf("v!=edges[i0 ].v[j0]")); 1040 } 1041 edges[i ].adj[ j ] =edges +i0; 1042 edges[i0].adj[ j0] =edges +i ; 1043 v->color = -3;} 1044 } 1045 // now reconstruct the sub domain info 1042 /*Color the vertices*/ 1043 1044 //initialize color of all vertices as 0 1045 for (i=0;i<nbv;i++) vertices[i].color =0; 1046 1047 //go through the edges and add a color to corresponding vertices 1048 //(A vertex in 4 edges will have a color 4) 1049 for (i=0;i<nbe;i++){ 1050 for (j=0;j<2;j++) edges[i].v[j]->color++; 1051 } 1052 1053 //change the color: if a vertex belongs to 2 edges -1, else -2 1054 for (i=0;i<nbv;i++) { 1055 vertices[i].color=(vertices[i].color ==2)? -1 : -2; 1056 } 1057 1058 /*Build edges[].adj*/ 1059 1060 for (i=0;i<nbe;i++){ 1061 for (j=0;j<2;j++){ 1062 //get current vertex 1063 Vertex* v=edges[i].v[j]; 1064 //get vertex color (i0) 1065 Int4 i0=v->color,j0; 1066 1067 //if color<0 (first time), no adjacent edge 1068 if(i0<0) edges[i].adj[j]=NULL; 1069 1070 //if color=-1 (corner),change the vertex color as 2*i+j (position of the vertex in edges) 1071 if(i0==-1) v->color=i*2+j; 1072 1073 //if color>=0 (i and i0 edge are adjacent by the vertex v) 1074 else if (i0>=0) { 1075 //get position of v in edges back 1076 j0 = i0%2; //column in edges 1077 i0 = i0/2; //line in edges 1078 1079 //check that we have the correct vertex 1080 if (v!=edges[i0 ].v[j0]){ 1081 throw ErrorException(__FUNCT__,exprintf("v!=edges[i0 ].v[j0]: this should not happen as the vertex belongs to this edge")); 1082 } 1083 1084 //Add adjacence 1085 edges[i ].adj[j ]=edges +i0; 1086 edges[i0].adj[j0]=edges +i ; 1087 1088 //change color to -3 1089 v->color = -3; 1090 } 1091 } 1092 } 1093 1094 /*Reconstruct subdomains info*/ 1095 1096 //check that NbSubDomains is empty 1046 1097 if (NbSubDomains){ 1047 1098 throw ErrorException(__FUNCT__,exprintf("NbSubDomains should be 0")); … … 1049 1100 NbSubDomains=0; 1050 1101 1051 { 1052 Int4 it; 1053 // find all the sub domain 1054 Int4 *colorT = new Int4[nbt]; 1055 Triangle *tt,*t; 1056 Int4 k; 1057 for ( it=0;it<nbt;it++) 1058 colorT[it]=-1; 1059 for (it=0;it<nbt;it++) 1060 { 1061 if (colorT[it]<0) 1062 { 1063 colorT[it]=NbSubDomains; 1064 Int4 level =1,j,jt,kolor=triangles[it].color; 1065 st[0]=it; // stack 1066 st[1]=0; 1067 k=1; 1068 while (level>0) 1069 if( ( j=st[level]++) <3) 1070 { 1071 t = &triangles[st[level-1]]; 1072 tt=t->TriangleAdj((int)j); 1073 1074 if ( ! t->Locked(j) && tt && (colorT[jt = Number(tt)] == -1) && ( tt->color==kolor)) 1075 { 1076 colorT[jt]=NbSubDomains; 1077 st[++level]=jt; 1078 st[++level]=0; 1079 k++; 1080 } 1102 //color the subdomains 1103 Int4* colorT= new Int4[nbt]; 1104 Triangle *tt,*t; 1105 1106 //initialize the color of each triangle as -1 1107 for (it=0;it<nbt;it++) colorT[it]=-1; 1108 1109 //loop over the triangles 1110 for (it=0;it<nbt;it++){ 1111 1112 //if the triangle has not been colored yet: 1113 if (colorT[it]<0){ 1114 1115 //color = number of subdomains 1116 colorT[it]=NbSubDomains; 1117 1118 //color all the adjacent triangles of T that share a non marked edge 1119 int level =1; 1120 int kolor=triangles[it].color; 1121 st[0]=it; // stack 1122 st[1]=0; 1123 k=1; 1124 while (level>0){ 1125 if( (j=st[level]++)<3 ){ 1126 t = &triangles[st[level-1]]; 1127 tt=t->TriangleAdj((int)j); 1128 1129 //color the adjacent triangle 1130 if ( ! t->Locked(j) && tt && (colorT[jt = Number(tt)] == -1) && ( tt->color==kolor)) { 1131 colorT[jt]=NbSubDomains; 1132 st[++level]=jt; 1133 st[++level]=0; 1134 k++; 1081 1135 } 1082 else1083 level-=2;1084 NbSubDomains++;1085 }1086 }1087 if (verbosity> 3) printf(" The Number of sub domain = %i\n",NbSubDomains);1088 1089 Int4 isd;1090 subdomains = new SubDomain[NbSubDomains];1091 for (isd=0;isd<NbSubDomains;isd++)1092 {1093 subdomains[isd].head =0;1094 }1095 k=0;1096 for (it=0;it<nbt;it++)1097 for (int j=0;j<3;j++)1098 {1099 tt=triangles[it].TriangleAdj(j);1100 if ((!tt || tt->color != triangles[it].color) && !subdomains[isd=colorT[it]].head)1101 {1102 subdomains[isd].head = triangles+it;1103 subdomains[isd].ref = triangles[it].color;1104 subdomains[isd].sens = j; // hack1105 subdomains[isd].edge = 0;1106 k++;1107 1136 } 1108 } 1109 if (k!= NbSubDomains){ 1110 throw ErrorException(__FUNCT__,exprintf("k!= NbSubDomains")); 1111 } 1112 1113 delete [] colorT; 1114 1115 1116 } 1137 else level-=2; 1138 } 1139 NbSubDomains++; 1140 } 1141 } 1142 if (verbosity> 3) printf(" The Number of sub domain = %i\n",NbSubDomains); 1143 1144 //build subdomains 1145 Int4 isd; 1146 subdomains = new SubDomain[NbSubDomains]; 1147 1148 //initialize subdomains[isd].head as 0 1149 for (isd=0;isd<NbSubDomains;isd++) subdomains[isd].head =0; 1150 1151 k=0; 1152 for (it=0;it<nbt;it++){ 1153 for (int j=0;j<3;j++){ 1154 tt=triangles[it].TriangleAdj(j); 1155 if ((!tt || tt->color != triangles[it].color) && !subdomains[isd=colorT[it]].head){ 1156 subdomains[isd].head = triangles+it; 1157 subdomains[isd].ref = triangles[it].color; 1158 subdomains[isd].sens = j; // hack 1159 subdomains[isd].edge = 0; 1160 k++; 1161 } 1162 } 1163 } 1164 //check that we have been through all subdomains 1165 if (k!= NbSubDomains){ 1166 throw ErrorException(__FUNCT__,exprintf("k!= NbSubDomains")); 1167 } 1168 //delete colorT and st 1169 delete [] colorT; 1117 1170 delete [] st; 1118 // now make the geometry 1119 // 1 compress the vertices 1120 Int4 * colorV = new Int4[nbv]; 1121 for (i=0;i<nbv;i++) 1122 colorV[i]=-1; 1123 for (i=0;i<nbe;i++) 1124 for ( j=0;j<2;j++) 1125 colorV[Number(edges[i][j])]=0; 1171 1172 /*Reconstruct Geometry Gh*/ 1173 1174 //build colorV -1 for all vertex and 0 for the vertices belonging to edges 1175 Int4* colorV = new Int4[nbv]; 1176 for (i=0;i<nbv;i++) colorV[i]=-1; 1177 for (i=0;i<nbe;i++){ 1178 for ( j=0;j<2;j++) colorV[Number(edges[i][j])]=0; 1179 } 1180 //number the vertices nelonging to edges 1126 1181 k=0; 1127 for (i=0;i<nbv;i++) 1128 if(!colorV[i]) 1129 colorV[i]=k++; 1130 1182 for (i=0;i<nbv;i++){ 1183 if(!colorV[i]) colorV[i]=k++; 1184 } 1185 1186 //Build Gh 1131 1187 Gh.nbv=k; 1132 1188 Gh.nbe = nbe; … … 1140 1196 NbVerticesOnGeomEdge =0; 1141 1197 VerticesOnGeomEdge =0; 1142 { 1143 Int4 j; 1144 for (i=0;i<nbv;i++) 1145 if((j=colorV[i])>=0) 1146 { 1147 1148 Vertex & v = Gh.vertices[j]; 1149 v = vertices[i]; 1150 v.color =0; 1151 VerticesOnGeomVertex[j] = VertexOnGeom(vertices[i], Gh.vertices[j]); 1152 } 1153 1154 } 1155 edge4= new SetOfEdges4(nbe,nbv); 1156 1157 Real4 * len = new Real4[Gh.nbv]; 1158 for(i=0;i<Gh.nbv;i++) 1159 len[i]=0; 1160 1198 1199 //Build VertexOnGeom 1200 for (i=0;i<nbv;i++){ 1201 if((j=colorV[i])>=0){ 1202 Vertex & v = Gh.vertices[j]; 1203 v = vertices[i]; 1204 v.color =0; 1205 VerticesOnGeomVertex[j] = VertexOnGeom(vertices[i], Gh.vertices[j]); 1206 } 1207 } 1208 1209 //Buid pmin and pmax of Gh (extrema coordinates) 1161 1210 Gh.pmin = Gh.vertices[0].r; 1162 1211 Gh.pmax = Gh.vertices[0].r; … … 1168 1217 Gh.pmax.y = Max(Gh.pmax.y,Gh.vertices[i].r.y); 1169 1218 } 1170 1171 1219 R2 DD05 = (Gh.pmax-Gh.pmin)*0.05; 1172 1220 Gh.pmin -= DD05; 1173 1221 Gh.pmax += DD05; 1174 1222 1223 //Build Gh.coefIcoor 1175 1224 Gh.coefIcoor= (MaxICoor)/(Max(Gh.pmax.x-Gh.pmin.x,Gh.pmax.y-Gh.pmin.y)); 1176 1225 if (Gh.coefIcoor<=0){ 1177 throw ErrorException(__FUNCT__,exprintf("Gh.coefIcoor<=0")); 1178 } 1179 1226 throw ErrorException(__FUNCT__,exprintf("Gh.coefIcoor<=0 in infered Geometry (this should not happen)")); 1227 } 1228 1229 /*Build Gh.edges*/ 1230 1231 //initialize len as 0 1232 Real4 * len = new Real4[Gh.nbv]; 1233 for(i=0;i<Gh.nbv;i++) len[i]=0; 1234 1235 //initialize edge4 again 1236 edge4= new SetOfEdges4(nbe,nbv); 1180 1237 Real8 hmin = HUGE_VAL; 1181 1238 int kreq=0; 1182 for (i=0;i<nbe;i++) 1183 { 1239 for (i=0;i<nbe;i++){ 1240 1184 1241 Int4 i0 = Number(edges[i][0]); 1185 1242 Int4 i1 = Number(edges[i][1]); 1186 Int4 j0 = 1187 Int4 j1 = 1243 Int4 j0 = colorV[i0]; 1244 Int4 j1 = colorV[i1]; 1188 1245 1189 1246 Gh.edges[i].v[0] = Gh.vertices + j0; 1190 1247 Gh.edges[i].v[1] = Gh.vertices + j1; 1248 1191 1249 Gh.edges[i].flag = 0; 1250 1192 1251 Gh.edges[i].tg[0]=R2(); 1193 1252 Gh.edges[i].tg[1]=R2(); 1194 bool requis= edges[i].onGeometry; 1195 if(requis) kreq++; 1253 1254 bool required= edges[i].onGeometry; 1255 if(required) kreq++; 1196 1256 edges[i].onGeometry = Gh.edges + i; 1197 if(requi s) { // correction fevr 2009 JYU ...1257 if(required){ 1198 1258 Gh.edges[i].v[0]->SetRequired(); 1199 1259 Gh.edges[i].v[1]->SetRequired(); 1200 Gh.edges[i].SetRequired(); // fin modif ... 1201 } 1260 Gh.edges[i].SetRequired(); 1261 } 1262 1202 1263 R2 x12 = Gh.vertices[j0].r-Gh.vertices[j1].r; 1203 1264 Real8 l12=Norme2(x12); … … 1212 1273 Gh.edges[i].ref = edges[i].ref; 1213 1274 1214 k = edge4-> addtrie(i0,i1);1275 k = edge4->SortAndAdd(i0,i1); 1215 1276 if (k != i){ 1216 throw ErrorException(__FUNCT__,exprintf("k != i")); 1217 } 1218 1219 } 1220 1221 1222 for (i=0;i<Gh.nbv;i++) 1277 throw ErrorException(__FUNCT__,exprintf("problem in Edge4 construction: k != i")); 1278 } 1279 } 1280 1281 //Build metric for all vertices of Gh 1282 for (i=0;i<Gh.nbv;i++){ 1223 1283 if (Gh.vertices[i].color > 0) 1224 1284 Gh.vertices[i].m= Metric(len[i] /(Real4) Gh.vertices[i].color); 1225 1285 else 1226 1286 Gh.vertices[i].m= Metric(hmin); 1287 } 1288 //delete len 1227 1289 delete [] len; 1228 for (i=0;i<NbSubDomains;i++) 1229 { 1230 Int4 it = Number(subdomains[i].head); 1231 int j = subdomains[i].sens; 1290 1291 //Build Gh.subdomains 1292 for (i=0;i<NbSubDomains;i++){ 1293 it = Number(subdomains[i].head); 1294 j = subdomains[i].sens; 1232 1295 Int4 i0 = Number(triangles[it][VerticesOfTriangularEdge[j][0]]); 1233 1296 Int4 i1 = Number(triangles[it][VerticesOfTriangularEdge[j][1]]); 1234 k = edge4->findtrie(i0,i1); 1235 if(k>=0) 1236 { 1297 k = edge4->SortAndFind(i0,i1); 1298 if(k>=0){ 1237 1299 subdomains[i].sens = (vertices + i0 == edges[k].v[0]) ? 1 : -1; 1238 1300 subdomains[i].edge = edges+k; … … 1240 1302 Gh.subdomains[i].sens = subdomains[i].sens; 1241 1303 Gh.subdomains[i].ref = subdomains[i].ref; 1242 1304 } 1243 1305 else 1244 1306 throw ErrorException(__FUNCT__,exprintf("%i should be >=0")); … … 1247 1309 delete edge4; 1248 1310 delete [] colorV; 1249 // -- unset adj 1250 for (i=0;i<nbt;i++) 1251 for ( j=0;j<3;j++) 1252 triangles[i].SetAdj2(j,0,triangles[i].GetAllflag(j)); 1311 1312 //unset adj 1313 for (i=0;i<nbt;i++){ 1314 for ( j=0;j<3;j++){ 1315 triangles[i].SetAdj2(j,0,triangles[i].GetAllflag(j)); 1316 } 1317 } 1253 1318 1254 1319 } … … 2006 2071 Int4 kk =0; 2007 2072 for (i=0;i<nbe;i++){ 2008 kk=kk+(i == edge4-> addtrie(Number(edges[i][0]),Number(edges[i][1])));2073 kk=kk+(i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1]))); 2009 2074 } 2010 2075 if (kk != nbe) { … … 2015 2080 for (i=0;i<nbt;i++){ 2016 2081 for (int j=0;j<3;j++) { 2017 Int4 k =edge4-> addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),2082 Int4 k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), 2018 2083 Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 2019 2084 Int4 invisible = triangles[i].Hidden(j); … … 2177 2242 Vertex *v0= ta.EdgeVertex(0); 2178 2243 Vertex *v1= ta.EdgeVertex(1); 2179 Int4 k =edge4-> addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv);2244 Int4 k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv); 2180 2245 if (st[k]<0){ 2181 2246 throw ErrorException(__FUNCT__,exprintf("st[k]<0")); … … 4504 4569 4505 4570 4506 kk += (i == edge4-> addtrie(Number(edges[i][0]),Number(edges[i][1])));4571 kk += (i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1]))); 4507 4572 if (ong) // a geometrical edges 4508 4573 { … … 4595 4660 const Vertex & v0 = t[VerticesOfTriangularEdge[j][0]]; 4596 4661 const Vertex & v1 = t[VerticesOfTriangularEdge[j][1]]; 4597 Int4 ke =edge4-> findtrie(Number(v0),Number(v1));4662 Int4 ke =edge4->SortAndFind(Number(v0),Number(v1)); 4598 4663 if (ke>0) 4599 4664 { … … 4646 4711 if ( kedge[3*i+j] < 0) 4647 4712 { 4648 Int4 ke =edge4-> findtrie(Number(v0),Number(v1));4713 Int4 ke =edge4->SortAndFind(Number(v0),Number(v1)); 4649 4714 if (ke<0) // new 4650 4715 {
Note:
See TracChangeset
for help on using the changeset viewer.
