Index: /issm/trunk/src/mex/Bamg/Bamg.cpp =================================================================== --- /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 3210) +++ /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 3211) @@ -15,40 +15,4 @@ BamgGeom bamggeom; - /*Mesh inputs*/ - int NumVerticesMesh; - double* VerticesMesh=NULL; - int NumTrianglesMesh; - double* TrianglesMesh=NULL; - double* hVerticesMesh=NULL; - int NumCrackedEdgesMesh; - double* CrackedEdgesMesh=NULL; - - /*Geom inputs: */ - int NumVerticesGeom; - double* VerticesGeom=NULL; - int NumEdgesGeom; - double* EdgesGeom=NULL; - int NumCrackedEdgesGeom; - double* CrackedEdgesGeom=NULL; - double* hVerticesGeom=NULL; - double MaximalAngleOfCorner; - int NumSubDomainsGeom; - double* SubDomainsGeom=NULL; - - /*Options inputs*/ - int maxnbv,verbose,splitcorners; - double hmin,hmax,anisomax; - double coef,maxsubdiv; - double power; - int Hessiantype,Metrictype,NbSmooth; - int nbjacobi,KeepVertices,renumber; - double omega; - double gradation,errg; - double cutoff; - double* metric=NULL; - double* field=NULL; - double* err; - int numfields=0,numerr=0; - /*Boot module: */ MODULEBOOT(); @@ -58,134 +22,70 @@ /*create bamg geometry input*/ - FetchData(&NumVerticesGeom,mxGetField(BAMGGEOMETRY,0,"NumVertices")); - bamggeom.NumVertices=NumVerticesGeom; - FetchData(&VerticesGeom,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Vertices")); - bamggeom.Vertices=VerticesGeom; - FetchData(&NumEdgesGeom,mxGetField(BAMGGEOMETRY,0,"NumEdges")); - bamggeom.NumEdges=NumEdgesGeom; - FetchData(&EdgesGeom,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Edges")); - bamggeom.Edges=EdgesGeom; - FetchData(&hVerticesGeom,&lines,&cols,mxGetField(BAMGGEOMETRY,0,"hVertices")); - if (hVerticesGeom && (cols!=1 || lines!=NumVerticesGeom)){throw ErrorException(__FUNCT__,exprintf("the size of 'hVertices' should be [%i %i]",NumVerticesGeom,1));} - bamggeom.hVertices=hVerticesGeom; - bamggeom.MetricVertices=NULL; - bamggeom.h1h2VpVertices=NULL; - bamggeom.NumTangentAtEdges=0; - bamggeom.TangentAtEdges=NULL; - bamggeom.NumCorners=0; - bamggeom.Corners=NULL; - bamggeom.NumRequiredVertices=0; - bamggeom.RequiredVertices=NULL; - bamggeom.NumRequiredEdges=0; - bamggeom.RequiredEdges=NULL; - FetchData(&NumCrackedEdgesGeom,mxGetField(BAMGGEOMETRY,0,"NumCrackedEdges")); - bamggeom.NumCrackedEdges=NumCrackedEdgesGeom; - FetchData(&CrackedEdgesGeom,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"CrackedEdges")); - bamggeom.CrackedEdges=CrackedEdgesGeom; - FetchData(&NumSubDomainsGeom,mxGetField(BAMGGEOMETRY,0,"NumSubDomains")); - bamggeom.NumSubDomains=NumSubDomainsGeom; - FetchData(&SubDomainsGeom,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"SubDomains")); - bamggeom.SubDomains=SubDomainsGeom; + BamgGeomInit(&bamggeom); + FetchData(&bamggeom.NumVertices,mxGetField(BAMGGEOMETRY,0,"NumVertices")); + FetchData(&bamggeom.Vertices,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Vertices")); + FetchData(&bamggeom.NumEdges,mxGetField(BAMGGEOMETRY,0,"NumEdges")); + FetchData(&bamggeom.Edges,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Edges")); + FetchData(&bamggeom.hVertices,&lines,&cols,mxGetField(BAMGGEOMETRY,0,"hVertices")); + FetchData(&bamggeom.NumCrackedEdges,mxGetField(BAMGGEOMETRY,0,"NumCrackedEdges")); + FetchData(&bamggeom.CrackedEdges,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"CrackedEdges")); + FetchData(&bamggeom.NumSubDomains,mxGetField(BAMGGEOMETRY,0,"NumSubDomains")); + FetchData(&bamggeom.SubDomains,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"SubDomains")); + if (bamggeom.hVertices && (cols!=1 || lines!=bamggeom.NumVertices)){throw ErrorException(__FUNCT__,exprintf("the size of 'hVertices' should be [%i %i]",bamggeom.NumVertices,1));} /*create bamg mesh input*/ - FetchData(&NumVerticesMesh,mxGetField(BAMGMESH,0,"NumVertices")); - bamgmesh.NumVertices=NumVerticesMesh; - FetchData(&VerticesMesh,NULL,NULL,mxGetField(BAMGMESH,0,"Vertices")); - bamgmesh.Vertices=VerticesMesh; - FetchData(&NumTrianglesMesh,mxGetField(BAMGMESH,0,"NumTriangles")); - bamgmesh.NumTriangles=NumTrianglesMesh; - FetchData(&TrianglesMesh,NULL,NULL,mxGetField(BAMGMESH,0,"Triangles")); - bamgmesh.Triangles=TrianglesMesh; - FetchData(&hVerticesMesh,&lines,&cols,mxGetField(BAMGMESH,0,"hVertices")); - if (hVerticesMesh && (cols!=1 || lines!=NumVerticesMesh)){throw ErrorException(__FUNCT__,exprintf("the size of 'hVertices' should be [%i %i]",NumVerticesMesh,1));} - bamgmesh.hVertices=hVerticesMesh; - bamgmesh.NumQuadrilaterals=0; - bamgmesh.Quadrilaterals=NULL; - bamgmesh.NumVerticesOnGeometricVertex=0; - bamgmesh.VerticesOnGeometricVertex=NULL; - bamgmesh.NumVerticesOnGeometricEdge=0; - bamgmesh.VerticesOnGeometricEdge=NULL; - bamgmesh.NumEdgesOnGeometricEdge=0; - bamgmesh.EdgesOnGeometricEdge=NULL; - bamgmesh.NumEdges=0; - bamgmesh.Edges=NULL; - bamgmesh.NumSegments=0; - bamgmesh.Segments=NULL; - bamgmesh.SegmentsMarkers=NULL; - FetchData(&NumCrackedEdgesMesh,mxGetField(BAMGMESH,0,"NumCrackedEdges")); - bamgmesh.NumCrackedEdges=NumCrackedEdgesMesh; - FetchData(&CrackedEdgesMesh,NULL,NULL,mxGetField(BAMGMESH,0,"CrackedEdges")); - bamgmesh.CrackedEdges=CrackedEdgesMesh; - bamgmesh.NumSubDomains=0; - bamgmesh.SubDomains=NULL; - bamgmesh.NumSubDomainsFromGeom=0; - bamgmesh.SubDomainsFromGeom=NULL; + BamgMeshInit(&bamgmesh); + FetchData(&bamgmesh.NumVertices,mxGetField(BAMGMESH,0,"NumVertices")); + FetchData(&bamgmesh.Vertices,NULL,NULL,mxGetField(BAMGMESH,0,"Vertices")); + FetchData(&bamgmesh.NumTriangles,mxGetField(BAMGMESH,0,"NumTriangles")); + FetchData(&bamgmesh.Triangles,NULL,NULL,mxGetField(BAMGMESH,0,"Triangles")); + FetchData(&bamgmesh.hVertices,&lines,&cols,mxGetField(BAMGMESH,0,"hVertices")); + FetchData(&bamgmesh.NumCrackedEdges,mxGetField(BAMGMESH,0,"NumCrackedEdges")); + FetchData(&bamgmesh.CrackedEdges,NULL,NULL,mxGetField(BAMGMESH,0,"CrackedEdges")); + if (bamgmesh.hVertices && (cols!=1 || lines!=bamgmesh.NumVertices)){throw ErrorException(__FUNCT__,exprintf("the size of 'hVertices' should be [%i %i]",bamgmesh.NumVertices,1));} /*create bamg options input*/ - FetchData(&coef,mxGetField(BAMGOPTIONS,0,"coef")); - if (coef==0) throw ErrorException(__FUNCT__,exprintf("'coef' should be positive")); - bamgopts.coef=coef; - FetchData(&maxsubdiv,mxGetField(BAMGOPTIONS,0,"maxsubdiv")); - if (maxsubdiv<1) throw ErrorException(__FUNCT__,exprintf("'maxsubdiv' should be >=1")); - bamgopts.maxsubdiv=maxsubdiv; - FetchData(&Hessiantype,mxGetField(BAMGOPTIONS,0,"Hessiantype")); - if (Hessiantype!=0 && Hessiantype!=1) throw ErrorException(__FUNCT__,exprintf("'Hessiantype' supported options are 0 and 1")); - bamgopts.Hessiantype=Hessiantype; - FetchData(&Metrictype,mxGetField(BAMGOPTIONS,0,"Metrictype")); - if (Metrictype!=0 && Metrictype!=1 && Metrictype!=2) throw ErrorException(__FUNCT__,exprintf("'Metrictype' supported options are 0, 1 and 2")); - bamgopts.Metrictype=Metrictype; - FetchData(&KeepVertices,mxGetField(BAMGOPTIONS,0,"KeepVertices")); - if (KeepVertices!=0 && KeepVertices!=1) throw ErrorException(__FUNCT__,exprintf("'KeepVertices' supported options are 0 and 1")); - bamgopts.KeepVertices=KeepVertices; - FetchData(&renumber,mxGetField(BAMGOPTIONS,0,"renumber")); - if (renumber!=0 && renumber!=1 && renumber!=2) throw ErrorException(__FUNCT__,exprintf("'renumber' supported options are 0, 1 and 2")); - bamgopts.renumber=renumber; - FetchData(&power,mxGetField(BAMGOPTIONS,0,"power")); - bamgopts.power=power; - FetchData(&errg,mxGetField(BAMGOPTIONS,0,"errg")); - if (errg<0) throw ErrorException(__FUNCT__,exprintf("'errg' option should be >0")); - bamgopts.errg=errg; - FetchData(&nbjacobi,mxGetField(BAMGOPTIONS,0,"nbjacobi")); - if (nbjacobi<=0) throw ErrorException(__FUNCT__,exprintf("'nbjacobi' option should be >0")); - bamgopts.nbjacobi=nbjacobi; - FetchData(&NbSmooth,mxGetField(BAMGOPTIONS,0,"NbSmooth")); - if (NbSmooth<=0) throw ErrorException(__FUNCT__,exprintf("'NbSmooth' option should be >0")); - bamgopts.NbSmooth=NbSmooth; - FetchData(&omega,mxGetField(BAMGOPTIONS,0,"omega")); - bamgopts.omega=omega; - FetchData(&maxnbv,mxGetField(BAMGOPTIONS,0,"maxnbv")); - if (maxnbv<3) throw ErrorException(__FUNCT__,exprintf("'maxnbv' option should be >3")); - bamgopts.maxnbv=maxnbv; - FetchData(&hmin,mxGetField(BAMGOPTIONS,0,"hmin")); - if (hmin<=0) throw ErrorException(__FUNCT__,exprintf("'hmin' option should be >0")); - bamgopts.hmin=hmin; - FetchData(&hmax,mxGetField(BAMGOPTIONS,0,"hmax")); - if (hmax<=0 || hmax=1")); - bamgopts.anisomax=anisomax; - FetchData(&gradation,mxGetField(BAMGOPTIONS,0,"gradation")); - if (anisomax<1) throw ErrorException(__FUNCT__,exprintf("'anisomax' option should be >=1")); - bamgopts.gradation=gradation; - FetchData(&cutoff,mxGetField(BAMGOPTIONS,0,"cutoff")); - bamgopts.cutoff=cutoff; - FetchData(&verbose,mxGetField(BAMGOPTIONS,0,"verbose")); - bamgopts.verbose=verbose; - FetchData(&splitcorners,mxGetField(BAMGOPTIONS,0,"splitcorners")); - bamgopts.splitcorners=splitcorners; - FetchData(&MaximalAngleOfCorner,mxGetField(BAMGOPTIONS,0,"MaximalAngleOfCorner")); - bamgopts.MaximalAngleOfCorner=MaximalAngleOfCorner; - FetchData(&metric,&lines,&cols,mxGetField(BAMGOPTIONS,0,"metric")); - if (metric && (cols!=3 || lines!=NumVerticesMesh)){throw ErrorException(__FUNCT__,exprintf("the size of 'metric' should be [%i %i]",NumVerticesMesh,3));} - bamgopts.metric=metric; - FetchData(&field,&lines,&numfields,mxGetField(BAMGOPTIONS,0,"field")); - if (field && lines!=NumVerticesMesh){throw ErrorException(__FUNCT__,exprintf("the size of 'field' should be [%i %i]",NumVerticesMesh,numfields));} - bamgopts.numfields=numfields; - bamgopts.field=field; - FetchData(&err,NULL,&cols,mxGetField(BAMGOPTIONS,0,"err")); - if (numfields!=0 && cols!=numfields){throw ErrorException(__FUNCT__,exprintf("the size of 'err' should be the same as 'field'"));} - for (i=0;i0"));}; - bamgopts.err=err; + BamgOptsInit(&bamgopts); + FetchData(&bamgopts.coef,mxGetField(BAMGOPTIONS,0,"coef")); + if (bamgopts.coef==0) throw ErrorException(__FUNCT__,exprintf("'coef' should be positive")); + FetchData(&bamgopts.maxsubdiv,mxGetField(BAMGOPTIONS,0,"maxsubdiv")); + if (bamgopts.maxsubdiv<1) throw ErrorException(__FUNCT__,exprintf("'maxsubdiv' should be >=1")); + FetchData(&bamgopts.Hessiantype,mxGetField(BAMGOPTIONS,0,"Hessiantype")); + if (bamgopts.Hessiantype!=0 && bamgopts.Hessiantype!=1) throw ErrorException(__FUNCT__,exprintf("'Hessiantype' supported options are 0 and 1")); + FetchData(&bamgopts.Metrictype,mxGetField(BAMGOPTIONS,0,"Metrictype")); + if (bamgopts.Metrictype!=0 && bamgopts.Metrictype!=1 && bamgopts.Metrictype!=2) throw ErrorException(__FUNCT__,exprintf("'Metrictype' supported options are 0, 1 and 2")); + FetchData(&bamgopts.KeepVertices,mxGetField(BAMGOPTIONS,0,"KeepVertices")); + if (bamgopts.KeepVertices!=0 && bamgopts.KeepVertices!=1) throw ErrorException(__FUNCT__,exprintf("'KeepVertices' supported options are 0 and 1")); + FetchData(&bamgopts.renumber,mxGetField(BAMGOPTIONS,0,"renumber")); + if (bamgopts.renumber!=0 && bamgopts.renumber!=1 && bamgopts.renumber!=2) throw ErrorException(__FUNCT__,exprintf("'renumber' supported options are 0, 1 and 2")); + FetchData(&bamgopts.power,mxGetField(BAMGOPTIONS,0,"power")); + FetchData(&bamgopts.errg,mxGetField(BAMGOPTIONS,0,"errg")); + if (bamgopts.errg<0) throw ErrorException(__FUNCT__,exprintf("'errg' option should be >0")); + FetchData(&bamgopts.nbjacobi,mxGetField(BAMGOPTIONS,0,"nbjacobi")); + if (bamgopts.nbjacobi<=0) throw ErrorException(__FUNCT__,exprintf("'nbjacobi' option should be >0")); + FetchData(&bamgopts.NbSmooth,mxGetField(BAMGOPTIONS,0,"NbSmooth")); + if (bamgopts.NbSmooth<=0) throw ErrorException(__FUNCT__,exprintf("'NbSmooth' option should be >0")); + FetchData(&bamgopts.omega,mxGetField(BAMGOPTIONS,0,"omega")); + FetchData(&bamgopts.maxnbv,mxGetField(BAMGOPTIONS,0,"maxnbv")); + if (bamgopts.maxnbv<3) throw ErrorException(__FUNCT__,exprintf("'maxnbv' option should be >3")); + FetchData(&bamgopts.hmin,mxGetField(BAMGOPTIONS,0,"hmin")); + if (bamgopts.hmin<=0) throw ErrorException(__FUNCT__,exprintf("'hmin' option should be >0")); + FetchData(&bamgopts.hmax,mxGetField(BAMGOPTIONS,0,"hmax")); + if (bamgopts.hmax<=0 || bamgopts.hmax=1")); + FetchData(&bamgopts.gradation,mxGetField(BAMGOPTIONS,0,"gradation")); + if (bamgopts.gradation<1) throw ErrorException(__FUNCT__,exprintf("'gradation' option should be >=1")); + FetchData(&bamgopts.cutoff,mxGetField(BAMGOPTIONS,0,"cutoff")); + FetchData(&bamgopts.verbose,mxGetField(BAMGOPTIONS,0,"verbose")); + FetchData(&bamgopts.splitcorners,mxGetField(BAMGOPTIONS,0,"splitcorners")); + FetchData(&bamgopts.MaximalAngleOfCorner,mxGetField(BAMGOPTIONS,0,"MaximalAngleOfCorner")); + FetchData(&bamgopts.metric,&lines,&cols,mxGetField(BAMGOPTIONS,0,"metric")); + if (bamgopts.metric && (cols!=3 || lines!=bamgmesh.NumVertices)){throw ErrorException(__FUNCT__,exprintf("the size of 'metric' should be [%i %i]",bamgmesh.NumVertices,3));} + FetchData(&bamgopts.field,&lines,&bamgopts.numfields,mxGetField(BAMGOPTIONS,0,"field")); + if (bamgopts.field && lines!=bamgmesh.NumVertices){throw ErrorException(__FUNCT__,exprintf("the size of 'field' should be [%i %i]",bamgmesh.NumVertices,bamgopts.numfields));} + FetchData(&bamgopts.err,NULL,&cols,mxGetField(BAMGOPTIONS,0,"err")); + if (bamgopts.numfields!=0 && cols!=bamgopts.numfields){throw ErrorException(__FUNCT__,exprintf("the size of 'err' should be the same as 'field'"));} + for (i=0;i0"));}; /*!Generate internal degree of freedom numbers: */