Index: /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp (revision 24112) +++ /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp (revision 24113) @@ -28,5 +28,5 @@ this->hminVertices=NULL; this->hminVerticesSize[0]=this->hminVerticesSize[1]=0; this->hmaxVertices=NULL; this->hmaxVerticesSize[0]=this->hmaxVerticesSize[1]=0; - this->hVertices=NULL; this->hVerticesSize[0]=this->hVerticesSize[1]=0; + this->hVertices=NULL; this->hVerticesLength=0; this->metric=NULL; this->metricSize[0]=this->metricSize[1]=0; this->field=NULL; this->fieldSize[0]=this->fieldSize[1]=0; @@ -70,5 +70,5 @@ if (this->hminVertices && this->hminVerticesSize[1]!=1) _error_("'hminVertices' should be a column"); if (this->hmaxVertices && this->hmaxVerticesSize[1]!=1) _error_("'hmaxVertices' should be a column"); - if (this->hVertices && this->hVerticesSize[1]!=1) _error_("'hVertices' should be a column"); + if (this->hVertices && this->hVerticesLength<2) _error_("'hVertices' should be a vector"); if (this->metric && (this->metricSize[1]!=1 && this->metricSize[1]!=3)) _error_("'metric' should have either 1 (iso) or 3 (aniso) columns."); if (this->field){ Index: /issm/trunk-jpl/src/c/bamg/BamgOpts.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgOpts.h (revision 24112) +++ /issm/trunk-jpl/src/c/bamg/BamgOpts.h (revision 24113) @@ -1,5 +1,5 @@ /*!\file: BamgOpts.h * \brief place holder for optimization function arguments - */ + */ #ifndef _BAMGOPTS_H_ @@ -38,5 +38,5 @@ int hmaxVerticesSize[2]; double* hmaxVertices; - int hVerticesSize[2]; + int hVerticesLength; double* hVertices; int metricSize[2]; Index: /issm/trunk-jpl/src/c/bamg/Geometry.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/Geometry.cpp (revision 24112) +++ /issm/trunk-jpl/src/c/bamg/Geometry.cpp (revision 24113) @@ -57,5 +57,5 @@ nbcurves=0; - double Hmin = HUGE_VAL;// the infinie value + double Hmin = HUGE_VAL;// the infinie value int i,j,n,i0,i1,i2,i3; @@ -96,5 +96,5 @@ /*coefIcoor is the coefficient used for integer coordinates: * (x-pmin.x) - * Icoor x = (2^30 -1) ------------ + * Icoor x = (2^30 -1) ------------ * D * where D is the longest side of the domain (direction x or y) @@ -155,9 +155,9 @@ } - // definition the default of the given mesh size + // definition the default of the given mesh size for (i=0;i 0) + if (vertices[i].color > 0) vertices[i].m=Metric(verticeslength[i] /(double) vertices[i].color); - else + else vertices[i].m=Metric(Hmin); } @@ -170,5 +170,5 @@ //hVertices - if(bamgopts->hVertices && bamgopts->hVerticesSize[0]==nbv){ + if(bamgopts->hVertices && bamgopts->hVerticesLength==nbv){ if(verbose>5) _printf_(" processing hVertices\n"); for (i=0;i< nbv;i++){ @@ -211,5 +211,5 @@ if (bamggeom->CornersSize[1]!=1) _error_("Corners should have 1 column"); n=bamggeom->CornersSize[0]; - for (i=0;iCorners[i]-1; //for C indexing if (j>nbv-1 || j<0) _error_("Bad corner definition: should in [0 " << nbv << "]"); @@ -225,5 +225,5 @@ if (bamggeom->RequiredVerticesSize[1]!=1) _error_("RequiredVertices should have 1 column"); n=bamggeom->RequiredVerticesSize[0]; - for (i=0;iRequiredVertices[i]-1; //for C indexing if (j>nbv-1 || j<0) _error_("Bad RequiredVerticess definition: should in [0 " << nbv << "]"); @@ -237,8 +237,8 @@ if (bamggeom->RequiredEdgesSize[1]!=1) _error_("RequiredEdges should have 1 column"); n=bamggeom->RequiredEdgesSize[0]; - for (i=0;iRequiredEdges[i]-1; //for C indexing if (j>nbe-1 || j<0) _error_("Bad RequiredEdges definition: should in [0 " << nbe << "]"); - edges[j].SetRequired(); + edges[j].SetRequired(); } } @@ -457,10 +457,10 @@ the vertices by groups of 5: All the coordinates are transformed to ]0,1[^2 - then, the integer coordinates are computed using + then, the integer coordinates are computed using the transformation ]0,1[^2 -> [0 2^30-1[^2 for a quadtree of depth 30*/ - vertices[i].i=R2ToI2(vertices[i].r); + vertices[i].i=R2ToI2(vertices[i].r); /*find nearest vertex already present in the quadtree (NULL if empty)*/ - BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); + BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); /*if there is a vertex found that is to close to vertices[i] -> error*/ @@ -506,5 +506,5 @@ * head_F[j]=i; * } - * + * * Then, we can go through all the elements that have for image j: * for(i=head_F[j]; i!=-1; i=next_F[i]) @@ -512,14 +512,14 @@ * stop the loop when i=-1 (end of the chain) * iterate using i=next_F[i] (next element that have for image j) - * + * * 2. How to use this algorithm here? - * + * * Here F is a function that associates two vertices v0 and v1 for a given edge E * We want to build the reciprocal function: what are the edges that contains * a vertex v? * To do so, we use the same chaining algorithm but there is a difficulty - * coming from the fact that for F we have a couple of vertices and not one + * coming from the fact that for F we have a couple of vertices and not one * vertex. - * To overcome this difficulty, we use a global indexing exactly like in + * To overcome this difficulty, we use a global indexing exactly like in * C/C++ so that * a member of a 2-column-table can be described by one index p=i*2+j @@ -566,11 +566,11 @@ //sort head_v by order of increasing edges angle for (i=0;iv1 - float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1]; + float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1]; //exchange if the current edge angle is smaller than the previous one @@ -593,5 +593,5 @@ exch=1; *pn=*po; // next_p[n] = n + 1 - *po=*p; // + *po=*p; // *p=n; // next_p[n+1] = n po=pn; // po now point toward pn (invert next and current) @@ -634,6 +634,6 @@ else{ /*all vertices provided in geometry are corners (ord = number of edges holding i)*/ - vertices[i].SetCorner() ; - if(ord==2){ + vertices[i].SetCorner() ; + if(ord==2){ long n1 = head_v[i]; long n2 = next_p[n1]; @@ -652,5 +652,5 @@ /*close the list around the vertex to have a circular loop*/ long no=-1, ne = head_v[i]; - while (ne >=0) ne = next_p[no=ne]; + while (ne >=0) ne = next_p[no=ne]; if(no>=0) next_p[no] = head_v[i]; } @@ -662,5 +662,5 @@ for (j=0;j<2;j++){ - long n1 = next_p[k++]; + long n1 = next_p[k++]; long i1 = n1/2 ,j1=n1%2; @@ -695,5 +695,5 @@ ltg= lAB; } - //else: a Corner no tangent => nothing to do + //else: a Corner no tangent => nothing to do } else{ @@ -710,5 +710,5 @@ if (ltg2[0]!=0) edges[i].SetTgA(); if (ltg2[1]!=0) edges[i].SetTgB(); - } + } /* generation of all curves (from corner to corner)*/ @@ -726,12 +726,12 @@ for (i=0;iSetMark(); nb_marked_edges++; @@ -765,10 +765,10 @@ } } - } + } } _assert_(nb_marked_edges && nbe); //allocate if first step if(step==0) curves=new Curve[nbcurves]; - } + } /*clean up*/ @@ -812,5 +812,5 @@ //Get edge direction and swap v0 and v1 if necessary - R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; + R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; int OppositeSens = (V01,AB)<0; int direction0=0,direction1=1; @@ -828,5 +828,5 @@ directionge[bge]=1; - while(eg0!=(GeomEdge*)vg0 && (*eg0)(direction0)!=(GeomVertex*)vg0){ + while(eg0!=(GeomEdge*)vg0 && (*eg0)(direction0)!=(GeomVertex*)vg0){ if (bge<=0) { if(NbTry) { @@ -846,6 +846,6 @@ direction0 = 1-( directionge[bge] = tmpge->AdjVertexIndex[direction0]); } - while (eg1 != (GeomEdge*) vg1 && (*eg1)(direction1) != (GeomVertex*) vg1) { - if(tge>=mxe ) { + while (eg1 != (GeomEdge*) vg1 && (*eg1)(direction1) != (GeomVertex*) vg1) { + if(tge>=mxe ) { _printf_("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)\n"); NbTry++; @@ -871,5 +871,5 @@ double sg; - if (eg0 == eg1) { + if (eg0 == eg1) { double s0=vg0,s1=vg1; sg = s0*(1.0-s) + s*s1; @@ -886,5 +886,5 @@ lge[i]=ll += Norme2(AA-BB); AA=BB ;} - lge[tge]=ll+=Norme2(AA-V1); + lge[tge]=ll+=Norme2(AA-V1); // search the geometrical edge _assert_(s<=1.0); @@ -900,10 +900,10 @@ } on=ge[i]; - if (i==tge) + if (i==tge) s1=vg1; s =(ls-l0)/(l1-l0); - sg =s0*(1.0-s)+s*s1; - } + sg =s0*(1.0-s)+s*s1; + } _assert_(on); V.r= on->F(sg); @@ -915,5 +915,5 @@ /*coefIcoor is the coefficient used for integer coordinates: * (x-pmin.x) - * Icoor x = (2^30 -1) ------------ + * Icoor x = (2^30 -1) ------------ * D * where D is the longest side of the domain (direction x or y) @@ -927,3 +927,3 @@ for (int i=0;iverbose; - // no metric -> no smoothing - if (bamgopts->metric==NULL) nbsmooth=0; + // no metric -> no smoothing + if (bamgopts->metric==NULL) nbsmooth=0; /*If no mesh in input, generate one*/ @@ -84,8 +84,8 @@ /*Anisotropic mesh adaptation {{{*/ - // read background mesh + // read background mesh if (verbosity>0) _printf_("Anisotropic mesh adaptation\n"); if (verbosity>1) _printf_(" Processing initial mesh and geometry...\n"); - Mesh BTh(bamggeom_in,bamgmesh_in,bamgopts); + Mesh BTh(bamggeom_in,bamgmesh_in,bamgopts); //Make Quadtree from background mesh @@ -101,5 +101,5 @@ BTh.ReadMetric(bamgopts); } - else { + else { if (verbosity>1) _printf_(" Generating initial Metric...\n"); Metric Mhmax(bamgopts->hmax); @@ -114,5 +114,5 @@ // change using hVertices if provided - if(bamgopts->hVertices && bamgopts->hVerticesSize[0]==BTh.nbv){ + if(bamgopts->hVertices && bamgopts->hVerticesLength==BTh.nbv){ if (verbosity>1) _printf_(" Merging Metric with hVertices...\n"); for (i=0;i