Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 21622) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 21623) @@ -20,6 +20,38 @@ #Core sources +#BAMG sources {{{ +issm_sources = +if BAMG +issm_sources += ./bamg/BamgGeom.cpp\ + ./bamg/BamgMesh.cpp\ + ./bamg/BamgOpts.cpp\ + ./bamg/CrackedEdge.cpp\ + ./bamg/Curve.cpp\ + ./bamg/Edge.cpp\ + ./bamg/GeomEdge.cpp\ + ./bamg/GeomSubDomain.cpp\ + ./bamg/GeomVertex.cpp\ + ./bamg/Geometry.cpp\ + ./bamg/ListofIntersectionTriangles.cpp\ + ./bamg/EigenMetric.cpp\ + ./bamg/Metric.cpp\ + ./bamg/BamgQuadtree.cpp\ + ./bamg/SetOfE4.cpp\ + ./bamg/SubDomain.cpp\ + ./bamg/AdjacentTriangle.cpp\ + ./bamg/Triangle.cpp\ + ./bamg/BamgVertex.cpp\ + ./bamg/VertexOnEdge.cpp\ + ./bamg/VertexOnGeom.cpp\ + ./bamg/VertexOnVertex.cpp\ + ./bamg/Mesh.cpp\ + ./shared/Bamg/BigPrimeNumber.cpp\ + ./modules/Bamgx/Bamgx.cpp\ + ./modules/BamgConvertMeshx/BamgConvertMeshx.cpp\ + ./modules/BamgTriangulatex/BamgTriangulatex.cpp +endif +#}}} #Core sources{{{ -issm_sources = ./datastructures/DataSet.cpp\ +issm_sources += ./datastructures/DataSet.cpp\ ./classes/gauss/GaussSeg.cpp\ ./classes/gauss/GaussTria.cpp\ @@ -282,36 +314,4 @@ endif #}}} -#BAMG sources {{{ -if BAMG -issm_sources += ./bamg/BamgGeom.cpp\ - ./bamg/BamgMesh.cpp\ - ./bamg/BamgOpts.cpp\ - ./bamg/CrackedEdge.cpp\ - ./bamg/Curve.cpp\ - ./bamg/Direction.cpp\ - ./bamg/Edge.cpp\ - ./bamg/GeomEdge.cpp\ - ./bamg/GeomSubDomain.cpp\ - ./bamg/GeomVertex.cpp\ - ./bamg/Geometry.cpp\ - ./bamg/ListofIntersectionTriangles.cpp\ - ./bamg/EigenMetric.cpp\ - ./bamg/Metric.cpp\ - ./bamg/BamgQuadtree.cpp\ - ./bamg/SetOfE4.cpp\ - ./bamg/SubDomain.cpp\ - ./bamg/AdjacentTriangle.cpp\ - ./bamg/Triangle.cpp\ - ./bamg/BamgVertex.cpp\ - ./bamg/VertexOnEdge.cpp\ - ./bamg/VertexOnGeom.cpp\ - ./bamg/VertexOnVertex.cpp\ - ./bamg/Mesh.cpp\ - ./shared/Bamg/BigPrimeNumber.cpp\ - ./modules/Bamgx/Bamgx.cpp\ - ./modules/BamgConvertMeshx/BamgConvertMeshx.cpp\ - ./modules/BamgTriangulatex/BamgTriangulatex.cpp -endif -#}}} #Petsc sources {{{ if PETSC Index: /issm/trunk-jpl/src/c/bamg/BamgMesh.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgMesh.cpp (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/BamgMesh.cpp (revision 21623) @@ -8,5 +8,4 @@ this->EdgesSize[0]=0, this->EdgesSize[1]=0; this->Edges=NULL; this->TrianglesSize[0]=0, this->TrianglesSize[1]=0; this->Triangles=NULL; - this->QuadrilateralsSize[0]=0, this->QuadrilateralsSize[1]=0; this->Quadrilaterals=NULL; this->SubDomainsSize[0]=0, this->SubDomainsSize[1]=0; this->SubDomains=NULL; @@ -33,5 +32,4 @@ xDelete(this->Edges); xDelete(this->Triangles); - xDelete(this->Quadrilaterals); xDelete(this->SubDomains); Index: /issm/trunk-jpl/src/c/bamg/BamgMesh.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgMesh.h (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/BamgMesh.h (revision 21623) @@ -16,6 +16,4 @@ int TrianglesSize[2]; double* Triangles; - int QuadrilateralsSize[2]; - double* Quadrilaterals; int VerticesOnGeomVertexSize[2]; Index: /issm/trunk-jpl/src/c/bamg/BamgQuadtree.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgQuadtree.cpp (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/BamgQuadtree.cpp (revision 21623) @@ -7,46 +7,5 @@ namespace bamg { - - /*MACROS {{{*/ - /* - * - * J j - * ^ ^ - * | | +--------+--------+ - * | | | | | - * 1X | | | 2 | 3 | - * | | | | | - * | | +--------+--------+ - * | | | | | - * 0X | | | 0 | 1 | - * | | | | | - * | | +--------+--------+ - * | +-----------------------> i - * | - * |----------------------------> I - * X0 X1 - * - * box 0 -> I=0 J=0 IJ=00 = 0 - * box 1 -> I=1 J=0 IJ=01 = 1 - * box 2 -> I=0 J=1 IJ=10 = 2 - * box 3 -> I=1 J=1 IJ=11 = 3 - */ -#define INTER_SEG(a,b,x,y) (((y) > (a)) && ((x) <(b))) #define ABS(i) ((i)<0 ?-(i) :(i)) -#define MAX1(i,j) ((i)>(j) ?(i) :(j)) -#define NORM(i1,j1,i2,j2) MAX1(ABS((i1)-(j1)),ABS((i2)-(j2))) - - //IJ(i,j,l) returns the box number of i and j with respect to l - //if !j&l and !i&l -> 0 (box zero: lower left ) - //if !j&l and i&l -> 1 (box one: lower right) - //if j&l and !i&l -> 2 (box two: upper left ) - //if j&l and i&l -> 3 (box three:upper right) -#define IJ(i,j,l) ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 )) - - //I_IJ(k,l) returns l if first bit of k is 1, else 0 -#define I_IJ(k,l) ((k&1) ? l:0) - //J_IJ(k,l) returns l if second bit of k is 1, else 0 -#define J_IJ(k,l) ((k&2) ? l:0) - /*}}}*/ /*DOCUMENTATION What is a BamgQuadtree? {{{ * A Quadtree is a very simple way to group vertices according @@ -100,34 +59,41 @@ BamgQuadtree::BamgQuadtree(){/*{{{*/ - /*Number of boxes and vertices*/ - NbBamgQuadtreeBox=0; - NbVertices=0; + /*Initialize fields*/ + this->MaxDepth = 30; + this->MaxISize = 1073741824; //2^30 + this->MaxICoord = 1073741823; //2^30 - 1 + this->NbQuadtreeBox = 0; + this->NbVertices = 0; /*Create container*/ - boxcontainer=new DataSet(); + this->boxcontainer=new DataSet(); /*Create Root, pointer toward the main box*/ - root=NewBamgQuadtreeBox(); - - } - /*}}}*/ + this->root=NewBamgQuadtreeBox(); + + } /*}}}*/ BamgQuadtree::BamgQuadtree(Mesh * t,long nbv){ /*{{{*/ - /*Number of boxes and vertices*/ - NbBamgQuadtreeBox=0; - NbVertices=0; + /*Initialize fields*/ + this->MaxDepth = 30; + this->MaxISize = 1073741824; //2^30 + this->MaxICoord = 1073741823; //2^30 - 1 + this->NbQuadtreeBox = 0; + this->NbVertices = 0; /*Create container*/ - boxcontainer=new DataSet(); + this->boxcontainer=new DataSet(); /*Create Root, pointer toward the main box*/ - root=NewBamgQuadtreeBox(); + this->root=NewBamgQuadtreeBox(); /*Check Sizes*/ - _assert_(MaxISize>MaxICoor); + _assert_(this->MaxISize>this->MaxICoord); /*Add all vertices of the mesh*/ - if (nbv==-1) nbv=t->nbv; - for (int i=0;ivertices[i]); + if(nbv==-1) nbv=t->nbv; + for(int i=0;iAdd(t->vertices[i]); + } } @@ -140,5 +106,5 @@ /*Methods*/ - void BamgQuadtree::Add(BamgVertex &w){/*{{{*/ + void BamgQuadtree::Add(BamgVertex &w){/*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, BamgQuadtree.cpp/Add)*/ BamgQuadtreeBox** pb=NULL; @@ -164,5 +130,5 @@ //Get next subbox according to the bit value (level) - pb = &b->b[IJ(i,j,level)]; + pb = &b->box[BoxNumber(i,j,level)]; } @@ -195,5 +161,5 @@ /*Initialize the 4 pointers toward the 4 subboxes*/ - b->b[0]=b->b[1]=b->b[2]=b->b[3]=NULL; + b->box[0]=b->box[1]=b->box[2]=b->box[3]=NULL; /*level = 0010000 -> 0001000*/ @@ -205,8 +171,8 @@ int ij; /*bb is the new "sub"box of b where v4[k] is located*/ - BamgQuadtreeBox *bb = b->b[ij=IJ(v4[k]->i.x,v4[k]->i.y,level)]; + BamgQuadtreeBox *bb = b->box[ij=BoxNumber(v4[k]->i.x,v4[k]->i.y,level)]; // alloc the BamgQuadtreeBox - if (!bb) bb=b->b[ij]=NewBamgQuadtreeBox(); + if (!bb) bb=b->box[ij]=NewBamgQuadtreeBox(); /*Copy the current vertex*/ @@ -215,5 +181,5 @@ /*Get the subbox where w (i,j) is located*/ - pb = &b->b[IJ(i,j,level)]; + pb = &b->box[BoxNumber(i,j,level)]; } @@ -228,48 +194,79 @@ } /*}}}*/ - BamgVertex* BamgQuadtree::NearestVertex(Icoor1 i,Icoor1 j) {/*{{{*/ - /*Original code from Frederic Hecht (BAMG v1.01, BamgQuadtree.cpp/NearestVertex)*/ - - /*Intermediaries*/ - BamgQuadtreeBox *pb[MaxDepth]; - int pi[MaxDepth]; - Icoor1 ii[MaxDepth]; - Icoor1 jj[MaxDepth]; - int level; - long n0; - BamgQuadtreeBox *b; - long h0; - long h = MaxISize; - long hb= MaxISize; - Icoor1 i0=0,j0=0; + int BamgQuadtree::BoxNumber(int i,int j,int l){/*{{{*/ + /* + * + * J j + * ^ ^ + * | | +--------+--------+ + * | | | | | + * 1X | | | 2 | 3 | + * | | | | | + * | | +--------+--------+ + * | | | | | + * 0X | | | 0 | 1 | + * | | | | | + * | | +--------+--------+ + * | +-----------------------> i + * | + * |----------------------------> I + * X0 X1 + * + * box 0 -> I=0 J=0 BoxNumber=00 = 0 + * box 1 -> I=1 J=0 BoxNumber=01 = 1 + * box 2 -> I=0 J=1 BoxNumber=10 = 2 + * box 3 -> I=1 J=1 BoxNumber=11 = 3 + */ + //BoxNumber(i,j,l) returns the box number of i and j with respect to l + //if !j&l and !i&l -> 0 (box zero: lower left ) + //if !j&l and i&l -> 1 (box one: lower right) + //if j&l and !i&l -> 2 (box two: upper left ) + //if j&l and i&l -> 3 (box three:upper right) + return ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 )); + }/*}}}*/ + bool BamgQuadtree::Intersection(int a,int b,int x,int y){/*{{{*/ + /*is [x y] in [a b]*/ + return ((y) > (a)) && ((x) <(b)); + }/*}}}*/ + BamgVertex* BamgQuadtree::NearestVertex(int xi,int yi) {/*{{{*/ /*initial output as NULL (no vertex found)*/ BamgVertex* nearest_v=NULL; - /*Project w coordinates (i,j) onto [0,MaxISize-1] x [0,MaxISize-1] -> (iplus,jplus)*/ - Icoor1 iplus( inbitems) return nearest_v; + if(!this->root->nbitems) return nearest_v; + + /*Project coordinates (xi,yi) onto [0,MaxICoord-1] x [0,MaxICoord-1]*/ + int xi2 = xi; + int yi2 = yi; + if(xi<0) xi2 = 0; + if(xi>MaxISize) xi2 = MaxICoord; + if(yi<0) yi2 = 0; + if(yi>MaxISize) yi2 = MaxICoord; + + /*Get inital box (the largest)*/ + BamgQuadtreeBox* b = this->root; + + /*Initialize level and sizes for largest box*/ + int levelbin = (1L<MaxDepth);// = 2^30 + int h = this->MaxISize; + int hb = this->MaxISize; + int i0 = 0; + int j0 = 0; /*else, find the smallest non-empty BamgQuadtreeBox containing the point (i,j)*/ - while((n0=b->nbitems)<0){ - - Icoor1 hb2 = hb >> 1; //size of the current box - int k = IJ(iplus,jplus,hb2); //box number (0,1,2 or 3) - BamgQuadtreeBox *b0 = b->b[k]; //pointer toward current box - - /* break if NULL box or empty (Keep previous box b)*/ - if (( b0 == NULL) || (b0->nbitems == 0)) break; + while(b->nbitems<0){ + + int hb2 = hb >> 1; //size of the current box + int k = BoxNumber(xi2,yi2,hb2); //box number (0,1,2 or 3) + BamgQuadtreeBox* b0 = b->box[k]; //pointer toward current box + + /* break if box is empty (Keep previous box b)*/ + if((b0==NULL) || (b0->nbitems==0)) break; /*Get next Quadtree box*/ - b=b0; - i0 += I_IJ(k,hb2); // i orign of BamgQuadtreeBox (macro) - j0 += J_IJ(k,hb2); // j orign of BamgQuadtreeBox - hb = hb2; // size of the box (in Int) + b = b0; + hb = hb2; + this->SubBoxCoords(&i0,&j0,k,hb2); } @@ -280,18 +277,18 @@ * - i0: x coordinate of the lower left corner * - j0: y coordinate of the lower left corner*/ + int n0 = b->nbitems; /* if the current subbox is holding vertices, we are almost done*/ - if (n0>0){ - //loop over the vertices of the box and find the closest vertex + if(n0>0){ + /*loop over the vertices of the box and find the closest vertex*/ for(int k=0;kv[k]->i; - - /*Compute norm with w*/ - h0=NORM(iplus,i2.x,jplus,i2.y); + int xiv = b->v[k]->i.x; + int yiv = b->v[k]->i.y; + + + int h0 = Norm(xi2,xiv,yi2,yiv); /*is it smaller than previous value*/ - if (h0v[k]; @@ -306,4 +303,8 @@ /*Initialize search variables*/ + BamgQuadtreeBox **pb = xNew(this->MaxDepth); + int* pi = xNew(this->MaxDepth); + int* ii = xNew(this->MaxDepth); + int* jj = xNew(this->MaxDepth); pb[0]=b; //pointer toward the box b pi[0]=b->nbitems>0?(int)b->nbitems:4;//number of boxes in b @@ -315,12 +316,11 @@ /*Main loop*/ - level=0; - do { - + int level=0; + do{ /*get current box*/ - b= pb[level]; + b = pb[level]; /*Loop over the items in current box (if not empty!)*/ - while (pi[level]){ + while(pi[level]){ /*We are looping now over the items of b. k is the current index (in [0 3])*/ @@ -330,6 +330,5 @@ /*if the current subbox is holding vertices (b->nbitems<0 is subboxes)*/ if (b->nbitems>0){ - I2 i2=b->v[k]->i; - h0 = NORM(iplus,i2.x,jplus,i2.y); + int h0 = Norm(xi2,b->v[k]->i.x,yi2,b->v[k]->i.y); if (h0b[k])){ + if((b=b->box[k])){ /*Get size (hb) and coordinates of the current sub-box lowest left corner*/ hb>>=1; - Icoor1 iii = ii[level]+I_IJ(k,hb); - Icoor1 jjj = jj[level]+J_IJ(k,hb); - - /*if the current point (iplus,jplus) is in b (modulo h), this box is good: + int iii = ii[level]; + int jjj = jj[level]; + this->SubBoxCoords(&iii,&jjj,k,hb); + + /*if the current point (xi2,yi2) is in b (modulo h), this box is good: * it is holding vertices that are close to w */ - if (INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)){ + if(this->Intersection(iii,iii+hb,xi2-h,xi2+h) && this->Intersection(jjj,jjj+hb,yi2-h,yi2+h)){ level++; - pb[level]= b; - pi[level]= b->nbitems>0 ?(int) b->nbitems : 4 ; - ii[level]= iii; - jj[level]= jjj; + pb[level] = b; + pi[level] = b->nbitems>0 ? b->nbitems:4 ; + ii[level] = iii; + jj[level] = jjj; } - //else go backwards + /*else go backwards*/ else{ - //shifted righ by one bit: hb=001000000 -> 01000000 + /*shifted righ by one bit: hb=001000000 -> 01000000*/ b=b0; hb<<=1; @@ -379,82 +379,107 @@ * in case there is a vertex closest to w that has not yet been tested*/ hb <<= 1; - } while (level--); - - /*return nearest_v, nearest vertex*/ + }while(level--); + + /*return nearest vertex pointer*/ + xDelete(pb); + xDelete(pi); + xDelete(ii); + xDelete(jj); return nearest_v; - } /*}}}*/ - BamgQuadtree::BamgQuadtreeBox* BamgQuadtree::NewBamgQuadtreeBox(void){/*{{{*/ - - /*Output*/ - BamgQuadtreeBox* newbox=NULL; - - /*Create and initialize a new box*/ - newbox=new BamgQuadtreeBox; - newbox->nbitems=0; - newbox->b[0]=NULL; - newbox->b[1]=NULL; - newbox->b[2]=NULL; - newbox->b[3]=NULL; - - /*Add root to the container*/ - boxcontainer->AddObject(newbox); - - /*Increase counter*/ - NbBamgQuadtreeBox++; - - /*currentbox now points toward next quadtree box*/ - return newbox; + int BamgQuadtree::Norm(int xi1,int xi2,int yi1,int yi2){/*{{{*/ + + int deltax = xi2 - xi1; + int deltay = yi2 - yi1; + + if(deltax<0) deltax = -deltax; + if(deltay<0) deltay = -deltay; + + if(deltax> deltay){ + return deltax; + } + else{ + return deltay; + } }/*}}}*/ - BamgVertex* BamgQuadtree::ToClose(BamgVertex & v,double seuil,Icoor1 hx,Icoor1 hy){/*{{{*/ - /*Original code from Frederic Hecht (BAMG v1.01, BamgQuadtree.cpp/ToClose)*/ - - const Icoor1 i=v.i.x; - const Icoor1 j=v.i.y; - const R2 X(v.r); - const Metric Mx(v.m); - - BamgQuadtreeBox * pb[ MaxDepth ]; - int pi[ MaxDepth ]; - Icoor1 ii[ MaxDepth ], jj [ MaxDepth]; + void BamgQuadtree::SubBoxCoords(int* pi,int*pj,int boxnumber,int length){/*{{{*/ + /* + * j (first bit) + * ^ + * | +--------+--------+ + * | | | | + * 1X | | 2 | 3 | + * | | | | + * | +--------+--------+ + * | | | | + * 0X | | 0 | 1 | + * | | | | + * | +--------+--------+ + * +-----------------------> i (second bit) + * X0 X1 + */ + + /*Add sub-box coordinate to i and j*/ + //_assert_(!(*pi & length)); + //_assert_(!(*pj & length)); + + /*length if first bit of boxnumber is 1, elengthse 0*/ + *pi += ((boxnumber & 1) ? length:0); + /*length if second bit of boxnumber is 1, elengthse 0*/ + *pj += ((boxnumber & 2) ? length:0); + + }/*}}}*/ + BamgVertex* BamgQuadtree::TooClose(BamgVertex* v,double threshold,int hx,int hy){/*{{{*/ + + const int i=v->i.x; + const int j=v->i.y; + const double Xx = v->r.x; + const double Xy = v->r.y; + Metric* Mx = new Metric(v->m); + + BamgQuadtreeBox *pb[MaxDepth]; + int pi[MaxDepth]; + int ii[MaxDepth], jj[MaxDepth]; int l=0; // level - BamgQuadtreeBox * b; - Icoor1 hb = MaxISize; - Icoor1 i0=0,j0=0; - - // BamgVertex *vn=0; - - if (!root->nbitems) - return 0; // empty tree - - // general case ----- + BamgQuadtreeBox* b; + int hb = MaxISize; + int i0=0,j0=0; + + // BamgVertex *vn=0; + + if(!root->nbitems) return 0; // empty tree + + // general case pb[0]=root; - pi[0]=root->nbitems>0 ?(int) root->nbitems : 4 ; + pi[0]=root->nbitems>0 ?(int)root->nbitems:4; ii[0]=i0; jj[0]=j0; - do { + do{ b= pb[l]; - while (pi[l]--){ + while(pi[l]--){ int k = pi[l]; - if (b->nbitems>0){ // BamgVertex BamgQuadtreeBox none empty - I2 i2 = b->v[k]->i; - if ( ABS(i-i2.x) v[k]->r); - if(LengthInterpole(Mx(XY), b->v[k]->m(XY)) < seuil){ + if(b->nbitems>0){ // BamgVertex BamgQuadtreeBox none empty + int i2x = b->v[k]->i.x; + int i2y = b->v[k]->i.y; + if (ABS(i-i2x)v[k]->r.x - Xx; + double XYy = b->v[k]->r.y - Xy; + if(LengthInterpole(Mx->Length(XYx,XYy),b->v[k]->m.Length(XYx,XYy)) < threshold){ + delete Mx; return b->v[k]; } - } + } } else{ // Pointer BamgQuadtreeBox BamgQuadtreeBox *b0=b; - if ((b=b->b[k])){ + if ((b=b->box[k])){ hb >>=1 ; // div by 2 - long iii = ii[l]+I_IJ(k,hb); - long jjj = jj[l]+J_IJ(k,hb); - - if (INTER_SEG(iii,iii+hb,i-hx,i+hx) && INTER_SEG(jjj,jjj+hb,j-hy,j+hy)){ + int iii = ii[l]; + int jjj = jj[l]; + this->SubBoxCoords(&iii,&jjj,k,hb); + + if(this->Intersection(iii,iii+hb,i-hx,i+hx) && this->Intersection(jjj,jjj+hb,j-hy,j+hy)){ pb[++l]= b; pi[l]= b->nbitems>0 ?(int) b->nbitems : 4 ; @@ -474,8 +499,30 @@ } hb <<= 1; // mul by 2 - } while (l--); - + }while(l--); + + delete Mx; return 0; } /*}}}*/ + + BamgQuadtree::BamgQuadtreeBox* BamgQuadtree::NewBamgQuadtreeBox(void){/*{{{*/ + + /*Output*/ + BamgQuadtreeBox* newbox=NULL; + + /*Create and initialize a new box*/ + newbox=new BamgQuadtreeBox; + newbox->nbitems=0; + newbox->box[0]=NULL; + newbox->box[1]=NULL; + newbox->box[2]=NULL; + newbox->box[3]=NULL; + + /*Add root to the container*/ + boxcontainer->AddObject(newbox); + this->NbQuadtreeBox++; + + /*currentbox now points toward next quadtree box*/ + return newbox; + }/*}}}*/ } Index: /issm/trunk-jpl/src/c/bamg/BamgQuadtree.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgQuadtree.h (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/BamgQuadtree.h (revision 21623) @@ -5,12 +5,7 @@ #include "./include.h" #include "../datastructures/datastructures.h" +class BamgVertex; namespace bamg { - - const int MaxDepth = 30; - const long MaxISize = ( 1L << MaxDepth); // = 2^30 : 010000000000..000 (bitwise operation) - - class BamgVertex; - class BamgQuadtree{ @@ -24,9 +19,7 @@ class BamgQuadtreeBox: public Object{ public: - int nbitems; // number of current vertices in the box - union{ - BamgQuadtreeBox* b[4]; - BamgVertex* v[4]; - }; + int nbitems; + BamgQuadtreeBox *box[4]; + BamgVertex *v[4]; /*Object functions*/ void Echo() {_error_("not implemented yet"); }; @@ -44,7 +37,10 @@ /*BamgQuadtree public Fields*/ - BamgQuadtreeBox* root; - long NbBamgQuadtreeBox; - long NbVertices; + int MaxDepth; // maximum number of subdivision + int MaxICoord; // maximum integer coordinate 2^MaxDepth -1 + int MaxISize; // maximum integer coordinate 2^MaxDepth + BamgQuadtreeBox *root; // main box + long NbQuadtreeBox; // total number of boxes + long NbVertices; // number of points BamgQuadtree(); @@ -52,8 +48,12 @@ ~BamgQuadtree(); - BamgVertex *NearestVertex(Icoor1 i,Icoor1 j); + void Add(BamgVertex &w); + int BoxNumber(int i,int j,int l); + bool Intersection(int a,int b,int x,int y); + BamgVertex *NearestVertex(int i,int j); BamgQuadtreeBox *NewBamgQuadtreeBox(void); - BamgVertex *ToClose(BamgVertex &,double ,Icoor1,Icoor1); - void Add(BamgVertex &w); + int Norm(int xi1,int yi1,int xi2,int yi2); + void SubBoxCoords(int* pi,int*pj,int boxcoord,int length); + BamgVertex *TooClose(BamgVertex*,double,int,int); }; } Index: /issm/trunk-jpl/src/c/bamg/BamgVertex.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgVertex.cpp (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/BamgVertex.cpp (revision 21623) @@ -115,5 +115,5 @@ ret = t->Optim(IndexInTriangle,koption); if(!i){ - t =0; // for no future optime + t =0; // for no future optim IndexInTriangle= 0; } @@ -196,14 +196,5 @@ BamgVertex *v0,*v1,*v2,*v3; if (tria->det<0) ok =1; - else if (tria->Quadrangle(v0,v1,v2,v3)) - { - vP = vPsave; - double qold =QuadQuality(*v0,*v1,*v2,*v3); - vP = vPnew; - double qnew =QuadQuality(*v0,*v1,*v2,*v3); - if (qnewdet < detold/2 ) ok=1; - } tria->SetUnMarkUnSwap(0); @@ -225,33 +216,3 @@ } /*}}}*/ - - /*Intermediary*/ - double QuadQuality(const BamgVertex & a,const BamgVertex &b,const BamgVertex &c,const BamgVertex &d) {/*{{{*/ - /*Original code from Frederic Hecht (BAMG v1.01, MeshQuad.cpp/QuadQuality)*/ - - // calcul de 4 angles -- - R2 A((R2)a),B((R2)b),C((R2)c),D((R2)d); - R2 AB(B-A),BC(C-B),CD(D-C),DA(A-D); - // Move(A),Line(B),Line(C),Line(D),Line(A); - const Metric & Ma = a; - const Metric & Mb = b; - const Metric & Mc = c; - const Metric & Md = d; - - double lAB=Norme2(AB); - double lBC=Norme2(BC); - double lCD=Norme2(CD); - double lDA=Norme2(DA); - AB /= lAB; BC /= lBC; CD /= lCD; DA /= lDA; - // version aniso - double cosDAB= Ma(DA,AB)/(Ma(DA)*Ma(AB)),sinDAB= Det(DA,AB); - double cosABC= Mb(AB,BC)/(Mb(AB)*Mb(BC)),sinABC= Det(AB,BC); - double cosBCD= Mc(BC,CD)/(Mc(BC)*Mc(CD)),sinBCD= Det(BC,CD); - double cosCDA= Md(CD,DA)/(Md(CD)*Md(DA)),sinCDA= Det(CD,DA); - double sinmin=Min(Min(sinDAB,sinABC),Min(sinBCD,sinCDA)); - if (sinmin<=0) return sinmin; - return 1.0-Max(Max(Abs(cosDAB),Abs(cosABC)),Max(Abs(cosBCD),Abs(cosCDA))); - } - /*}}}*/ - } Index: /issm/trunk-jpl/src/c/bamg/BamgVertex.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgVertex.h (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/BamgVertex.h (revision 21623) @@ -4,5 +4,4 @@ #include "./include.h" #include "./Metric.h" -#include "./Direction.h" #include "./BamgOpts.h" @@ -24,5 +23,4 @@ long ReferenceNumber; long PreviousNumber; - Direction DirOfSearch; short IndexInTriangle; // the vertex number in triangle; varies between 0 and 2 in t @@ -40,5 +38,5 @@ operator const R2 & () const {return r;} // Cast operator operator Metric () const {return m;} // Cast operator - double operator()(R2 x) const { return m(x);} // Get x in the metric m + double operator()(R2 x) const { return m.Length(x.x,x.y);} // Get x in the metric m /*methods (No constructor and no destructors...)*/ @@ -54,7 +52,4 @@ inline void Set(const BamgVertex &rec,const Mesh & ,Mesh & ){*this=rec;} }; - - //Intermediary - double QuadQuality(const BamgVertex &,const BamgVertex &,const BamgVertex &,const BamgVertex &); } #endif Index: sm/trunk-jpl/src/c/bamg/Direction.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/Direction.cpp (revision 21622) +++ (revision ) @@ -1,33 +1,0 @@ -#include -#include -#include -#include - -#include "Direction.h" - -namespace bamg { - - /*Constructors/Destructors*/ - Direction::Direction():/*{{{*/ - dir(MaxICoor){ - - }/*}}}*/ - Direction::Direction(Icoor1 i,Icoor1 j) {/*{{{*/ - Icoor2 n2 = 2*(Abs(i)+Abs(j)); - Icoor2 r = MaxICoor* (Icoor2) i; - Icoor1 r1 = (Icoor1) (2*(r/ n2)); // odd number - dir = (j>0) ? r1 : r1+1; // odd-> j>0 even-> j<0 - }/*}}}*/ - - /*Methods*/ - int Direction::direction(Icoor1 i,Icoor1 j) {/*{{{*/ - int r =1; - if (dir!= MaxICoor) { - Icoor2 x(dir/2),y1(MaxICoor/2-Abs(x)),y((dir%2)?-y1:y1); - r = (x*i + y*j) >=0; - } - return r; - } - /*}}}*/ - -} Index: sm/trunk-jpl/src/c/bamg/Direction.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/Direction.h (revision 21622) +++ (revision ) @@ -1,20 +1,0 @@ -#ifndef _DIRECTION_H_ -#define _DIRECTION_H_ - -#include "./include.h" -#include "../shared/Bamg/shared.h" - -namespace bamg { - - class Direction { - private: - Icoor1 dir; - - public: - //Methods - Direction(); - Direction(Icoor1 i,Icoor1 j); - int direction(Icoor1 i,Icoor1 j); - }; -} -#endif Index: sm/trunk-jpl/src/c/bamg/DoubleAndInt.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/DoubleAndInt.h (revision 21622) +++ (revision ) @@ -1,19 +1,0 @@ -#ifndef _DOUBLEANDINT_H_ -#define _DOUBLEANDINT_H_ - -#include "./include.h" - -namespace bamg { - - class DoubleAndInt { - //class used by Mesh::MakeQuadrangles - - public: - double q; - long i3j; - - //Operators - int operator<(DoubleAndInt a){return q > a.q;} - }; -} -#endif Index: /issm/trunk-jpl/src/c/bamg/EigenMetric.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/EigenMetric.cpp (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/EigenMetric.cpp (revision 21623) @@ -37,6 +37,6 @@ lambda1=0; lambda2=0; - v.x=1; - v.y=0; + vx=1; + vy=0; } /*2: delta is small -> double root*/ @@ -44,6 +44,6 @@ lambda1=-b/2; lambda2=-b/2; - v.x=1; - v.y=0; + vx=1; + vy=0; } /*3: general case -> two roots*/ @@ -76,11 +76,11 @@ if (norm2lambda1 = r1; + this->lambda2 = r2; + this->vx = vp1.x; + this->vy = vp1.y; }/*}}}*/ @@ -104,6 +107,6 @@ _printf_(" lambda1: " << lambda1 << "\n"); _printf_(" lambda2: " << lambda2 << "\n"); - _printf_(" v.x: " << v.x << "\n"); - _printf_(" v.y: " << v.y << "\n"); + _printf_(" vx: " << vx << "\n"); + _printf_(" vy: " << vy << "\n"); return; Index: /issm/trunk-jpl/src/c/bamg/Geometry.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/Geometry.cpp (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/Geometry.cpp (revision 21623) @@ -8,6 +8,4 @@ namespace bamg { - - static const Direction NoDirOfSearch=Direction(); /*Constructors/Destructors*/ @@ -80,5 +78,4 @@ vertices[i].r.y=(double)bamggeom->Vertices[i*3+1]; vertices[i].ReferenceNumber=(long)bamggeom->Vertices[i*3+2]; - vertices[i].DirOfSearch=NoDirOfSearch; vertices[i].color =0; vertices[i].type=0; @@ -772,4 +769,5 @@ /*Original code from Frederic Hecht (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/ /*Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge*/ + /*FIXME: should go away*/ double save_s=s; @@ -807,5 +805,5 @@ //Compute metric of new vertex as a linear interpolation of the two others - V.m=Metric(1.0-s,v0,s,v1); + V.m=Metric(1.0-s,v0.m,s,v1.m); const int mxe=100; Index: /issm/trunk-jpl/src/c/bamg/Mesh.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/Mesh.cpp (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/Mesh.cpp (revision 21623) @@ -9,6 +9,4 @@ namespace bamg { - - static const Direction NoDirOfSearch=Direction(); /*Constructors/Destructors*/ @@ -167,5 +165,4 @@ nbsubdomains = Th.nbsubdomains; nbtout = Th.nbtout; - nbq = Th.nbq ; NbVerticesOnGeomVertex = Th.NbVerticesOnGeomVertex; if(NbVerticesOnGeomVertex) @@ -269,5 +266,4 @@ vertices[i].r.y=y[i]; vertices[i].ReferenceNumber=1; - vertices[i].DirOfSearch =NoDirOfSearch; vertices[i].m=M1; vertices[i].color=0; @@ -336,5 +332,4 @@ vertices[i].r.y=bamgmesh->Vertices[i*3+1]; vertices[i].ReferenceNumber=(long)bamgmesh->Vertices[i*3+2]; - vertices[i].DirOfSearch =NoDirOfSearch; vertices[i].m=M1; vertices[i].color=0; @@ -362,26 +357,4 @@ else{ if(verbose>5) _error_("no Triangles found in the initial mesh"); - } - - //Quadrilaterals - if(bamgmesh->Quadrilaterals){ - if(verbose>5) _printf_(" processing Quadrilaterals\n"); - long i1,i2,i3,i4; - triangles =new Triangle[nbt]; - for (i=0;iQuadrilateralsSize[0];i++){ - //divide the quad into two triangles - Triangle & t1 = triangles[2*i]; - Triangle & t2 = triangles[2*i+1]; - i1=(long)bamgmesh->Quadrilaterals[i*5+0]-1; //for C indexing - i2=(long)bamgmesh->Quadrilaterals[i*5+1]-1; //for C indexing - i3=(long)bamgmesh->Quadrilaterals[i*5+2]-1; //for C indexing - i4=(long)bamgmesh->Quadrilaterals[i*5+3]-1; //for C indexing - t1=Triangle(this,i1,i2,i3); - t2=Triangle(this,i3,i4,i1); - t1.color=(long)bamgmesh->Quadrilaterals[i*5+4]; - t2.color=(long)bamgmesh->Quadrilaterals[i*5+4]; - t1.SetHidden(OppositeEdge[1]); // two times because the adj was not created - t2.SetHidden(OppositeEdge[1]); - } } @@ -698,5 +671,5 @@ /*Triangles*/ if(verbose>5) _printf_(" writing Triangles\n"); - k=nbInT-nbq*2; + k=nbInT; num=0; bamgmesh->TrianglesSize[0]=k; @@ -713,25 +686,4 @@ bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ReferenceNumber; num=num+1; - } - } - } - - /*Quadrilaterals*/ - if(verbose>5) _printf_(" writing Quadrilaterals\n"); - bamgmesh->QuadrilateralsSize[0]=nbq; - bamgmesh->QuadrilateralsSize[1]=5; - if (nbq){ - bamgmesh->Quadrilaterals=xNew(5*nbq); - for (i=0;iQuadrilaterals[i*5+0]=GetId(v0)+1; //back to M indexing - bamgmesh->Quadrilaterals[i*5+1]=GetId(v1)+1; //back to M indexing - bamgmesh->Quadrilaterals[i*5+2]=GetId(v2)+1; //back to M indexing - bamgmesh->Quadrilaterals[i*5+3]=GetId(v3)+1; //back to M indexing - bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]].ReferenceNumber; } } @@ -1113,5 +1065,5 @@ // s0 tt2 s1 //------------------------------- - + /*Intermediaries*/ Triangle* tt[3]; //the three triangles @@ -1230,6 +1182,6 @@ } } - } - /*}}}*/ + + }/*}}}*/ void Mesh::BoundAnisotropy(double anisomax,double hminaniso) {/*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Metric.cpp/BoundAnisotropy)*/ @@ -2371,15 +2323,14 @@ _printf_(" nbt = " << nbt << "\n"); _printf_(" nbe = " << nbe << "\n"); - _printf_(" nbq = " << nbq << "\n"); _printf_(" index:\n"); for (i=0;iToClose(vi,seuil,hi,hj)){ + if(!quadtree->TooClose(&vi,seuil,hi,hj)){ // a good new point BamgVertex &vj = vertices[iv]; @@ -2999,45 +2949,4 @@ } /*}}}*/ - void Mesh::MakeQuadrangles(double costheta){/*{{{*/ - /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/MakeQuadrangles)*/ - - long int verbose=0; - - if (verbose>2) _printf_("MakeQuadrangles costheta = " << costheta << "\n"); - - if (costheta >1) { - if (verbose>5) _printf_(" do nothing: costheta > 1\n"); - } - - long nbqq = (nbt*3)/2; - DoubleAndInt *qq = new DoubleAndInt[nbqq]; - - long i,ij; - int j; - long k=0; - for (i=0;i=costheta) - qq[k++].i3j=i*3+j; - // sort qq - HeapSort(qq,k); - - long kk=0; - for (ij=0;ij=costheta) - triangles[i].SetHidden(j),kk++; - } - nbq = kk; - if (verbose>2){ - _printf_(" number of quadrilaterals = " << nbq << "\n"); - _printf_(" number of triangles = " << nbt-nbtout- nbq*2 << "\n"); - _printf_(" number of outside triangles = " << nbtout << "\n"); - } - delete [] qq; - } - /*}}}*/ void Mesh::MakeBamgQuadtree() { /*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/MakeBamgQuadtree)*/ @@ -3073,5 +2982,5 @@ lmax = Max(lmax,l); if(l> maxsubdiv2){ - R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M + R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M double lc = M(AC,AC); D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; @@ -3086,5 +2995,5 @@ lmax = Max(lmax,l); if(l> maxsubdiv2){ - R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M + R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M double lc = M(AC,AC); D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; @@ -3857,8 +3766,7 @@ if (tstart[i] != &vide) // not a boundary vertex delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega)); - if (!nbq) - for ( i=0;i3) _printf_(" move max = " << pow(delta,0.5) << ", iteration = " << k << ", nb of swap = " << NbSwap << "\n"); } @@ -3914,6 +3822,6 @@ double ll = Norme2(Aij); if (0) { - double hi = ll/vi.m(Aij); - double hj = ll/vj.m(Aij); + double hi = ll/vi.m.Length(Aij.x,Aij.y); + double hj = ll/vj.m.Length(Aij.x,Aij.y); if(hi < hj) { @@ -3928,5 +3836,5 @@ else { - double li = vi.m(Aij); + double li = vi.m.Length(Aij.x,Aij.y); if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) ) if(first_np_or_next_t1[j]<0) // if the metrix change @@ -3973,5 +3881,4 @@ vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m); vertices[nbv].ReferenceNumber=0; - vertices[nbv].DirOfSearch = NoDirOfSearch ; } NbSplitEdge++; @@ -3992,5 +3899,4 @@ // a good new point vi.ReferenceNumber=0; - vi.DirOfSearch =NoDirOfSearch; Triangle *tcvi = TriangleFindFromCoord(vi.i,det3); if (tcvi && !tcvi->link) { @@ -4181,5 +4087,4 @@ vertices[i].r.y=y[i]; vertices[i].ReferenceNumber=1; - vertices[i].DirOfSearch =NoDirOfSearch; vertices[i].m=M1; vertices[i].color=0; @@ -4331,5 +4236,5 @@ Metric MA = background ? BTh.MetricAt(a->r) :a->m ; //Get metric associated to A Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B - double ledge = (MA(AB) + MB(AB))/2; //Get edge length in metric + double ledge = (MA.Length(AB.x,AB.y) + MB.Length(AB.x,AB.y))/2; //Get edge length in metric /* We are now creating the mesh edges from the geometrical edge selected above. @@ -4367,5 +4272,5 @@ MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB); AB = A-B; - lSubEdge[kk]=(ledge+=(MAs(AB)+MBs(AB))/2); + lSubEdge[kk]=(ledge+=(MAs.Length(AB.x,AB.y)+MBs.Length(AB.x,AB.y))/2); } } @@ -4412,5 +4317,4 @@ vb->m = Metric(aa,a->m,bb,b->m); vb->ReferenceNumber = e->ReferenceNumber; - vb->DirOfSearch =NoDirOfSearch; double abcisse = k ? bb : aa; vb->r = e->F(abcisse); @@ -4735,5 +4639,4 @@ ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); A1->ReferenceNumber = eeequi.ReferenceNumber; - A1->DirOfSearch =NoDirOfSearch; e->GeomEdgeHook = ongequi; e->v[0]=A0; Index: /issm/trunk-jpl/src/c/bamg/Mesh.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/Mesh.h (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/Mesh.h (revision 21623) @@ -33,5 +33,5 @@ long NbRef; // counter of ref on the this class if 0 we can delete long maxnbv,maxnbt; // nombre max de sommets , de triangles - long nbv,nbt,nbe,nbq; // nb of vertices, of triangles, of edges and quadrilaterals + long nbv,nbt,nbe; // nb of vertices, of triangles and edges long nbsubdomains; long nbtout; // Nb of oudeside triangle @@ -86,9 +86,8 @@ void SmoothMetric(double raisonmax) ; void BoundAnisotropy(double anisomax,double hminaniso= 1e-100) ; - void MaxSubDivision(double maxsubdiv); Edge** MakeGeomEdgeToEdge(); long SplitInternalEdgeWithBorderVertices(); - void MakeQuadrangles(double costheta); void MakeBamgQuadtree(); + void MaxSubDivision(double maxsubdiv); void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1); long InsertNewPoints(long nbvold,long & NbTSwap,bool random); Index: /issm/trunk-jpl/src/c/bamg/Metric.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/Metric.cpp (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/Metric.cpp (revision 21623) @@ -13,32 +13,49 @@ /*Constructor/Destructor*/ - Metric::Metric(double a): a11(1/(a*a)),a21(0),a22(1/(a*a)){/*{{{*/ + Metric::Metric(double a){/*{{{*/ + + /*Isotropic metric for a length of a as unit*/ + this->a11 = 1./(a*a); + this->a21 = 0.; + this->a22 = 1./(a*a); }/*}}}*/ - Metric::Metric(double a,double b,double c) :a11(a),a21(b),a22(c){/*{{{*/ + Metric::Metric(double a11_in,double a21_in,double a22_in){/*{{{*/ + + this->a11 = a11_in; + this->a21 = a21_in; + this->a22 = a22_in; }/*}}}*/ - Metric::Metric(const double a[3],const Metric& m0, const Metric& m1,const Metric& m2 ){/*{{{*/ - /*Original code from Frederic Hecht (BAMG v1.01, Metric.cpp/Metric)*/ - + Metric::Metric(const double a[3],const Metric m0,const Metric m1,const Metric m2 ){/*{{{*/ + + /*Create metric using 3 inputs*/ Metric mab(a[0]*m0.a11 + a[1]*m1.a11 + a[2]*m2.a11, a[0]*m0.a21 + a[1]*m1.a21 + a[2]*m2.a21, a[0]*m0.a22 + a[1]*m1.a22 + a[2]*m2.a22); + /*Convert to eigen metric*/ EigenMetric vab(mab); - - R2 v1(vab.v.x,vab.v.y); - R2 v2(-v1.y,v1.x); - - double h1 = a[0] / m0(v1) + a[1] / m1(v1) + a[2] / m2(v1); - double h2 = a[0] / m0(v2) + a[1] / m1(v2) + a[2] / m2(v2); - - vab.lambda1 = 1 / (h1*h1); - vab.lambda2 = 1 / (h2*h2); - *this = vab; - } - /*}}}*/ - Metric::Metric(double a,const Metric& ma, double b,const Metric& mb) { /*{{{*/ - /*Original code from Frederic Hecht (BAMG v1.01, Metric.cpp/EigenMetric)*/ + double v1x = + vab.vx; + double v1y = + vab.vy; + double v2x = - vab.vy; + double v2y = + vab.vx; + + double h1=a[0] / m0.Length(v1x,v1y) + a[1] / m1.Length(v1x,v1y) + a[2] / m2.Length(v1x,v1y); + double h2=a[0] / m0.Length(v2x,v2y) + a[1] / m1.Length(v2x,v2y) + a[2] / m2.Length(v2x,v2y); + + vab.lambda1 = 1. / (h1*h1); + vab.lambda2 = 1. / (h2*h2); + + /*Build metric from vab*/ + double v00=vab.vx*vab.vx; + double v11=vab.vy*vab.vy; + double v01=vab.vx*vab.vy; + this->a11=v00*vab.lambda1+v11*vab.lambda2; + this->a21=v01*(vab.lambda1-vab.lambda2); + this->a22=v00*vab.lambda2+v11*vab.lambda1; + + } /*}}}*/ + Metric::Metric(double a,const Metric ma,double b,const Metric mb) { /*{{{*/ /*Compute metric (linear combination of ma and mb)*/ @@ -47,13 +64,22 @@ /*Get Eigen values and vectors*/ EigenMetric vab(mab); - R2 v1(vab.v.x,vab.v.y); - R2 v2(-v1.y,v1.x); + double v1x = + vab.vx; + double v1y = + vab.vy; + double v2x = - vab.vy; + double v2y = + vab.vx; /*Modify eigen values (a+b=1)*/ - double h1 = a/ma(v1) + b/mb(v1); - double h2 = a/ma(v2) + b/mb(v2); - vab.lambda1 = 1/(h1*h1); - vab.lambda2 = 1/(h2*h2); - *this=vab; + double h1 = a/ma.Length(v1x,v1y) + b/mb.Length(v1x,v1y); + double h2 = a/ma.Length(v2x,v2y) + b/mb.Length(v2x,v2y); + vab.lambda1 = 1./(h1*h1); + vab.lambda2 = 1./(h2*h2); + + /*Build metric from vab*/ + double v00=vab.vx*vab.vx; + double v11=vab.vy*vab.vy; + double v01=vab.vx*vab.vy; + this->a11=v00*vab.lambda1+v11*vab.lambda2; + this->a21=v01*(vab.lambda1-vab.lambda2); + this->a22=v00*vab.lambda2+v11*vab.lambda1; } /*}}}*/ @@ -123,7 +149,9 @@ } /*}}}*/ - R2 Metric::mul(const R2 x)const {/*{{{*/ - return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y); - }/*}}}*/ + double Metric::Length(double Ax,double Ay) const{/*{{{*/ + /*Length of A in the current metric*/ + return sqrt(Ax*Ax*a11+2*Ax*Ay*a21+Ay*Ay*a22); + } + /*}}}*/ /*Intermediary*/ @@ -140,6 +168,6 @@ Ms1[level]=Ma; Ms2[level]=Mb; - double sa = Ma(AB); - double sb = Mb(AB); + double sa = Ma.Length(AB.x,AB.y); + double sb = Mb.Length(AB.x,AB.y); lMs1[level]=sa; lMs2[level]=sb; @@ -160,5 +188,5 @@ if( s > sstop && level < 30 && i < 500-level ) { Metric Mi(0.5,M1,0.5,M2); - double si = Mi(AB); + double si = Mi.Length(AB.x,AB.y); if( Abs((s1+s2)-(si+si)) > s1*0.001) { Index: /issm/trunk-jpl/src/c/bamg/Metric.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/Metric.h (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/Metric.h (revision 21623) @@ -30,8 +30,7 @@ Metric(double a); Metric(double a,double b,double c); - Metric( double a,const Metric& ma, double b,const Metric& mb); - Metric(const double a[3],const Metric& m0,const Metric& m1,const Metric& m2 ); + Metric(double a,const Metric ma,double b,const Metric mb); + Metric(const double a[3],const Metric m0,const Metric m1,const Metric m2 ); void Echo(); - R2 mul(const R2 x)const; double det() const; int IntersectWith(const Metric& M2); @@ -42,4 +41,5 @@ R2 Orthogonal(const R2 x){ return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y); } R2 Orthogonal(const I2 x){ return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y); } + double Length(double Ax,double Ay) const; //operators @@ -47,5 +47,5 @@ Metric operator/(double c) const {double c2=1/(c*c);return Metric(a11*c2,a21*c2,a22*c2);} operator D2xD2(){ return D2xD2(a11,a21,a21,a22);} - 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() + //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() FIXME: replace by Length! 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;}; @@ -57,5 +57,5 @@ //fields double lambda1,lambda2; - D2 v; + double vx,vy; //friends @@ -114,7 +114,7 @@ } inline Metric::Metric(const EigenMetric& M) { - double v00=M.v.x*M.v.x; - double v11=M.v.y*M.v.y; - double v01=M.v.x*M.v.y; + double v00=M.vx*M.vx; + double v11=M.vy*M.vy; + double v01=M.vx*M.vy; a11=v00*M.lambda1+v11*M.lambda2; a21=v01*(M.lambda1-M.lambda2); Index: /issm/trunk-jpl/src/c/bamg/Triangle.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/Triangle.cpp (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/Triangle.cpp (revision 21623) @@ -145,46 +145,4 @@ } /*}}}*/ - Triangle* Triangle::Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const{/*{{{*/ - // return the other triangle of the quad if a quad or 0 if not a quat - Triangle * t =0; - if (link) { - int a=-1; - if (AdjEdgeIndex[0] & 16 ) a=0; - if (AdjEdgeIndex[1] & 16 ) a=1; - if (AdjEdgeIndex[2] & 16 ) a=2; - if (a>=0) { - t = adj[a]; - // if (t-this<0) return 0; - v2 = vertices[VerticesOfTriangularEdge[a][0]]; - v0 = vertices[VerticesOfTriangularEdge[a][1]]; - v1 = vertices[OppositeEdge[a]]; - v3 = t->vertices[OppositeEdge[AdjEdgeIndex[a]&3]]; - } - } - return t; - } - /*}}}*/ - double Triangle::QualityQuad(int a,int option) const{/*{{{*/ - double q; - if (!link || AdjEdgeIndex[a] &4) - q= -1; - else { - Triangle * t = adj[a]; - if (t-this<0) q= -1;// because we do 2 times - else if (!t->link ) q= -1; - else if (AdjEdgeIndex[0] & 16 || AdjEdgeIndex[1] & 16 || AdjEdgeIndex[2] & 16 || t->AdjEdgeIndex[0] & 16 || t->AdjEdgeIndex[1] & 16 || t->AdjEdgeIndex[2] & 16 ) - q= -1; - else if(option){ - const BamgVertex & v2 = *vertices[VerticesOfTriangularEdge[a][0]]; - const BamgVertex & v0 = *vertices[VerticesOfTriangularEdge[a][1]]; - const BamgVertex & v1 = *vertices[OppositeEdge[a]]; - const BamgVertex & v3 = * t->vertices[OppositeEdge[AdjEdgeIndex[a]&3]]; - q = QuadQuality(v0,v1,v2,v3); // do the float part - } - else q= 1; - } - return q; - } - /*}}}*/ void Triangle::Renumbering(Triangle *tb,Triangle *te, long *renu){/*{{{*/ @@ -262,4 +220,7 @@ t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] &=55; // 23 + 32 AdjEdgeIndex[a] &=55 ; + }/*}}}*/ + Triangle* Triangle::TriangleAdj(int i) const {/*{{{*/ + return adj[i&3]; }/*}}}*/ int Triangle::swap(short a,int koption){/*{{{*/ @@ -305,11 +266,11 @@ // critere de Delaunay pure isotrope Icoor2 xb1 = sb->i.x - s1->i.x, - x21 = s2->i.x - s1->i.x, - yb1 = sb->i.y - s1->i.y, - y21 = s2->i.y - s1->i.y, - xba = sb->i.x - sa->i.x, - x2a = s2->i.x - sa->i.x, - yba = sb->i.y - sa->i.y, - y2a = s2->i.y - sa->i.y; + x21 = s2->i.x - s1->i.x, + yb1 = sb->i.y - s1->i.y, + y21 = s2->i.y - s1->i.y, + xba = sb->i.x - sa->i.x, + x2a = s2->i.x - sa->i.x, + yba = sb->i.y - sa->i.y, + y2a = s2->i.y - sa->i.y; double cosb12 = double(xb1*x21 + yb1*y21), @@ -342,5 +303,5 @@ if (Abs(d) > dd*1.e-3) { R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); - som = M(C - S2)/M(C - S1); + som = M.Length(C.x-S2.x,C.y-S2.y) / M.Length(C.x-S1.x,C.y-S1.y); } else {kopt=1;continue;} @@ -357,5 +318,5 @@ if(Abs(d) > dd*1.e-3) { R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); - som += M(C - S2)/M(C - S1); + som += M.Length(C.x-S2.x,C.y-S2.y) / M.Length(C.x-S1.x,C.y-S1.y); } else {kopt=1;continue;} @@ -376,7 +337,4 @@ } /*}}}*/ - Triangle* Triangle::TriangleAdj(int i) const {/*{{{*/ - return adj[i&3]; - }/*}}}*/ } Index: /issm/trunk-jpl/src/c/bamg/Triangle.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/Triangle.h (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/Triangle.h (revision 21623) @@ -46,9 +46,7 @@ int Hidden(int a)const; int GetAllflag(int a); - double QualityQuad(int a,int option=1) const; short NuEdgeTriangleAdj(int i) const; AdjacentTriangle Adj(int i) const; Triangle *TriangleAdj(int i) const; - Triangle *Quadrangle(BamgVertex *& v0,BamgVertex *& v1,BamgVertex *& v2,BamgVertex *& v3) const; void Renumbering(Triangle *tb,Triangle *te, long *renu); void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu); Index: /issm/trunk-jpl/src/c/bamg/bamgobjects.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/bamgobjects.h (revision 21622) +++ /issm/trunk-jpl/src/c/bamg/bamgobjects.h (revision 21623) @@ -11,6 +11,4 @@ #include "./BamgMesh.h" #include "./Metric.h" -#include "./DoubleAndInt.h" -#include "./Direction.h" #include "./BamgVertex.h" #include "./AdjacentTriangle.h" Index: /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp (revision 21622) +++ /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp (revision 21623) @@ -21,5 +21,4 @@ int i; int noerr=1; - double costheta=2; double hminaniso=1e-100; Mesh* Thr=NULL; @@ -168,8 +167,4 @@ if (Thr != &BTh) delete Thr; - //Make quadrangle if requested - if(costheta<=1.0) Th.MakeQuadrangles(costheta); - //if () Th.SplitElement(2); - //Split corners if requested if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices(); @@ -183,9 +178,6 @@ //display info if(verbosity>0) { - if (Th.nbt-Th.nbtout-Th.nbq*2){ - _printf_(" new number of triangles = " << (Th.nbt-Th.nbtout-Th.nbq*2) << "\n"); - } - if (Th.nbq ){ - _printf_(" new number of quads = " << Th.nbq << "\n"); + if (Th.nbt-Th.nbtout){ + _printf_(" new number of triangles = " << (Th.nbt-Th.nbtout) << "\n"); } }