Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 3300) +++ /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 3301) @@ -16,16 +16,12 @@ /*Bamg options*/ - double anisomax; int maxnbv; - double hmin,hmax; - double coef,maxsubdiv; + double coef; int verbosity; int NbSmooth; - double omega; /*intermediary*/ int noerr=1; int i,j,num; - int allquad=0; double costheta=2; double hminaniso=1e-100; @@ -35,12 +31,7 @@ /*Bamg options*/ - anisomax=bamgopts->anisomax; NbSmooth=bamgopts->NbSmooth; - omega=bamgopts->omega; coef=bamgopts->coef; maxnbv=bamgopts->maxnbv; - maxsubdiv=bamgopts->maxsubdiv; - hmin=bamgopts->hmin; - hmax=bamgopts->hmax; verbosity=bamgopts->verbose; @@ -55,22 +46,20 @@ /*Mesh generation {{{1*/ + //Step1: generate geometry Gh if (verbosity>0) printf("Construction of a mesh from a given geometry\n"); if (verbosity>1) printf(" Processing geometry...\n"); Geometry Gh(bamggeom_in,bamgopts); - if (verbosity>10){ - Gh.Echo(); - } //get hmin and hmax from geometry to generate the metric - hmin = Max(hmin,Gh.MinimalHmin()); - hmax = Min(hmax,Gh.MaximalHmax()); - - //build metric if not given in input + bamgopts->hmin = Max(bamgopts->hmin,Gh.MinimalHmin()); + bamgopts->hmax = Min(bamgopts->hmax,Gh.MaximalHmax()); + + //build metric using geometry if (verbosity>1) printf(" Generating Metric...\n"); for(i=0;ihmax); + Vp.Minh(bamgopts->hmin); Gh.vertices[i].m = Vp; } @@ -79,5 +68,9 @@ if (verbosity>1) printf(" Generating Mesh...\n"); Triangles Th(maxnbv,Gh,bamgopts); + + //Split corners if requested if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices(); + + //Renumbering Th.ReNumberingTheTriangleBySubDomain(); @@ -101,54 +94,95 @@ if (verbosity>1) printf(" Processing initial mesh and geometry...\n"); Triangles BTh(bamggeom_in,bamgmesh_in,bamgopts); - hmin=Max(hmin,BTh.MinimalHmin()); - hmax=Min(hmax,BTh.MaximalHmax()); //Make Quadtree from background mesh BTh.MakeQuadTree(); - //build metric if not given in input + //Bound hmin and hmax + bamgopts->hmin=Max(bamgopts->hmin,BTh.MinimalHmin()); + bamgopts->hmax=Min(bamgopts->hmax,BTh.MaximalHmax()); + + //Generate initial metric if (bamgopts->metric){ if (verbosity>1) printf(" Processing Metric...\n"); - BTh.ReadMetric(bamgopts,hmin,hmax,coef); + BTh.ReadMetric(bamgopts); } else { - - // init with hmax - Metric Mhmax(hmax); + if (verbosity>1) printf(" Generating initial Metric...\n"); + Metric Mhmax(bamgopts->hmax); for (int iv=0;ivhVertices){ - for (i=0;ihVertices[i])){ - BTh[i].m=Metric((float)bamgmesh_in->hVertices[i]); - } + } + + //use present fields to generate metric if present + if (bamgopts->field){ + if (verbosity>1) printf(" Merge metric with field provided...\n"); + BTh.AddMetric(bamgopts); + } + + // change using hVertices if provided + if (bamgmesh_in->hVertices){ + if (verbosity>1) printf(" Merging Metric with hVertices...\n"); + for (i=0;ihVertices[i])){ + BTh[i].m=Metric((float)bamgmesh_in->hVertices[i]); } - - } - - //use present fields to generate metric if present - if (bamgopts->field){ - if (verbosity>1) printf(" Generating Metric from solution field...\n"); - BTh.AddMetric(bamgopts); - } - } - - BTh.AddGeometryMetric(bamgopts); - if (verbosity>1) printf(" Smoothing Metric..."); - if(bamgopts->gradation) BTh.SmoothMetric(bamgopts->gradation); - if (verbosity>1) printf(" done\n"); - BTh.MaxSubDivision(maxsubdiv); - BTh.BoundAnisotropy(anisomax,hminaniso); - - //Build new mesh Thr + } + } + + // change using hminVertices if provided + if (bamgopts->hminVertices){ + if (verbosity>1) printf(" Merging Metric with hminVertices...\n"); + for (i=0;ihminVertices[i])){ + Metric M=BTh.vertices[i].m; + MatVVP2x2 Vp(M/coef); + Vp.Minh(bamgopts->hminVertices[i]); + BTh.vertices[i].m=Vp; + } + } + } + + // change using hmaxVertices if provided + if (bamgopts->hmaxVertices){ + if (verbosity>1) printf(" Merging Metric with hmaxVertices...\n"); + for (i=0;ihmaxVertices[i])){ + Metric M=BTh.vertices[i].m; + MatVVP2x2 Vp(M/coef); + Vp.Minh(bamgopts->hmaxVertices[i]); + BTh.vertices[i].m=Vp; + } + } + } + + //Add geometry metric if provided + if(bamgopts->geometricalmetric) BTh.AddGeometryMetric(bamgopts); + + //Smoothe metric + BTh.SmoothMetric(bamgopts->gradation); + + //Control element subdivision + BTh.MaxSubDivision(bamgopts->maxsubdiv); + + //Bound anisotropy + BTh.BoundAnisotropy(bamgopts->anisomax,hminaniso); + + //Build new mesh if (verbosity>1) printf(" Generating Mesh...\n"); Thr=&BTh,Thb=0; Triangles & Th( *(0 ? new Triangles(*Thr,&Thr->Gh,Thb,maxnbv) : new Triangles(maxnbv,BTh,bamgopts,bamgopts->KeepVertices))); if (Thr != &BTh) delete Thr; + + //Make quadrangle if requested if(costheta<=1.0) Th.MakeQuadrangles(costheta); - if (allquad) Th.SplitElement(allquad==2); + //if () Th.SplitElement(2); + + //Split corners if requested if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices(); + + //Renumber by subdomain Th.ReNumberingTheTriangleBySubDomain(); - if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega); + + //Smooth vertices + if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,bamgopts->omega); //display info Index: /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp (revision 3300) +++ /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp (revision 3301) @@ -25,71 +25,67 @@ double dca,dcb,dcta,dctb; double ddca,ddcb,ddcta,ddctb; - // double t1 = 1 -theta; - // double t1t1 = t1*t1; double tt = theta*theta; - if ( theta<0){ - throw ErrorException(__FUNCT__,exprintf("theta<0")); + + //check theta + assert(theta>=0 && theta<=1); + + if (TgA()){ + if (TgB()){ + // Tangent A and B provided: + // interpolation d'hermite + dcb = 6*theta*(1-theta); + ddcb = 6*(1-2*theta); + dca = -dcb; + ddca = -ddcb; + dcta = (3*theta - 4)*theta + 1; + ddcta=6*theta-4; + dctb = 3*tt - 2*theta; + ddctb = 6*theta-2; + } + else { + //Tangent A provided but tangent B not provided + // 1-t*t, t-t*t, t*t + double t = theta; + dcb = 2*t; + ddcb = 2; + dca = -dcb; + ddca = -2; + dcta = 1-dcb; + ddcta = -ddcb; + dctb=0; + ddctb=0; + } } - if ( theta>1){ - throw ErrorException(__FUNCT__,exprintf("theta>1")); + else{ + if (TgB()){ + //Tangent B provided but tangent A not provided + double t = 1-theta; + dca = -2*t; + ddca = 2; + dcb = -dca; + ddcb = -2; + dctb = 1+dca; + ddctb= ddca; + dcta =0; + ddcta =0; + } + else { + //Neither thangent A nor tangent B provided + // lagrange P1 + t=B-A; + return 0; + } } - if (TgA()) - if (TgB()) // interpolation d'hermite - { //cb = theta*theta*(3-2*theta); - dcb = 6*theta*(1-theta); - ddcb = 6*(1-2*theta); - //ca = 1-cb; - dca = -dcb; - ddca = -ddcb; - - // cta = (1-theta)*(1-theta)*theta; - dcta = (3*theta - 4)*theta + 1; - ddcta=6*theta-4; - - //ctb = (theta-1)*theta*theta ; - dctb = 3*tt - 2*theta; - ddctb = 6*theta-2; - } - else { // 1-t*t, t-t*t, t*t - double t = theta; - // cb = t*t; - dcb = 2*t; - ddcb = 2; - //ca = 1-cb; - dca = -dcb; - ddca = -2; - // cta= t-cb; - dcta = 1-dcb; - ddcta = -ddcb; - // ctb =0; - dctb=0; - ddctb=0; - } - else - if (TgB()){ - double t = 1-theta; - //ca = t*t; - dca = -2*t; - ddca = 2; - //cb = 1-ca; - dcb = -dca; - ddcb = -2; - //ctb= -t+ca; - dctb = 1+dca; - ddctb= ddca; - //cta=0; - dcta =0; - ddcta =0; - } - else {t=B-A;return 0;} // lagrange P1 - R2 d = A*dca + B*dcb + tg[0]* dcta + tg[1] * dctb; - + R2 d = A*dca + B*dcb + tg[0]* dcta + tg[1] * dctb; R2 dd = A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb; double d2=(d,d); double sd2 = sqrt(d2); t=d; - if(d2>1.0e-20) {t/=sd2;return Abs(Det(d,dd))/(d2*sd2);} + if(d2>1.0e-20){ + t/=sd2; + return Abs(Det(d,dd))/(d2*sd2); + } else return 0; - } + } /*}}}1*/ /*FUNCTION GeometricalEdge::F{{{1*/ Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 3300) +++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 3301) @@ -742,5 +742,4 @@ // generation of all the tangents for (i=0;ivertex2 R2 AB =edges[i].v[1]->r -edges[i].v[0]->r; @@ -748,10 +747,9 @@ double lAB=Norme2(AB); //initialize tangent - double ltg2[2]; - ltg2[0]=0;ltg2[1]=0; + double ltg2[2]={0.0}; //loop over the 2 vertices of the edge for (jj=0;jj<2;jj++) { - R2 tg =edges[i].tg[jj]; + R2 tg =edges[i].tg[jj]; double ltg=Norme2(tg); @@ -760,15 +758,22 @@ //if the current vertex of the edge is not a corner if(!edges[i].v[jj]->Corner()){ - tg = edges[i].v[1-jj]->r - edges[i].Adj[jj]->v[1-edges[i].DirAdj[jj]]->r; - ltg= Norme2(tg); - tg = tg *(lAB/ltg),ltg=lAB; + tg = edges[i].v[1-jj]->r - edges[i].Adj[jj]->v[1-edges[i].DirAdj[jj]]->r; //previous point and next point on curve + ltg= Norme2(tg); + tg = tg *(lAB/ltg); + ltg= lAB; } //else: a Corner with no tangent => nothing to do } else{ - tg = tg *(lAB/ltg),ltg=lAB; - } + //tangent has already been computed + tg = tg*(lAB/ltg),ltg=lAB; + } + ltg2[jj] = ltg; - if ((tg,AB)<0) tg = -tg; + + if ((tg,AB)<0){ + tg = -tg; + } + edges[i].tg[jj] = tg; } Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3300) +++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3301) @@ -822,11 +822,13 @@ /*}}}1*/ /*FUNCTION Triangles::ReadMetric{{{1*/ - void Triangles::ReadMetric(BamgOpts* bamgopts,const double hmin1=1.0e-30,const double hmax1=1.0e30,const double coef=1) { + void Triangles::ReadMetric(const BamgOpts* bamgopts) { + + /*Intermediary*/ int i,j; if(bamgopts->verbose>3) printf(" processing metric\n"); - - double hmin = Max(hmin1,MinimalHmin()); - double hmax = Min(hmax1,MaximalHmax()); + double hmin = Max(bamgopts->hmin,MinimalHmin()); + double hmax = Min(bamgopts->hmax,MaximalHmax()); + double coef = bamgopts->coef; //for now we only use j==3 @@ -904,5 +906,5 @@ Gh.ProjectOnCurve(edges[i],ss[j],V,GV); - GeometricalEdge * eg = GV; + GeometricalEdge* eg = GV; double s = GV; R2 tg; @@ -914,4 +916,5 @@ } double hn=Min(hmax,ht*anisomax); + if (ht<=0 || hn<=0){ throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0")); @@ -4418,10 +4421,10 @@ } } - if (quadtree) {delete quadtree; + if (quadtree){ + delete quadtree; quadtree = new QuadTree(this); } for ( it=0;it -#include -#include -#include #include "../Bamgx/objects/BamgObjects.h" Index: /issm/trunk/src/c/objects/BamgOpts.cpp =================================================================== --- /issm/trunk/src/c/objects/BamgOpts.cpp (revision 3300) +++ /issm/trunk/src/c/objects/BamgOpts.cpp (revision 3301) @@ -19,7 +19,10 @@ bamgopts->hmin=0; bamgopts->hmax=0; + bamgopts->hminVertices=NULL; + bamgopts->hmaxVertices=NULL; bamgopts->gradation=0; bamgopts->cutoff=0; bamgopts->splitcorners=0; + bamgopts->geometricalmetric=0; bamgopts->verbose=0; bamgopts->err=NULL; @@ -43,4 +46,5 @@ if (bamgopts->errg<0) throw ErrorException(__FUNCT__,exprintf("'errg' option should be >0")); if (bamgopts->nbjacobi<=0) throw ErrorException(__FUNCT__,exprintf("'nbjacobi' option should be >0")); + if (bamgopts->geometricalmetric!=0 && bamgopts->geometricalmetric!=1) throw ErrorException(__FUNCT__,exprintf("'geometricalmetric' supported options are 0 and 1")); if (bamgopts->NbSmooth<=0) throw ErrorException(__FUNCT__,exprintf("'NbSmooth' option should be >0")); if (bamgopts->maxnbv<3) throw ErrorException(__FUNCT__,exprintf("'maxnbv' option should be >3")); Index: /issm/trunk/src/c/objects/BamgOpts.h =================================================================== --- /issm/trunk/src/c/objects/BamgOpts.h (revision 3300) +++ /issm/trunk/src/c/objects/BamgOpts.h (revision 3301) @@ -22,10 +22,13 @@ double hmin; double hmax; + double* hminVertices; + double* hmaxVertices; double gradation; double cutoff; int splitcorners; + int geometricalmetric; int verbose; - double* err; - double errg; + double* err; + double errg; double coef; double* metric; Index: /issm/trunk/src/m/classes/public/bamg.m =================================================================== --- /issm/trunk/src/m/classes/public/bamg.m (revision 3300) +++ /issm/trunk/src/m/classes/public/bamg.m (revision 3301) @@ -4,8 +4,10 @@ % Available options (for more details see ISSM website http://issm.jpl.nasa.gov/): % -% - domain : followed by an ARGUS file that prescribes the domain outline (rifts and holes) -% - hmin : minimum edge length (default is 10^-100) -% - hmax : maximum esge length (default is 10^100) -% - hVertices : imposed edge length for each vertex (geometry or mesh) +% - domain: followed by an ARGUS file that prescribes the domain outline (rifts and holes) +% - hmin : minimum edge length (default is 10^-100) +% - hmax : maximum esge length (default is 10^100) +% - hVertices : imposed edge length for each vertex (geometry or mesh) +% - hminVertices: minimum edge length for each vertex (mesh) +% - hmaxVertices: maximum edge length for each vertex (mesh) % % - anisomax : maximum ration between the smallest and largest edges (default is 10^30) @@ -31,5 +33,6 @@ % - omega : relaxation parameter of the smoothing procedure (default is 1.8) % - power : power applied to the metric (default is 1) -% - splitcorner : split triangles whuch have 3 vertices on the outline (default is 1) +% - splitcorners : split triangles whuch have 3 vertices on the outline (default is 1) +% - geometricalmetric : Take the geometry into account to generate the metric (default is 0) % - verbose : level of verbosity (default is 1) % @@ -227,4 +230,6 @@ bamg_options.hmin=getfieldvalue(options,'hmin',10^-100); bamg_options.hmax=getfieldvalue(options,'hmax',10^100); +bamg_options.hminVertices=getfieldvalue(options,'hminVertices',[]); +bamg_options.hmaxVertices=getfieldvalue(options,'hmaxVertices',[]); bamg_options.KeepVertices=getfieldvalue(options,'KeepVertices',1); bamg_options.MaximalAngleOfCorner=getfieldvalue(options,'MaximalAngleOfCorner',10); @@ -238,17 +243,7 @@ bamg_options.power=getfieldvalue(options,'power',1); bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1); +bamg_options.geometricalmetric=getfieldvalue(options,'geometricalmetric',0); bamg_options.verbose=getfieldvalue(options,'verbose',1); %}}} - -%check bamg options -fieldsize=size(bamg_options.field); -if fieldsize(1), - if fieldsize(1)~=md.numberofgrids, - error(['bamg error message, ''field'' should have ' num2str(md.numberofgrids) ' lines (numberofgrids)']); - end - if (fieldsize(2)~=0 & length(bamg_options.err)~=fieldsize(2)), - error(['bamg error message, ''err'' should be of length ' num2str(fieldsize(2)) ' (As many as fields)']); - end -end %call Bamg Index: /issm/trunk/src/mex/Bamg/Bamg.cpp =================================================================== --- /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 3300) +++ /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 3301) @@ -82,5 +82,10 @@ FetchData(&bamgopts.verbose,mxGetField(BAMGOPTIONS,0,"verbose")); FetchData(&bamgopts.splitcorners,mxGetField(BAMGOPTIONS,0,"splitcorners")); + FetchData(&bamgopts.geometricalmetric,mxGetField(BAMGOPTIONS,0,"geometricalmetric")); FetchData(&bamgopts.MaximalAngleOfCorner,mxGetField(BAMGOPTIONS,0,"MaximalAngleOfCorner")); + FetchData(&bamgopts.hminVertices,&lines,&cols,mxGetField(BAMGOPTIONS,0,"hminVertices")); + if (bamgopts.hminVertices && (cols!=1 || lines!=bamgmesh_in.NumVertices)){throw ErrorException(__FUNCT__,exprintf("the size of 'hminVertices' should be [%i %i]",bamgmesh_in.NumVertices,1));} + FetchData(&bamgopts.hmaxVertices,&lines,&cols,mxGetField(BAMGOPTIONS,0,"hmaxVertices")); + if (bamgopts.hmaxVertices && (cols!=1 || lines!=bamgmesh_in.NumVertices)){throw ErrorException(__FUNCT__,exprintf("the size of 'hmaxVertices' should be [%i %i]",bamgmesh_in.NumVertices,1));} FetchData(&bamgopts.metric,&lines,&cols,mxGetField(BAMGOPTIONS,0,"metric")); if (bamgopts.metric && (cols!=3 || lines!=bamgmesh_in.NumVertices)){throw ErrorException(__FUNCT__,exprintf("the size of 'metric' should be [%i %i]",bamgmesh_in.NumVertices,3));}