Changeset 3301

Timestamp:

03/18/10 14:19:34 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added 2 VERY usefull options in bamg hminVertices and hmaxVertices

Location:

issm/trunk/src

Files:

• c/Bamgx/Bamgx.cpp (modified) (5 diffs)
• c/Bamgx/objects/GeometricalEdge.cpp (modified) (1 diff)
• c/Bamgx/objects/Geometry.cpp (modified) (3 diffs)
• c/Bamgx/objects/Triangles.cpp (modified) (4 diffs)
• c/Bamgx/objects/Triangles.h (modified) (1 diff)
• c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (modified) (1 diff)
• c/objects/BamgOpts.cpp (modified) (2 diffs)
• c/objects/BamgOpts.h (modified) (1 diff)
• m/classes/public/bamg.m (modified) (4 diffs)
• mex/Bamg/Bamg.cpp (modified) (1 diff)

issm/trunk/src/c/Bamgx/Bamgx.cpp

@@ -16 +16 @@
 
         /*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 +31 @@
 
         /*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 +46 @@
                 /*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;i<Gh.nbv;i++){
                         Metric M=Gh[i];
                         MatVVP2x2 Vp(M/coef);
-                        Vp.Maxh(hmax);
-                        Vp.Minh(hmin);
+                        Vp.Maxh(bamgopts->hmax);
+                        Vp.Minh(bamgopts->hmin);
                         Gh.vertices[i].m = Vp;
                 }
@@ -79 +68 @@
                 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 +94 @@
                 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;iv<BTh.nbv;iv++) BTh[iv].m = Mhmax;
-                        // change using hVertices if required
-                        if (bamgmesh_in->hVertices){
-                                for (i=0;i<BTh.nbv;i++){
-                                        if (!isnan(bamgmesh_in->hVertices[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;i<BTh.nbv;i++){
+                                if (!isnan(bamgmesh_in->hVertices[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;i<BTh.nbv;i++){
+                                if (!isnan(bamgopts->hminVertices[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;i<BTh.nbv;i++){
+                                if (!isnan(bamgopts->hmaxVertices[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

issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp

@@ -25 +25 @@
                 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*/

issm/trunk/src/c/Bamgx/objects/Geometry.cpp

@@ -742 +742 @@
                 // generation of  all the tangents
                 for (i=0;i<nbe;i++) {
-
                         //Get AB = vertex1->vertex2
                         R2    AB =edges[i].v[1]->r -edges[i].v[0]->r;        
@@ -748 +747 @@
                         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 +758 @@
                                         //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;
                         }

issm/trunk/src/c/Bamgx/objects/Triangles.cpp

@@ -822 +822 @@
         /*}}}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 +906 @@
                                 Gh.ProjectOnCurve(edges[i],ss[j],V,GV);
 
-                                GeometricalEdge * eg = GV;
+                                GeometricalEdge* eg = GV;
                                 double s = GV;
                                 R2 tg;
@@ -914 +916 @@
                                 }
                                 double hn=Min(hmax,ht*anisomax);
+
                                 if (ht<=0 || hn<=0){
                                         throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0"));
@@ -4418 +4421 @@
                                 }  
                         }
-                if (quadtree) {delete quadtree;
+                if (quadtree){
+                        delete quadtree;
                         quadtree = new QuadTree(this);
                 }
 
                 for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;
-                printf("end renumbervertex\n");
         }
         /*}}}1*/

issm/trunk/src/c/Bamgx/objects/Triangles.h

@@ -136 +136 @@
                         void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
                         void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
-                        void ReadMetric(BamgOpts* bamgopts,const double hmin,const double hmax,const double coef);
+                        void ReadMetric(const BamgOpts* bamgopts);
                         void WriteMetric(BamgOpts* bamgopts);
                         void AddMetric(BamgOpts* bamgopts);

issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp

@@ -11 +11 @@
 #include "../toolkits/toolkits.h"
 
-/*Bamg: */
-#include <unistd.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
 #include "../Bamgx/objects/BamgObjects.h"
 

issm/trunk/src/c/objects/BamgOpts.cpp

@@ -19 +19 @@
         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 +46 @@
         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"));

issm/trunk/src/c/objects/BamgOpts.h

@@ -22 +22 @@
         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;

issm/trunk/src/m/classes/public/bamg.m

@@ -4 +4 @@
 %   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 +33 @@
 %   - 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 +230 @@
 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 +243 @@
 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

issm/trunk/src/mex/Bamg/Bamg.cpp

@@ -82 +82 @@
         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));}