![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 21623
- Timestamp:
- 03/23/17 14:07:02 (8 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 3 deleted
- 17 edited
-
Makefile.am (modified) (2 diffs)
-
bamg/BamgMesh.cpp (modified) (2 diffs)
-
bamg/BamgMesh.h (modified) (1 diff)
-
bamg/BamgQuadtree.cpp (modified) (15 diffs)
-
bamg/BamgQuadtree.h (modified) (4 diffs)
-
bamg/BamgVertex.cpp (modified) (3 diffs)
-
bamg/BamgVertex.h (modified) (4 diffs)
-
bamg/Direction.cpp (deleted)
-
bamg/Direction.h (deleted)
-
bamg/DoubleAndInt.h (deleted)
-
bamg/EigenMetric.cpp (modified) (5 diffs)
-
bamg/Geometry.cpp (modified) (4 diffs)
-
bamg/Mesh.cpp (modified) (26 diffs)
-
bamg/Mesh.h (modified) (2 diffs)
-
bamg/Metric.cpp (modified) (5 diffs)
-
bamg/Metric.h (modified) (5 diffs)
-
bamg/Triangle.cpp (modified) (6 diffs)
-
bamg/Triangle.h (modified) (1 diff)
-
bamg/bamgobjects.h (modified) (1 diff)
-
modules/Bamgx/Bamgx.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/Makefile.am
r21583 r21623 20 20 21 21 #Core sources 22 #BAMG sources {{{ 23 issm_sources = 24 if BAMG 25 issm_sources += ./bamg/BamgGeom.cpp\ 26 ./bamg/BamgMesh.cpp\ 27 ./bamg/BamgOpts.cpp\ 28 ./bamg/CrackedEdge.cpp\ 29 ./bamg/Curve.cpp\ 30 ./bamg/Edge.cpp\ 31 ./bamg/GeomEdge.cpp\ 32 ./bamg/GeomSubDomain.cpp\ 33 ./bamg/GeomVertex.cpp\ 34 ./bamg/Geometry.cpp\ 35 ./bamg/ListofIntersectionTriangles.cpp\ 36 ./bamg/EigenMetric.cpp\ 37 ./bamg/Metric.cpp\ 38 ./bamg/BamgQuadtree.cpp\ 39 ./bamg/SetOfE4.cpp\ 40 ./bamg/SubDomain.cpp\ 41 ./bamg/AdjacentTriangle.cpp\ 42 ./bamg/Triangle.cpp\ 43 ./bamg/BamgVertex.cpp\ 44 ./bamg/VertexOnEdge.cpp\ 45 ./bamg/VertexOnGeom.cpp\ 46 ./bamg/VertexOnVertex.cpp\ 47 ./bamg/Mesh.cpp\ 48 ./shared/Bamg/BigPrimeNumber.cpp\ 49 ./modules/Bamgx/Bamgx.cpp\ 50 ./modules/BamgConvertMeshx/BamgConvertMeshx.cpp\ 51 ./modules/BamgTriangulatex/BamgTriangulatex.cpp 52 endif 53 #}}} 22 54 #Core sources{{{ 23 issm_sources = ./datastructures/DataSet.cpp\55 issm_sources += ./datastructures/DataSet.cpp\ 24 56 ./classes/gauss/GaussSeg.cpp\ 25 57 ./classes/gauss/GaussTria.cpp\ … … 282 314 endif 283 315 #}}} 284 #BAMG sources {{{285 if BAMG286 issm_sources += ./bamg/BamgGeom.cpp\287 ./bamg/BamgMesh.cpp\288 ./bamg/BamgOpts.cpp\289 ./bamg/CrackedEdge.cpp\290 ./bamg/Curve.cpp\291 ./bamg/Direction.cpp\292 ./bamg/Edge.cpp\293 ./bamg/GeomEdge.cpp\294 ./bamg/GeomSubDomain.cpp\295 ./bamg/GeomVertex.cpp\296 ./bamg/Geometry.cpp\297 ./bamg/ListofIntersectionTriangles.cpp\298 ./bamg/EigenMetric.cpp\299 ./bamg/Metric.cpp\300 ./bamg/BamgQuadtree.cpp\301 ./bamg/SetOfE4.cpp\302 ./bamg/SubDomain.cpp\303 ./bamg/AdjacentTriangle.cpp\304 ./bamg/Triangle.cpp\305 ./bamg/BamgVertex.cpp\306 ./bamg/VertexOnEdge.cpp\307 ./bamg/VertexOnGeom.cpp\308 ./bamg/VertexOnVertex.cpp\309 ./bamg/Mesh.cpp\310 ./shared/Bamg/BigPrimeNumber.cpp\311 ./modules/Bamgx/Bamgx.cpp\312 ./modules/BamgConvertMeshx/BamgConvertMeshx.cpp\313 ./modules/BamgTriangulatex/BamgTriangulatex.cpp314 endif315 #}}}316 316 #Petsc sources {{{ 317 317 if PETSC -
issm/trunk-jpl/src/c/bamg/BamgMesh.cpp
r18455 r21623 8 8 this->EdgesSize[0]=0, this->EdgesSize[1]=0; this->Edges=NULL; 9 9 this->TrianglesSize[0]=0, this->TrianglesSize[1]=0; this->Triangles=NULL; 10 this->QuadrilateralsSize[0]=0, this->QuadrilateralsSize[1]=0; this->Quadrilaterals=NULL;11 10 12 11 this->SubDomainsSize[0]=0, this->SubDomainsSize[1]=0; this->SubDomains=NULL; … … 33 32 xDelete<double>(this->Edges); 34 33 xDelete<double>(this->Triangles); 35 xDelete<double>(this->Quadrilaterals);36 34 37 35 xDelete<double>(this->SubDomains); -
issm/trunk-jpl/src/c/bamg/BamgMesh.h
r18455 r21623 16 16 int TrianglesSize[2]; 17 17 double* Triangles; 18 int QuadrilateralsSize[2];19 double* Quadrilaterals;20 18 21 19 int VerticesOnGeomVertexSize[2]; -
issm/trunk-jpl/src/c/bamg/BamgQuadtree.cpp
r18064 r21623 7 7 8 8 namespace bamg { 9 10 /*MACROS {{{*/11 /*12 *13 * J j14 * ^ ^15 * | | +--------+--------+16 * | | | | |17 * 1X | | | 2 | 3 |18 * | | | | |19 * | | +--------+--------+20 * | | | | |21 * 0X | | | 0 | 1 |22 * | | | | |23 * | | +--------+--------+24 * | +-----------------------> i25 * |26 * |----------------------------> I27 * X0 X128 *29 * box 0 -> I=0 J=0 IJ=00 = 030 * box 1 -> I=1 J=0 IJ=01 = 131 * box 2 -> I=0 J=1 IJ=10 = 232 * box 3 -> I=1 J=1 IJ=11 = 333 */34 #define INTER_SEG(a,b,x,y) (((y) > (a)) && ((x) <(b)))35 9 #define ABS(i) ((i)<0 ?-(i) :(i)) 36 #define MAX1(i,j) ((i)>(j) ?(i) :(j))37 #define NORM(i1,j1,i2,j2) MAX1(ABS((i1)-(j1)),ABS((i2)-(j2)))38 39 //IJ(i,j,l) returns the box number of i and j with respect to l40 //if !j&l and !i&l -> 0 (box zero: lower left )41 //if !j&l and i&l -> 1 (box one: lower right)42 //if j&l and !i&l -> 2 (box two: upper left )43 //if j&l and i&l -> 3 (box three:upper right)44 #define IJ(i,j,l) ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 ))45 46 //I_IJ(k,l) returns l if first bit of k is 1, else 047 #define I_IJ(k,l) ((k&1) ? l:0)48 //J_IJ(k,l) returns l if second bit of k is 1, else 049 #define J_IJ(k,l) ((k&2) ? l:0)50 /*}}}*/51 10 /*DOCUMENTATION What is a BamgQuadtree? {{{ 52 11 * A Quadtree is a very simple way to group vertices according … … 100 59 BamgQuadtree::BamgQuadtree(){/*{{{*/ 101 60 102 /*Number of boxes and vertices*/ 103 NbBamgQuadtreeBox=0; 104 NbVertices=0; 61 /*Initialize fields*/ 62 this->MaxDepth = 30; 63 this->MaxISize = 1073741824; //2^30 64 this->MaxICoord = 1073741823; //2^30 - 1 65 this->NbQuadtreeBox = 0; 66 this->NbVertices = 0; 105 67 106 68 /*Create container*/ 107 boxcontainer=new DataSet();69 this->boxcontainer=new DataSet(); 108 70 109 71 /*Create Root, pointer toward the main box*/ 110 root=NewBamgQuadtreeBox(); 111 112 } 113 /*}}}*/ 72 this->root=NewBamgQuadtreeBox(); 73 74 } /*}}}*/ 114 75 BamgQuadtree::BamgQuadtree(Mesh * t,long nbv){ /*{{{*/ 115 76 116 /*Number of boxes and vertices*/ 117 NbBamgQuadtreeBox=0; 118 NbVertices=0; 77 /*Initialize fields*/ 78 this->MaxDepth = 30; 79 this->MaxISize = 1073741824; //2^30 80 this->MaxICoord = 1073741823; //2^30 - 1 81 this->NbQuadtreeBox = 0; 82 this->NbVertices = 0; 119 83 120 84 /*Create container*/ 121 boxcontainer=new DataSet();85 this->boxcontainer=new DataSet(); 122 86 123 87 /*Create Root, pointer toward the main box*/ 124 root=NewBamgQuadtreeBox();88 this->root=NewBamgQuadtreeBox(); 125 89 126 90 /*Check Sizes*/ 127 _assert_( MaxISize>MaxICoor);91 _assert_(this->MaxISize>this->MaxICoord); 128 92 129 93 /*Add all vertices of the mesh*/ 130 if (nbv==-1) nbv=t->nbv; 131 for (int i=0;i<nbv;i++) Add(t->vertices[i]); 94 if(nbv==-1) nbv=t->nbv; 95 for(int i=0;i<nbv;i++){ 96 this->Add(t->vertices[i]); 97 } 132 98 133 99 } … … 140 106 141 107 /*Methods*/ 142 void BamgQuadtree::Add(BamgVertex &w){/*{{{*/108 void BamgQuadtree::Add(BamgVertex &w){/*{{{*/ 143 109 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, BamgQuadtree.cpp/Add)*/ 144 110 BamgQuadtreeBox** pb=NULL; … … 164 130 165 131 //Get next subbox according to the bit value (level) 166 pb = &b->b [IJ(i,j,level)];132 pb = &b->box[BoxNumber(i,j,level)]; 167 133 } 168 134 … … 195 161 196 162 /*Initialize the 4 pointers toward the 4 subboxes*/ 197 b->b [0]=b->b[1]=b->b[2]=b->b[3]=NULL;163 b->box[0]=b->box[1]=b->box[2]=b->box[3]=NULL; 198 164 199 165 /*level = 0010000 -> 0001000*/ … … 205 171 int ij; 206 172 /*bb is the new "sub"box of b where v4[k] is located*/ 207 BamgQuadtreeBox *bb = b->b [ij=IJ(v4[k]->i.x,v4[k]->i.y,level)];173 BamgQuadtreeBox *bb = b->box[ij=BoxNumber(v4[k]->i.x,v4[k]->i.y,level)]; 208 174 209 175 // alloc the BamgQuadtreeBox 210 if (!bb) bb=b->b [ij]=NewBamgQuadtreeBox();176 if (!bb) bb=b->box[ij]=NewBamgQuadtreeBox(); 211 177 212 178 /*Copy the current vertex*/ … … 215 181 216 182 /*Get the subbox where w (i,j) is located*/ 217 pb = &b->b [IJ(i,j,level)];183 pb = &b->box[BoxNumber(i,j,level)]; 218 184 } 219 185 … … 228 194 } 229 195 /*}}}*/ 230 BamgVertex* BamgQuadtree::NearestVertex(Icoor1 i,Icoor1 j) {/*{{{*/ 231 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, BamgQuadtree.cpp/NearestVertex)*/ 232 233 /*Intermediaries*/ 234 BamgQuadtreeBox *pb[MaxDepth]; 235 int pi[MaxDepth]; 236 Icoor1 ii[MaxDepth]; 237 Icoor1 jj[MaxDepth]; 238 int level; 239 long n0; 240 BamgQuadtreeBox *b; 241 long h0; 242 long h = MaxISize; 243 long hb= MaxISize; 244 Icoor1 i0=0,j0=0; 196 int BamgQuadtree::BoxNumber(int i,int j,int l){/*{{{*/ 197 /* 198 * 199 * J j 200 * ^ ^ 201 * | | +--------+--------+ 202 * | | | | | 203 * 1X | | | 2 | 3 | 204 * | | | | | 205 * | | +--------+--------+ 206 * | | | | | 207 * 0X | | | 0 | 1 | 208 * | | | | | 209 * | | +--------+--------+ 210 * | +-----------------------> i 211 * | 212 * |----------------------------> I 213 * X0 X1 214 * 215 * box 0 -> I=0 J=0 BoxNumber=00 = 0 216 * box 1 -> I=1 J=0 BoxNumber=01 = 1 217 * box 2 -> I=0 J=1 BoxNumber=10 = 2 218 * box 3 -> I=1 J=1 BoxNumber=11 = 3 219 */ 220 //BoxNumber(i,j,l) returns the box number of i and j with respect to l 221 //if !j&l and !i&l -> 0 (box zero: lower left ) 222 //if !j&l and i&l -> 1 (box one: lower right) 223 //if j&l and !i&l -> 2 (box two: upper left ) 224 //if j&l and i&l -> 3 (box three:upper right) 225 return ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 )); 226 }/*}}}*/ 227 bool BamgQuadtree::Intersection(int a,int b,int x,int y){/*{{{*/ 228 /*is [x y] in [a b]*/ 229 return ((y) > (a)) && ((x) <(b)); 230 }/*}}}*/ 231 BamgVertex* BamgQuadtree::NearestVertex(int xi,int yi) {/*{{{*/ 245 232 246 233 /*initial output as NULL (no vertex found)*/ 247 234 BamgVertex* nearest_v=NULL; 248 235 249 /*Project w coordinates (i,j) onto [0,MaxISize-1] x [0,MaxISize-1] -> (iplus,jplus)*/250 Icoor1 iplus( i<MaxISize ? (i<0?0:i) : MaxISize-1);251 Icoor1 jplus( j<MaxISize ? (j<0?0:j) : MaxISize-1);252 253 /*Get initial Quadtree box (largest)*/254 b = root;255 256 236 /*if the tree is empty, return NULL pointer*/ 257 if (!root->nbitems) return nearest_v; 237 if(!this->root->nbitems) return nearest_v; 238 239 /*Project coordinates (xi,yi) onto [0,MaxICoord-1] x [0,MaxICoord-1]*/ 240 int xi2 = xi; 241 int yi2 = yi; 242 if(xi<0) xi2 = 0; 243 if(xi>MaxISize) xi2 = MaxICoord; 244 if(yi<0) yi2 = 0; 245 if(yi>MaxISize) yi2 = MaxICoord; 246 247 /*Get inital box (the largest)*/ 248 BamgQuadtreeBox* b = this->root; 249 250 /*Initialize level and sizes for largest box*/ 251 int levelbin = (1L<<this->MaxDepth);// = 2^30 252 int h = this->MaxISize; 253 int hb = this->MaxISize; 254 int i0 = 0; 255 int j0 = 0; 258 256 259 257 /*else, find the smallest non-empty BamgQuadtreeBox containing the point (i,j)*/ 260 while( (n0=b->nbitems)<0){261 262 Icoor1 hb2 = hb >> 1;//size of the current box263 int k = IJ(iplus,jplus,hb2);//box number (0,1,2 or 3)264 BamgQuadtreeBox *b0 = b->b[k];//pointer toward current box265 266 /* break if NULL box orempty (Keep previous box b)*/267 if (( b0 == NULL) || (b0->nbitems ==0)) break;258 while(b->nbitems<0){ 259 260 int hb2 = hb >> 1; //size of the current box 261 int k = BoxNumber(xi2,yi2,hb2); //box number (0,1,2 or 3) 262 BamgQuadtreeBox* b0 = b->box[k]; //pointer toward current box 263 264 /* break if box is empty (Keep previous box b)*/ 265 if((b0==NULL) || (b0->nbitems==0)) break; 268 266 269 267 /*Get next Quadtree box*/ 270 b=b0; 271 i0 += I_IJ(k,hb2); // i orign of BamgQuadtreeBox (macro) 272 j0 += J_IJ(k,hb2); // j orign of BamgQuadtreeBox 273 hb = hb2; // size of the box (in Int) 268 b = b0; 269 hb = hb2; 270 this->SubBoxCoords(&i0,&j0,k,hb2); 274 271 } 275 272 … … 280 277 * - i0: x coordinate of the lower left corner 281 278 * - j0: y coordinate of the lower left corner*/ 279 int n0 = b->nbitems; 282 280 283 281 /* if the current subbox is holding vertices, we are almost done*/ 284 if 285 / /loop over the vertices of the box and find the closest vertex282 if(n0>0){ 283 /*loop over the vertices of the box and find the closest vertex*/ 286 284 for(int k=0;k<n0;k++){ 287 288 /*get integer coordinates of current vertex*/ 289 I2 i2=b->v[k]->i; 290 291 /*Compute norm with w*/ 292 h0=NORM(iplus,i2.x,jplus,i2.y); 285 int xiv = b->v[k]->i.x; 286 int yiv = b->v[k]->i.y; 287 288 289 int h0 = Norm(xi2,xiv,yi2,yiv); 293 290 294 291 /*is it smaller than previous value*/ 295 if 292 if(h0<h){ 296 293 h = h0; 297 294 nearest_v = b->v[k]; … … 306 303 307 304 /*Initialize search variables*/ 305 BamgQuadtreeBox **pb = xNew<BamgQuadtreeBox*>(this->MaxDepth); 306 int* pi = xNew<int>(this->MaxDepth); 307 int* ii = xNew<int>(this->MaxDepth); 308 int* jj = xNew<int>(this->MaxDepth); 308 309 pb[0]=b; //pointer toward the box b 309 310 pi[0]=b->nbitems>0?(int)b->nbitems:4;//number of boxes in b … … 315 316 316 317 /*Main loop*/ 317 level=0; 318 do { 319 318 int level=0; 319 do{ 320 320 /*get current box*/ 321 b = pb[level];321 b = pb[level]; 322 322 323 323 /*Loop over the items in current box (if not empty!)*/ 324 while 324 while(pi[level]){ 325 325 326 326 /*We are looping now over the items of b. k is the current index (in [0 3])*/ … … 330 330 /*if the current subbox is holding vertices (b->nbitems<0 is subboxes)*/ 331 331 if (b->nbitems>0){ 332 I2 i2=b->v[k]->i; 333 h0 = NORM(iplus,i2.x,jplus,i2.y); 332 int h0 = Norm(xi2,b->v[k]->i.x,yi2,b->v[k]->i.y); 334 333 if (h0<h){ 335 334 h=h0; … … 345 344 346 345 /*if the next box exists:*/ 347 if((b=b->b [k])){346 if((b=b->box[k])){ 348 347 349 348 /*Get size (hb) and coordinates of the current sub-box lowest left corner*/ 350 349 hb>>=1; 351 Icoor1 iii = ii[level]+I_IJ(k,hb); 352 Icoor1 jjj = jj[level]+J_IJ(k,hb); 353 354 /*if the current point (iplus,jplus) is in b (modulo h), this box is good: 350 int iii = ii[level]; 351 int jjj = jj[level]; 352 this->SubBoxCoords(&iii,&jjj,k,hb); 353 354 /*if the current point (xi2,yi2) is in b (modulo h), this box is good: 355 355 * it is holding vertices that are close to w */ 356 if (INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)){356 if(this->Intersection(iii,iii+hb,xi2-h,xi2+h) && this->Intersection(jjj,jjj+hb,yi2-h,yi2+h)){ 357 357 level++; 358 pb[level] = b;359 pi[level] = b->nbitems>0 ?(int) b->nbitems :4 ;360 ii[level] = iii;361 jj[level] = jjj;358 pb[level] = b; 359 pi[level] = b->nbitems>0 ? b->nbitems:4 ; 360 ii[level] = iii; 361 jj[level] = jjj; 362 362 } 363 363 364 / /else go backwards364 /*else go backwards*/ 365 365 else{ 366 / /shifted righ by one bit: hb=001000000 -> 01000000366 /*shifted righ by one bit: hb=001000000 -> 01000000*/ 367 367 b=b0; 368 368 hb<<=1; … … 379 379 * in case there is a vertex closest to w that has not yet been tested*/ 380 380 hb <<= 1; 381 } while (level--); 382 383 /*return nearest_v, nearest vertex*/ 381 }while(level--); 382 383 /*return nearest vertex pointer*/ 384 xDelete<BamgQuadtreeBox*>(pb); 385 xDelete<int>(pi); 386 xDelete<int>(ii); 387 xDelete<int>(jj); 384 388 return nearest_v; 385 386 389 } 387 390 /*}}}*/ 388 BamgQuadtree::BamgQuadtreeBox* BamgQuadtree::NewBamgQuadtreeBox(void){/*{{{*/ 389 390 /*Output*/ 391 BamgQuadtreeBox* newbox=NULL; 392 393 /*Create and initialize a new box*/ 394 newbox=new BamgQuadtreeBox; 395 newbox->nbitems=0; 396 newbox->b[0]=NULL; 397 newbox->b[1]=NULL; 398 newbox->b[2]=NULL; 399 newbox->b[3]=NULL; 400 401 /*Add root to the container*/ 402 boxcontainer->AddObject(newbox); 403 404 /*Increase counter*/ 405 NbBamgQuadtreeBox++; 406 407 /*currentbox now points toward next quadtree box*/ 408 return newbox; 391 int BamgQuadtree::Norm(int xi1,int xi2,int yi1,int yi2){/*{{{*/ 392 393 int deltax = xi2 - xi1; 394 int deltay = yi2 - yi1; 395 396 if(deltax<0) deltax = -deltax; 397 if(deltay<0) deltay = -deltay; 398 399 if(deltax> deltay){ 400 return deltax; 401 } 402 else{ 403 return deltay; 404 } 409 405 }/*}}}*/ 410 BamgVertex* BamgQuadtree::ToClose(BamgVertex & v,double seuil,Icoor1 hx,Icoor1 hy){/*{{{*/ 411 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, BamgQuadtree.cpp/ToClose)*/ 412 413 const Icoor1 i=v.i.x; 414 const Icoor1 j=v.i.y; 415 const R2 X(v.r); 416 const Metric Mx(v.m); 417 418 BamgQuadtreeBox * pb[ MaxDepth ]; 419 int pi[ MaxDepth ]; 420 Icoor1 ii[ MaxDepth ], jj [ MaxDepth]; 406 void BamgQuadtree::SubBoxCoords(int* pi,int*pj,int boxnumber,int length){/*{{{*/ 407 /* 408 * j (first bit) 409 * ^ 410 * | +--------+--------+ 411 * | | | | 412 * 1X | | 2 | 3 | 413 * | | | | 414 * | +--------+--------+ 415 * | | | | 416 * 0X | | 0 | 1 | 417 * | | | | 418 * | +--------+--------+ 419 * +-----------------------> i (second bit) 420 * X0 X1 421 */ 422 423 /*Add sub-box coordinate to i and j*/ 424 //_assert_(!(*pi & length)); 425 //_assert_(!(*pj & length)); 426 427 /*length if first bit of boxnumber is 1, elengthse 0*/ 428 *pi += ((boxnumber & 1) ? length:0); 429 /*length if second bit of boxnumber is 1, elengthse 0*/ 430 *pj += ((boxnumber & 2) ? length:0); 431 432 }/*}}}*/ 433 BamgVertex* BamgQuadtree::TooClose(BamgVertex* v,double threshold,int hx,int hy){/*{{{*/ 434 435 const int i=v->i.x; 436 const int j=v->i.y; 437 const double Xx = v->r.x; 438 const double Xy = v->r.y; 439 Metric* Mx = new Metric(v->m); 440 441 BamgQuadtreeBox *pb[MaxDepth]; 442 int pi[MaxDepth]; 443 int ii[MaxDepth], jj[MaxDepth]; 421 444 int l=0; // level 422 BamgQuadtreeBox * b; 423 Icoor1 hb = MaxISize; 424 Icoor1 i0=0,j0=0; 425 426 // BamgVertex *vn=0; 427 428 if (!root->nbitems) 429 return 0; // empty tree 430 431 // general case ----- 445 BamgQuadtreeBox* b; 446 int hb = MaxISize; 447 int i0=0,j0=0; 448 449 // BamgVertex *vn=0; 450 451 if(!root->nbitems) return 0; // empty tree 452 453 // general case 432 454 pb[0]=root; 433 pi[0]=root->nbitems>0 ?(int) root->nbitems : 4;455 pi[0]=root->nbitems>0 ?(int)root->nbitems:4; 434 456 ii[0]=i0; 435 457 jj[0]=j0; 436 do {458 do{ 437 459 b= pb[l]; 438 while 460 while(pi[l]--){ 439 461 int k = pi[l]; 440 462 441 if (b->nbitems>0){ // BamgVertex BamgQuadtreeBox none empty 442 I2 i2 = b->v[k]->i; 443 if ( ABS(i-i2.x) <hx && ABS(j-i2.y) <hy ) 444 { 445 R2 XY(X,b->v[k]->r); 446 if(LengthInterpole(Mx(XY), b->v[k]->m(XY)) < seuil){ 463 if(b->nbitems>0){ // BamgVertex BamgQuadtreeBox none empty 464 int i2x = b->v[k]->i.x; 465 int i2y = b->v[k]->i.y; 466 if (ABS(i-i2x)<hx && ABS(j-i2y) <hy ){ 467 double XYx = b->v[k]->r.x - Xx; 468 double XYy = b->v[k]->r.y - Xy; 469 if(LengthInterpole(Mx->Length(XYx,XYy),b->v[k]->m.Length(XYx,XYy)) < threshold){ 470 delete Mx; 447 471 return b->v[k]; 448 472 } 449 473 } 450 474 } 451 475 else{ // Pointer BamgQuadtreeBox 452 476 BamgQuadtreeBox *b0=b; 453 if ((b=b->b [k])){477 if ((b=b->box[k])){ 454 478 hb >>=1 ; // div by 2 455 long iii = ii[l]+I_IJ(k,hb); 456 long jjj = jj[l]+J_IJ(k,hb); 457 458 if (INTER_SEG(iii,iii+hb,i-hx,i+hx) && INTER_SEG(jjj,jjj+hb,j-hy,j+hy)){ 479 int iii = ii[l]; 480 int jjj = jj[l]; 481 this->SubBoxCoords(&iii,&jjj,k,hb); 482 483 if(this->Intersection(iii,iii+hb,i-hx,i+hx) && this->Intersection(jjj,jjj+hb,j-hy,j+hy)){ 459 484 pb[++l]= b; 460 485 pi[l]= b->nbitems>0 ?(int) b->nbitems : 4 ; … … 474 499 } 475 500 hb <<= 1; // mul by 2 476 } while (l--); 477 501 }while(l--); 502 503 delete Mx; 478 504 return 0; 479 505 } 480 506 /*}}}*/ 507 508 BamgQuadtree::BamgQuadtreeBox* BamgQuadtree::NewBamgQuadtreeBox(void){/*{{{*/ 509 510 /*Output*/ 511 BamgQuadtreeBox* newbox=NULL; 512 513 /*Create and initialize a new box*/ 514 newbox=new BamgQuadtreeBox; 515 newbox->nbitems=0; 516 newbox->box[0]=NULL; 517 newbox->box[1]=NULL; 518 newbox->box[2]=NULL; 519 newbox->box[3]=NULL; 520 521 /*Add root to the container*/ 522 boxcontainer->AddObject(newbox); 523 this->NbQuadtreeBox++; 524 525 /*currentbox now points toward next quadtree box*/ 526 return newbox; 527 }/*}}}*/ 481 528 } -
issm/trunk-jpl/src/c/bamg/BamgQuadtree.h
r19198 r21623 5 5 #include "./include.h" 6 6 #include "../datastructures/datastructures.h" 7 class BamgVertex; 7 8 8 9 namespace bamg { 9 10 const int MaxDepth = 30;11 const long MaxISize = ( 1L << MaxDepth); // = 2^30 : 010000000000..000 (bitwise operation)12 13 class BamgVertex;14 15 10 class BamgQuadtree{ 16 11 … … 24 19 class BamgQuadtreeBox: public Object{ 25 20 public: 26 int nbitems; // number of current vertices in the box 27 union{ 28 BamgQuadtreeBox* b[4]; 29 BamgVertex* v[4]; 30 }; 21 int nbitems; 22 BamgQuadtreeBox *box[4]; 23 BamgVertex *v[4]; 31 24 /*Object functions*/ 32 25 void Echo() {_error_("not implemented yet"); }; … … 44 37 45 38 /*BamgQuadtree public Fields*/ 46 BamgQuadtreeBox* root; 47 long NbBamgQuadtreeBox; 48 long NbVertices; 39 int MaxDepth; // maximum number of subdivision 40 int MaxICoord; // maximum integer coordinate 2^MaxDepth -1 41 int MaxISize; // maximum integer coordinate 2^MaxDepth 42 BamgQuadtreeBox *root; // main box 43 long NbQuadtreeBox; // total number of boxes 44 long NbVertices; // number of points 49 45 50 46 BamgQuadtree(); … … 52 48 ~BamgQuadtree(); 53 49 54 BamgVertex *NearestVertex(Icoor1 i,Icoor1 j); 50 void Add(BamgVertex &w); 51 int BoxNumber(int i,int j,int l); 52 bool Intersection(int a,int b,int x,int y); 53 BamgVertex *NearestVertex(int i,int j); 55 54 BamgQuadtreeBox *NewBamgQuadtreeBox(void); 56 BamgVertex *ToClose(BamgVertex &,double ,Icoor1,Icoor1); 57 void Add(BamgVertex &w); 55 int Norm(int xi1,int yi1,int xi2,int yi2); 56 void SubBoxCoords(int* pi,int*pj,int boxcoord,int length); 57 BamgVertex *TooClose(BamgVertex*,double,int,int); 58 58 }; 59 59 } -
issm/trunk-jpl/src/c/bamg/BamgVertex.cpp
r18455 r21623 115 115 ret = t->Optim(IndexInTriangle,koption); 116 116 if(!i){ 117 t =0; // for no future optim e117 t =0; // for no future optim 118 118 IndexInTriangle= 0; 119 119 } … … 196 196 BamgVertex *v0,*v1,*v2,*v3; 197 197 if (tria->det<0) ok =1; 198 else if (tria->Quadrangle(v0,v1,v2,v3))199 {200 vP = vPsave;201 double qold =QuadQuality(*v0,*v1,*v2,*v3);202 vP = vPnew;203 double qnew =QuadQuality(*v0,*v1,*v2,*v3);204 if (qnew<qold) ok = 1;205 }206 198 else if ( (double)tria->det < detold/2 ) ok=1; 207 208 199 } 209 200 tria->SetUnMarkUnSwap(0); … … 225 216 } 226 217 /*}}}*/ 227 228 /*Intermediary*/229 double QuadQuality(const BamgVertex & a,const BamgVertex &b,const BamgVertex &c,const BamgVertex &d) {/*{{{*/230 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/QuadQuality)*/231 232 // calcul de 4 angles --233 R2 A((R2)a),B((R2)b),C((R2)c),D((R2)d);234 R2 AB(B-A),BC(C-B),CD(D-C),DA(A-D);235 // Move(A),Line(B),Line(C),Line(D),Line(A);236 const Metric & Ma = a;237 const Metric & Mb = b;238 const Metric & Mc = c;239 const Metric & Md = d;240 241 double lAB=Norme2(AB);242 double lBC=Norme2(BC);243 double lCD=Norme2(CD);244 double lDA=Norme2(DA);245 AB /= lAB; BC /= lBC; CD /= lCD; DA /= lDA;246 // version aniso247 double cosDAB= Ma(DA,AB)/(Ma(DA)*Ma(AB)),sinDAB= Det(DA,AB);248 double cosABC= Mb(AB,BC)/(Mb(AB)*Mb(BC)),sinABC= Det(AB,BC);249 double cosBCD= Mc(BC,CD)/(Mc(BC)*Mc(CD)),sinBCD= Det(BC,CD);250 double cosCDA= Md(CD,DA)/(Md(CD)*Md(DA)),sinCDA= Det(CD,DA);251 double sinmin=Min(Min(sinDAB,sinABC),Min(sinBCD,sinCDA));252 if (sinmin<=0) return sinmin;253 return 1.0-Max(Max(Abs(cosDAB),Abs(cosABC)),Max(Abs(cosBCD),Abs(cosCDA)));254 }255 /*}}}*/256 257 218 } -
issm/trunk-jpl/src/c/bamg/BamgVertex.h
r19293 r21623 4 4 #include "./include.h" 5 5 #include "./Metric.h" 6 #include "./Direction.h"7 6 #include "./BamgOpts.h" 8 7 … … 24 23 long ReferenceNumber; 25 24 long PreviousNumber; 26 Direction DirOfSearch;27 25 short IndexInTriangle; // the vertex number in triangle; varies between 0 and 2 in t 28 26 … … 40 38 operator const R2 & () const {return r;} // Cast operator 41 39 operator Metric () const {return m;} // Cast operator 42 double operator()(R2 x) const { return m (x);} // Get x in the metric m40 double operator()(R2 x) const { return m.Length(x.x,x.y);} // Get x in the metric m 43 41 44 42 /*methods (No constructor and no destructors...)*/ … … 54 52 inline void Set(const BamgVertex &rec,const Mesh & ,Mesh & ){*this=rec;} 55 53 }; 56 57 //Intermediary58 double QuadQuality(const BamgVertex &,const BamgVertex &,const BamgVertex &,const BamgVertex &);59 54 } 60 55 #endif -
issm/trunk-jpl/src/c/bamg/EigenMetric.cpp
r18064 r21623 37 37 lambda1=0; 38 38 lambda2=0; 39 v .x=1;40 v .y=0;39 vx=1; 40 vy=0; 41 41 } 42 42 /*2: delta is small -> double root*/ … … 44 44 lambda1=-b/2; 45 45 lambda2=-b/2; 46 v .x=1;47 v .y=0;46 vx=1; 47 vy=0; 48 48 } 49 49 /*3: general case -> two roots*/ … … 76 76 if (norm2<norm1){ 77 77 norm1=sqrt(norm1); 78 v .x = - a21/norm1;79 v .y = (a11-lambda1)/norm1;78 vx = - a21/norm1; 79 vy = (a11-lambda1)/norm1; 80 80 } 81 81 else{ 82 82 norm2=sqrt(norm2); 83 v .x = - (a22-lambda1)/norm2;84 v .y = a21/norm2;83 vx = - (a22-lambda1)/norm2; 84 vy = a21/norm2; 85 85 } 86 86 } … … 88 88 } 89 89 /*}}}*/ 90 EigenMetric::EigenMetric(double r1,double r2,const D2& vp1): lambda1(r1),lambda2(r2),v(vp1){/*{{{*/ 91 90 EigenMetric::EigenMetric(double r1,double r2,const D2& vp1){/*{{{*/ 91 this->lambda1 = r1; 92 this->lambda2 = r2; 93 this->vx = vp1.x; 94 this->vy = vp1.y; 92 95 }/*}}}*/ 93 96 … … 104 107 _printf_(" lambda1: " << lambda1 << "\n"); 105 108 _printf_(" lambda2: " << lambda2 << "\n"); 106 _printf_(" v .x: " << v.x << "\n");107 _printf_(" v .y: " << v.y << "\n");109 _printf_(" vx: " << vx << "\n"); 110 _printf_(" vy: " << vy << "\n"); 108 111 109 112 return; -
issm/trunk-jpl/src/c/bamg/Geometry.cpp
r18064 r21623 8 8 9 9 namespace bamg { 10 11 static const Direction NoDirOfSearch=Direction();12 10 13 11 /*Constructors/Destructors*/ … … 80 78 vertices[i].r.y=(double)bamggeom->Vertices[i*3+1]; 81 79 vertices[i].ReferenceNumber=(long)bamggeom->Vertices[i*3+2]; 82 vertices[i].DirOfSearch=NoDirOfSearch;83 80 vertices[i].color =0; 84 81 vertices[i].type=0; … … 772 769 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/ 773 770 /*Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge*/ 771 /*FIXME: should go away*/ 774 772 775 773 double save_s=s; … … 807 805 808 806 //Compute metric of new vertex as a linear interpolation of the two others 809 V.m=Metric(1.0-s,v0 ,s,v1);807 V.m=Metric(1.0-s,v0.m,s,v1.m); 810 808 811 809 const int mxe=100; -
issm/trunk-jpl/src/c/bamg/Mesh.cpp
r21299 r21623 9 9 10 10 namespace bamg { 11 12 static const Direction NoDirOfSearch=Direction();13 11 14 12 /*Constructors/Destructors*/ … … 167 165 nbsubdomains = Th.nbsubdomains; 168 166 nbtout = Th.nbtout; 169 nbq = Th.nbq ;170 167 NbVerticesOnGeomVertex = Th.NbVerticesOnGeomVertex; 171 168 if(NbVerticesOnGeomVertex) … … 269 266 vertices[i].r.y=y[i]; 270 267 vertices[i].ReferenceNumber=1; 271 vertices[i].DirOfSearch =NoDirOfSearch;272 268 vertices[i].m=M1; 273 269 vertices[i].color=0; … … 336 332 vertices[i].r.y=bamgmesh->Vertices[i*3+1]; 337 333 vertices[i].ReferenceNumber=(long)bamgmesh->Vertices[i*3+2]; 338 vertices[i].DirOfSearch =NoDirOfSearch;339 334 vertices[i].m=M1; 340 335 vertices[i].color=0; … … 362 357 else{ 363 358 if(verbose>5) _error_("no Triangles found in the initial mesh"); 364 }365 366 //Quadrilaterals367 if(bamgmesh->Quadrilaterals){368 if(verbose>5) _printf_(" processing Quadrilaterals\n");369 long i1,i2,i3,i4;370 triangles =new Triangle[nbt];371 for (i=0;i<bamgmesh->QuadrilateralsSize[0];i++){372 //divide the quad into two triangles373 Triangle & t1 = triangles[2*i];374 Triangle & t2 = triangles[2*i+1];375 i1=(long)bamgmesh->Quadrilaterals[i*5+0]-1; //for C indexing376 i2=(long)bamgmesh->Quadrilaterals[i*5+1]-1; //for C indexing377 i3=(long)bamgmesh->Quadrilaterals[i*5+2]-1; //for C indexing378 i4=(long)bamgmesh->Quadrilaterals[i*5+3]-1; //for C indexing379 t1=Triangle(this,i1,i2,i3);380 t2=Triangle(this,i3,i4,i1);381 t1.color=(long)bamgmesh->Quadrilaterals[i*5+4];382 t2.color=(long)bamgmesh->Quadrilaterals[i*5+4];383 t1.SetHidden(OppositeEdge[1]); // two times because the adj was not created384 t2.SetHidden(OppositeEdge[1]);385 }386 359 } 387 360 … … 698 671 /*Triangles*/ 699 672 if(verbose>5) _printf_(" writing Triangles\n"); 700 k=nbInT -nbq*2;673 k=nbInT; 701 674 num=0; 702 675 bamgmesh->TrianglesSize[0]=k; … … 713 686 bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ReferenceNumber; 714 687 num=num+1; 715 }716 }717 }718 719 /*Quadrilaterals*/720 if(verbose>5) _printf_(" writing Quadrilaterals\n");721 bamgmesh->QuadrilateralsSize[0]=nbq;722 bamgmesh->QuadrilateralsSize[1]=5;723 if (nbq){724 bamgmesh->Quadrilaterals=xNew<double>(5*nbq);725 for (i=0;i<nbt;i++){726 Triangle &t =triangles[i];727 Triangle* ta;728 BamgVertex *v0,*v1,*v2,*v3;729 if (reft[i]<0) continue;730 if ((ta=t.Quadrangle(v0,v1,v2,v3)) !=0 && &t<ta) {731 bamgmesh->Quadrilaterals[i*5+0]=GetId(v0)+1; //back to M indexing732 bamgmesh->Quadrilaterals[i*5+1]=GetId(v1)+1; //back to M indexing733 bamgmesh->Quadrilaterals[i*5+2]=GetId(v2)+1; //back to M indexing734 bamgmesh->Quadrilaterals[i*5+3]=GetId(v3)+1; //back to M indexing735 bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]].ReferenceNumber;736 688 } 737 689 } … … 1113 1065 // s0 tt2 s1 1114 1066 //------------------------------- 1115 1067 1116 1068 /*Intermediaries*/ 1117 1069 Triangle* tt[3]; //the three triangles … … 1230 1182 } 1231 1183 } 1232 } 1233 /*}}}*/1184 1185 }/*}}}*/ 1234 1186 void Mesh::BoundAnisotropy(double anisomax,double hminaniso) {/*{{{*/ 1235 1187 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/BoundAnisotropy)*/ … … 2371 2323 _printf_(" nbt = " << nbt << "\n"); 2372 2324 _printf_(" nbe = " << nbe << "\n"); 2373 _printf_(" nbq = " << nbq << "\n");2374 2325 _printf_(" index:\n"); 2375 2326 for (i=0;i<nbt;i++){ 2376 2327 _printf_(" " << setw(4) << i+1 << ": [" 2377 << setw(4) << (((BamgVertex 2378 << setw(4) << (((BamgVertex *)triangles[i](0))?GetId(triangles[i][1])+1:0) << " "2379 << setw(4) << (((BamgVertex *)triangles[i](0))?GetId(triangles[i][2])+1:0) << "]");2328 << setw(4) << (((BamgVertex*)triangles[i](0))?GetId(triangles[i][0])+1:0) << " " 2329 << setw(4) << (((BamgVertex*)triangles[i](1))?GetId(triangles[i][1])+1:0) << " " 2330 << setw(4) << (((BamgVertex*)triangles[i](2))?GetId(triangles[i][2])+1:0) << "]\n"); 2380 2331 } 2381 2332 _printf_(" coordinates:\n"); 2382 2333 for (i=0;i<nbv;i++){ 2383 _printf_(" " << setw(4) << i+1 << ": [" << vertices[i].r.x << " " << vertices[i].r.y << "] \n");2334 _printf_(" " << setw(4) << i+1 << ": [" << vertices[i].r.x << " " << vertices[i].r.y << "] and [" << vertices[i].i.x << " " << vertices[i].i.y << "]\n"); 2384 2335 } 2385 2336 … … 2707 2658 nbe=0; 2708 2659 edges=NULL; 2709 nbq=0;2710 2660 nbsubdomains=0; 2711 2661 subdomains=NULL; … … 2770 2720 * let's show that: 2771 2721 * 2772 * for all j in [0 nbv[, ∃!i in [0 nbv[ such that k_i=j2722 * for all j in [0 nbv[, there is a unique i in [0 nbv[ such that k_i=j 2773 2723 * 2774 2724 * Let's assume that there are 2 distinct j1 and j2 such that … … 2914 2864 vi.m.Box(hx,hy); 2915 2865 Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor); 2916 if(!quadtree->To Close(vi,seuil,hi,hj)){2866 if(!quadtree->TooClose(&vi,seuil,hi,hj)){ 2917 2867 // a good new point 2918 2868 BamgVertex &vj = vertices[iv]; … … 2999 2949 } 3000 2950 /*}}}*/ 3001 void Mesh::MakeQuadrangles(double costheta){/*{{{*/3002 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadrangles)*/3003 3004 long int verbose=0;3005 3006 if (verbose>2) _printf_("MakeQuadrangles costheta = " << costheta << "\n");3007 3008 if (costheta >1) {3009 if (verbose>5) _printf_(" do nothing: costheta > 1\n");3010 }3011 3012 long nbqq = (nbt*3)/2;3013 DoubleAndInt *qq = new DoubleAndInt[nbqq];3014 3015 long i,ij;3016 int j;3017 long k=0;3018 for (i=0;i<nbt;i++)3019 for (j=0;j<3;j++)3020 if ((qq[k].q=triangles[i].QualityQuad(j))>=costheta)3021 qq[k++].i3j=i*3+j;3022 // sort qq3023 HeapSort(qq,k);3024 3025 long kk=0;3026 for (ij=0;ij<k;ij++) {3027 i=qq[ij].i3j/3;3028 j=(int) (qq[ij].i3j%3);3029 // optisamition no float computation3030 if (triangles[i].QualityQuad(j,0) >=costheta)3031 triangles[i].SetHidden(j),kk++;3032 }3033 nbq = kk;3034 if (verbose>2){3035 _printf_(" number of quadrilaterals = " << nbq << "\n");3036 _printf_(" number of triangles = " << nbt-nbtout- nbq*2 << "\n");3037 _printf_(" number of outside triangles = " << nbtout << "\n");3038 }3039 delete [] qq;3040 }3041 /*}}}*/3042 2951 void Mesh::MakeBamgQuadtree() { /*{{{*/ 3043 2952 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeBamgQuadtree)*/ … … 3073 2982 lmax = Max(lmax,l); 3074 2983 if(l> maxsubdiv2){ 3075 2984 R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M 3076 2985 double lc = M(AC,AC); 3077 2986 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; … … 3086 2995 lmax = Max(lmax,l); 3087 2996 if(l> maxsubdiv2){ 3088 2997 R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M 3089 2998 double lc = M(AC,AC); 3090 2999 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; … … 3857 3766 if (tstart[i] != &vide) // not a boundary vertex 3858 3767 delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega)); 3859 if (!nbq) 3860 for ( i=0;i<nbv;i++) 3861 if (tstart[i] != &vide) // not a boundary vertex 3862 NbSwap += vertices[i].Optim(1); 3768 for ( i=0;i<nbv;i++) 3769 if (tstart[i] != &vide) // not a boundary vertex 3770 NbSwap += vertices[i].Optim(1); 3863 3771 if (verbose>3) _printf_(" move max = " << pow(delta,0.5) << ", iteration = " << k << ", nb of swap = " << NbSwap << "\n"); 3864 3772 } … … 3914 3822 double ll = Norme2(Aij); 3915 3823 if (0) { 3916 double hi = ll/vi.m (Aij);3917 double hj = ll/vj.m (Aij);3824 double hi = ll/vi.m.Length(Aij.x,Aij.y); 3825 double hj = ll/vj.m.Length(Aij.x,Aij.y); 3918 3826 if(hi < hj) 3919 3827 { … … 3928 3836 else 3929 3837 { 3930 double li = vi.m (Aij);3838 double li = vi.m.Length(Aij.x,Aij.y); 3931 3839 if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) ) 3932 3840 if(first_np_or_next_t1[j]<0) // if the metrix change … … 3973 3881 vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m); 3974 3882 vertices[nbv].ReferenceNumber=0; 3975 vertices[nbv].DirOfSearch = NoDirOfSearch ;3976 3883 } 3977 3884 NbSplitEdge++; … … 3992 3899 // a good new point 3993 3900 vi.ReferenceNumber=0; 3994 vi.DirOfSearch =NoDirOfSearch;3995 3901 Triangle *tcvi = TriangleFindFromCoord(vi.i,det3); 3996 3902 if (tcvi && !tcvi->link) { … … 4181 4087 vertices[i].r.y=y[i]; 4182 4088 vertices[i].ReferenceNumber=1; 4183 vertices[i].DirOfSearch =NoDirOfSearch;4184 4089 vertices[i].m=M1; 4185 4090 vertices[i].color=0; … … 4331 4236 Metric MA = background ? BTh.MetricAt(a->r) :a->m ; //Get metric associated to A 4332 4237 Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B 4333 double ledge = (MA (AB) + MB(AB))/2; //Get edge length in metric4238 double ledge = (MA.Length(AB.x,AB.y) + MB.Length(AB.x,AB.y))/2; //Get edge length in metric 4334 4239 4335 4240 /* We are now creating the mesh edges from the geometrical edge selected above. … … 4367 4272 MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB); 4368 4273 AB = A-B; 4369 lSubEdge[kk]=(ledge+=(MAs (AB)+MBs(AB))/2);4274 lSubEdge[kk]=(ledge+=(MAs.Length(AB.x,AB.y)+MBs.Length(AB.x,AB.y))/2); 4370 4275 } 4371 4276 } … … 4412 4317 vb->m = Metric(aa,a->m,bb,b->m); 4413 4318 vb->ReferenceNumber = e->ReferenceNumber; 4414 vb->DirOfSearch =NoDirOfSearch;4415 4319 double abcisse = k ? bb : aa; 4416 4320 vb->r = e->F(abcisse); … … 4735 4639 ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); 4736 4640 A1->ReferenceNumber = eeequi.ReferenceNumber; 4737 A1->DirOfSearch =NoDirOfSearch;4738 4641 e->GeomEdgeHook = ongequi; 4739 4642 e->v[0]=A0; -
issm/trunk-jpl/src/c/bamg/Mesh.h
r16967 r21623 33 33 long NbRef; // counter of ref on the this class if 0 we can delete 34 34 long maxnbv,maxnbt; // nombre max de sommets , de triangles 35 long nbv,nbt,nbe ,nbq; // nb of vertices, of triangles, of edges and quadrilaterals35 long nbv,nbt,nbe; // nb of vertices, of triangles and edges 36 36 long nbsubdomains; 37 37 long nbtout; // Nb of oudeside triangle … … 86 86 void SmoothMetric(double raisonmax) ; 87 87 void BoundAnisotropy(double anisomax,double hminaniso= 1e-100) ; 88 void MaxSubDivision(double maxsubdiv);89 88 Edge** MakeGeomEdgeToEdge(); 90 89 long SplitInternalEdgeWithBorderVertices(); 91 void MakeQuadrangles(double costheta);92 90 void MakeBamgQuadtree(); 91 void MaxSubDivision(double maxsubdiv); 93 92 void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1); 94 93 long InsertNewPoints(long nbvold,long & NbTSwap,bool random); -
issm/trunk-jpl/src/c/bamg/Metric.cpp
r18064 r21623 13 13 14 14 /*Constructor/Destructor*/ 15 Metric::Metric(double a): a11(1/(a*a)),a21(0),a22(1/(a*a)){/*{{{*/ 15 Metric::Metric(double a){/*{{{*/ 16 17 /*Isotropic metric for a length of a as unit*/ 18 this->a11 = 1./(a*a); 19 this->a21 = 0.; 20 this->a22 = 1./(a*a); 16 21 17 22 }/*}}}*/ 18 Metric::Metric(double a,double b,double c) :a11(a),a21(b),a22(c){/*{{{*/ 23 Metric::Metric(double a11_in,double a21_in,double a22_in){/*{{{*/ 24 25 this->a11 = a11_in; 26 this->a21 = a21_in; 27 this->a22 = a22_in; 19 28 20 29 }/*}}}*/ 21 Metric::Metric(const double a[3],const Metric& m0, const Metric& m1,const Metric&m2 ){/*{{{*/22 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/Metric)*/ 23 30 Metric::Metric(const double a[3],const Metric m0,const Metric m1,const Metric m2 ){/*{{{*/ 31 32 /*Create metric using 3 inputs*/ 24 33 Metric mab(a[0]*m0.a11 + a[1]*m1.a11 + a[2]*m2.a11, 25 34 a[0]*m0.a21 + a[1]*m1.a21 + a[2]*m2.a21, 26 35 a[0]*m0.a22 + a[1]*m1.a22 + a[2]*m2.a22); 27 36 37 /*Convert to eigen metric*/ 28 38 EigenMetric vab(mab); 29 30 R2 v1(vab.v.x,vab.v.y); 31 R2 v2(-v1.y,v1.x); 32 33 double h1 = a[0] / m0(v1) + a[1] / m1(v1) + a[2] / m2(v1); 34 double h2 = a[0] / m0(v2) + a[1] / m1(v2) + a[2] / m2(v2); 35 36 vab.lambda1 = 1 / (h1*h1); 37 vab.lambda2 = 1 / (h2*h2); 38 *this = vab; 39 } 40 /*}}}*/ 41 Metric::Metric(double a,const Metric& ma, double b,const Metric& mb) { /*{{{*/ 42 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/EigenMetric)*/ 39 double v1x = + vab.vx; 40 double v1y = + vab.vy; 41 double v2x = - vab.vy; 42 double v2y = + vab.vx; 43 44 double h1=a[0] / m0.Length(v1x,v1y) + a[1] / m1.Length(v1x,v1y) + a[2] / m2.Length(v1x,v1y); 45 double h2=a[0] / m0.Length(v2x,v2y) + a[1] / m1.Length(v2x,v2y) + a[2] / m2.Length(v2x,v2y); 46 47 vab.lambda1 = 1. / (h1*h1); 48 vab.lambda2 = 1. / (h2*h2); 49 50 /*Build metric from vab*/ 51 double v00=vab.vx*vab.vx; 52 double v11=vab.vy*vab.vy; 53 double v01=vab.vx*vab.vy; 54 this->a11=v00*vab.lambda1+v11*vab.lambda2; 55 this->a21=v01*(vab.lambda1-vab.lambda2); 56 this->a22=v00*vab.lambda2+v11*vab.lambda1; 57 58 } /*}}}*/ 59 Metric::Metric(double a,const Metric ma,double b,const Metric mb) { /*{{{*/ 43 60 44 61 /*Compute metric (linear combination of ma and mb)*/ … … 47 64 /*Get Eigen values and vectors*/ 48 65 EigenMetric vab(mab); 49 R2 v1(vab.v.x,vab.v.y); 50 R2 v2(-v1.y,v1.x); 66 double v1x = + vab.vx; 67 double v1y = + vab.vy; 68 double v2x = - vab.vy; 69 double v2y = + vab.vx; 51 70 52 71 /*Modify eigen values (a+b=1)*/ 53 double h1 = a/ma(v1) + b/mb(v1); 54 double h2 = a/ma(v2) + b/mb(v2); 55 vab.lambda1 = 1/(h1*h1); 56 vab.lambda2 = 1/(h2*h2); 57 *this=vab; 72 double h1 = a/ma.Length(v1x,v1y) + b/mb.Length(v1x,v1y); 73 double h2 = a/ma.Length(v2x,v2y) + b/mb.Length(v2x,v2y); 74 vab.lambda1 = 1./(h1*h1); 75 vab.lambda2 = 1./(h2*h2); 76 77 /*Build metric from vab*/ 78 double v00=vab.vx*vab.vx; 79 double v11=vab.vy*vab.vy; 80 double v01=vab.vx*vab.vy; 81 this->a11=v00*vab.lambda1+v11*vab.lambda2; 82 this->a21=v01*(vab.lambda1-vab.lambda2); 83 this->a22=v00*vab.lambda2+v11*vab.lambda1; 58 84 } 59 85 /*}}}*/ … … 123 149 } 124 150 /*}}}*/ 125 R2 Metric::mul(const R2 x)const {/*{{{*/ 126 return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y); 127 }/*}}}*/ 151 double Metric::Length(double Ax,double Ay) const{/*{{{*/ 152 /*Length of A in the current metric*/ 153 return sqrt(Ax*Ax*a11+2*Ax*Ay*a21+Ay*Ay*a22); 154 } 155 /*}}}*/ 128 156 129 157 /*Intermediary*/ … … 140 168 Ms1[level]=Ma; 141 169 Ms2[level]=Mb; 142 double sa = Ma (AB);143 double sb = Mb (AB);170 double sa = Ma.Length(AB.x,AB.y); 171 double sb = Mb.Length(AB.x,AB.y); 144 172 lMs1[level]=sa; 145 173 lMs2[level]=sb; … … 160 188 if( s > sstop && level < 30 && i < 500-level ) { 161 189 Metric Mi(0.5,M1,0.5,M2); 162 double si = Mi (AB);190 double si = Mi.Length(AB.x,AB.y); 163 191 if( Abs((s1+s2)-(si+si)) > s1*0.001) 164 192 { -
issm/trunk-jpl/src/c/bamg/Metric.h
r16237 r21623 30 30 Metric(double a); 31 31 Metric(double a,double b,double c); 32 Metric( double a,const Metric& ma, double b,const Metric&mb);33 Metric(const double a[3],const Metric& m0,const Metric& m1,const Metric&m2 );32 Metric(double a,const Metric ma,double b,const Metric mb); 33 Metric(const double a[3],const Metric m0,const Metric m1,const Metric m2 ); 34 34 void Echo(); 35 R2 mul(const R2 x)const;36 35 double det() const; 37 36 int IntersectWith(const Metric& M2); … … 42 41 R2 Orthogonal(const R2 x){ return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y); } 43 42 R2 Orthogonal(const I2 x){ return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y); } 43 double Length(double Ax,double Ay) const; 44 44 45 45 //operators … … 47 47 Metric operator/(double c) const {double c2=1/(c*c);return Metric(a11*c2,a21*c2,a22*c2);} 48 48 operator D2xD2(){ return D2xD2(a11,a21,a21,a22);} 49 double operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; // length of x in metric sqrt(<Mx,x>)49 //double operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; // length of x in metric sqrt(<Mx,x>) FIXME: replace by Length! 50 50 double operator()(R2 x,R2 y) const { return x.x*y.x*a11+(x.x*x.y+x.y*y.x)*a21+x.y*y.y*a22;}; 51 51 … … 57 57 //fields 58 58 double lambda1,lambda2; 59 D2 v;59 double vx,vy; 60 60 61 61 //friends … … 114 114 } 115 115 inline Metric::Metric(const EigenMetric& M) { 116 double v00=M.v .x*M.v.x;117 double v11=M.v .y*M.v.y;118 double v01=M.v .x*M.v.y;116 double v00=M.vx*M.vx; 117 double v11=M.vy*M.vy; 118 double v01=M.vx*M.vy; 119 119 a11=v00*M.lambda1+v11*M.lambda2; 120 120 a21=v01*(M.lambda1-M.lambda2); -
issm/trunk-jpl/src/c/bamg/Triangle.cpp
r18064 r21623 145 145 } 146 146 /*}}}*/ 147 Triangle* Triangle::Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const{/*{{{*/148 // return the other triangle of the quad if a quad or 0 if not a quat149 Triangle * t =0;150 if (link) {151 int a=-1;152 if (AdjEdgeIndex[0] & 16 ) a=0;153 if (AdjEdgeIndex[1] & 16 ) a=1;154 if (AdjEdgeIndex[2] & 16 ) a=2;155 if (a>=0) {156 t = adj[a];157 // if (t-this<0) return 0;158 v2 = vertices[VerticesOfTriangularEdge[a][0]];159 v0 = vertices[VerticesOfTriangularEdge[a][1]];160 v1 = vertices[OppositeEdge[a]];161 v3 = t->vertices[OppositeEdge[AdjEdgeIndex[a]&3]];162 }163 }164 return t;165 }166 /*}}}*/167 double Triangle::QualityQuad(int a,int option) const{/*{{{*/168 double q;169 if (!link || AdjEdgeIndex[a] &4)170 q= -1;171 else {172 Triangle * t = adj[a];173 if (t-this<0) q= -1;// because we do 2 times174 else if (!t->link ) q= -1;175 else if (AdjEdgeIndex[0] & 16 || AdjEdgeIndex[1] & 16 || AdjEdgeIndex[2] & 16 || t->AdjEdgeIndex[0] & 16 || t->AdjEdgeIndex[1] & 16 || t->AdjEdgeIndex[2] & 16 )176 q= -1;177 else if(option){178 const BamgVertex & v2 = *vertices[VerticesOfTriangularEdge[a][0]];179 const BamgVertex & v0 = *vertices[VerticesOfTriangularEdge[a][1]];180 const BamgVertex & v1 = *vertices[OppositeEdge[a]];181 const BamgVertex & v3 = * t->vertices[OppositeEdge[AdjEdgeIndex[a]&3]];182 q = QuadQuality(v0,v1,v2,v3); // do the float part183 }184 else q= 1;185 }186 return q;187 }188 /*}}}*/189 147 void Triangle::Renumbering(Triangle *tb,Triangle *te, long *renu){/*{{{*/ 190 148 … … 262 220 t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] &=55; // 23 + 32 263 221 AdjEdgeIndex[a] &=55 ; 222 }/*}}}*/ 223 Triangle* Triangle::TriangleAdj(int i) const {/*{{{*/ 224 return adj[i&3]; 264 225 }/*}}}*/ 265 226 int Triangle::swap(short a,int koption){/*{{{*/ … … 305 266 // critere de Delaunay pure isotrope 306 267 Icoor2 xb1 = sb->i.x - s1->i.x, 307 308 309 310 311 312 313 268 x21 = s2->i.x - s1->i.x, 269 yb1 = sb->i.y - s1->i.y, 270 y21 = s2->i.y - s1->i.y, 271 xba = sb->i.x - sa->i.x, 272 x2a = s2->i.x - sa->i.x, 273 yba = sb->i.y - sa->i.y, 274 y2a = s2->i.y - sa->i.y; 314 275 double 315 276 cosb12 = double(xb1*x21 + yb1*y21), … … 342 303 if (Abs(d) > dd*1.e-3) { 343 304 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); 344 som = M (C - S2)/M(C - S1);305 som = M.Length(C.x-S2.x,C.y-S2.y) / M.Length(C.x-S1.x,C.y-S1.y); 345 306 } else 346 307 {kopt=1;continue;} … … 357 318 if(Abs(d) > dd*1.e-3) { 358 319 R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); 359 som += M(C - S2)/M(C - S1);320 som += M.Length(C.x-S2.x,C.y-S2.y) / M.Length(C.x-S1.x,C.y-S1.y); 360 321 } else 361 322 {kopt=1;continue;} … … 376 337 } 377 338 /*}}}*/ 378 Triangle* Triangle::TriangleAdj(int i) const {/*{{{*/379 return adj[i&3];380 }/*}}}*/381 339 382 340 } -
issm/trunk-jpl/src/c/bamg/Triangle.h
r19293 r21623 46 46 int Hidden(int a)const; 47 47 int GetAllflag(int a); 48 double QualityQuad(int a,int option=1) const;49 48 short NuEdgeTriangleAdj(int i) const; 50 49 AdjacentTriangle Adj(int i) const; 51 50 Triangle *TriangleAdj(int i) const; 52 Triangle *Quadrangle(BamgVertex *& v0,BamgVertex *& v1,BamgVertex *& v2,BamgVertex *& v3) const;53 51 void Renumbering(Triangle *tb,Triangle *te, long *renu); 54 52 void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu); -
issm/trunk-jpl/src/c/bamg/bamgobjects.h
r12832 r21623 11 11 #include "./BamgMesh.h" 12 12 #include "./Metric.h" 13 #include "./DoubleAndInt.h"14 #include "./Direction.h"15 13 #include "./BamgVertex.h" 16 14 #include "./AdjacentTriangle.h" -
issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp
r18500 r21623 21 21 int i; 22 22 int noerr=1; 23 double costheta=2;24 23 double hminaniso=1e-100; 25 24 Mesh* Thr=NULL; … … 168 167 if (Thr != &BTh) delete Thr; 169 168 170 //Make quadrangle if requested171 if(costheta<=1.0) Th.MakeQuadrangles(costheta);172 //if () Th.SplitElement(2);173 174 169 //Split corners if requested 175 170 if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices(); … … 183 178 //display info 184 179 if(verbosity>0) { 185 if (Th.nbt-Th.nbtout-Th.nbq*2){ 186 _printf_(" new number of triangles = " << (Th.nbt-Th.nbtout-Th.nbq*2) << "\n"); 187 } 188 if (Th.nbq ){ 189 _printf_(" new number of quads = " << Th.nbq << "\n"); 180 if (Th.nbt-Th.nbtout){ 181 _printf_(" new number of triangles = " << (Th.nbt-Th.nbtout) << "\n"); 190 182 } 191 183 }
Note:
See TracChangeset
for help on using the changeset viewer.
