Index: /issm/trunk/src/c/modules/BamgConvertMeshx/BamgConvertMeshx.cpp =================================================================== --- /issm/trunk/src/c/modules/BamgConvertMeshx/BamgConvertMeshx.cpp (revision 5176) +++ /issm/trunk/src/c/modules/BamgConvertMeshx/BamgConvertMeshx.cpp (revision 5177) @@ -20,5 +20,5 @@ /*Options*/ - BamgOpts bamgopts; + BamgOpts* bamgopts=NULL; // read mesh @@ -28,7 +28,10 @@ //write mesh and geometry if (verbose) printf("Write Geometry\n"); - Th.Gh.WriteGeometry(bamggeom,&bamgopts); + Th.Gh.WriteGeometry(bamggeom,bamgopts); if (verbose) printf("Write Mesh\n"); - Th.WriteMesh(bamgmesh,&bamgopts); + Th.WriteMesh(bamgmesh,bamgopts); + + //clean up + delete bamgopts; /*No error return*/ Index: /issm/trunk/src/c/objects/Bamg/BamgGeom.cpp =================================================================== --- /issm/trunk/src/c/objects/Bamg/BamgGeom.cpp (revision 5176) +++ /issm/trunk/src/c/objects/Bamg/BamgGeom.cpp (revision 5177) @@ -20,4 +20,27 @@ } /*}}}*/ +/*FUNCTION BamgGeom::BamgGeom(mxArray* matlab_struct){{{1*/ +#ifdef _SERIAL_ +BamgGeom::BamgGeom(mxArray* matlab_struct){ + + int lines,cols; + + FetchData(&this->Vertices, &this->VerticesSize[0], &this->VerticesSize[1], mxGetField(matlab_struct,0,"Vertices")); + FetchData(&this->Edges, &this->EdgesSize[0], &this->EdgesSize[1], mxGetField(matlab_struct,0,"Edges")); + FetchData(&this->hVertices,&lines,&cols,mxGetField(matlab_struct,0,"hVertices")); + this->MetricVertices=NULL; + this->TangentAtEdgesSize[0]=0, this->TangentAtEdgesSize[1]=0; this->TangentAtEdges=NULL; + FetchData(&this->Corners, &this->CornersSize[0], &this->CornersSize[1], mxGetField(matlab_struct,0,"Corners")); + FetchData(&this->RequiredVertices,&this->RequiredVerticesSize[0],&this->RequiredVerticesSize[1],mxGetField(matlab_struct,0,"RequiredVertices")); + FetchData(&this->RequiredEdges, &this->RequiredEdgesSize[0], &this->RequiredEdgesSize[1], mxGetField(matlab_struct,0,"RequiredEdges")); + FetchData(&this->CrackedEdges, &this->CrackedEdgesSize[0], &this->CrackedEdgesSize[1], mxGetField(matlab_struct,0,"CrackedEdges")); + FetchData(&this->SubDomains, &this->SubDomainsSize[0], &this->SubDomainsSize[1], mxGetField(matlab_struct,0,"SubDomains")); + + /*Some checks*/ + if (this->hVertices && (cols!=1 || lines!=this->VerticesSize[0])){ISSMERROR("the size of 'hVertices' should be [%i %i]",this->VerticesSize[0],1);} + +} +#endif +/*}}}*/ /*FUNCTION BamgGeom::~BamgGeom(){{{1*/ BamgGeom::~BamgGeom(){ @@ -38,23 +61,4 @@ /*Methods*/ -/*FUNCTION BamgGeom::GetMatlabStructureFields{{{1*/ -#ifdef _SERIAL_ -void BamgGeom::GetMatlabStructureFields(mxArray* matlab_struct){ - - int lines,cols; - - FetchData(&this->Vertices, &this->VerticesSize[0], &this->VerticesSize[1], mxGetField(matlab_struct,0,"Vertices")); - FetchData(&this->Edges, &this->EdgesSize[0], &this->EdgesSize[1], mxGetField(matlab_struct,0,"Edges")); - FetchData(&this->Corners, &this->CornersSize[0], &this->CornersSize[1], mxGetField(matlab_struct,0,"Corners")); - FetchData(&this->RequiredVertices,&this->RequiredVerticesSize[0],&this->RequiredVerticesSize[1],mxGetField(matlab_struct,0,"RequiredVertices")); - FetchData(&this->RequiredEdges, &this->RequiredEdgesSize[0], &this->RequiredEdgesSize[1], mxGetField(matlab_struct,0,"RequiredEdges")); - FetchData(&this->CrackedEdges, &this->CrackedEdgesSize[0], &this->CrackedEdgesSize[1], mxGetField(matlab_struct,0,"CrackedEdges")); - FetchData(&this->SubDomains, &this->SubDomainsSize[0], &this->SubDomainsSize[1], mxGetField(matlab_struct,0,"SubDomains")); - FetchData(&this->hVertices,&lines,&cols,mxGetField(matlab_struct,0,"hVertices")); - if (this->hVertices && (cols!=1 || lines!=this->VerticesSize[0])){ISSMERROR("the size of 'hVertices' should be [%i %i]",this->VerticesSize[0],1);} - -} -#endif -/*}}}*/ /*FUNCTION BamgGeom::SetMatlabStructureFields{{{1*/ #ifdef _SERIAL_ @@ -62,51 +66,35 @@ /*Intermediary*/ - int i,i1,i2; - mxArray* pfield=NULL; - mxArray* pfield2=NULL; + int i; mxArray* output=NULL; int numfields=7; const char* fnames[numfields]; - int fsize[numfields][2]; - double** fpointer[numfields]; mwSize ndim=2; mwSize dimensions[2]={1,1}; - /*Build fnames and fsize names and sizes of each field*/ - i=-1; - fnames[++i] = "Vertices"; fsize[i][0]=this->VerticesSize[0]; fsize[i][1]=this->VerticesSize[1]; fpointer[i]=&this->Vertices; - fnames[++i] = "Edges"; fsize[i][0]=this->EdgesSize[0]; fsize[i][1]=this->EdgesSize[1]; fpointer[i]=&this->Edges; - fnames[++i] = "TangentAtEdges"; fsize[i][0]=this->TangentAtEdgesSize[0]; fsize[i][1]=this->TangentAtEdgesSize[1]; fpointer[i]=&this->TangentAtEdges; - fnames[++i] = "RequiredVertices";fsize[i][0]=this->RequiredVerticesSize[0];fsize[i][1]=this->RequiredVerticesSize[1];fpointer[i]=&this->RequiredVertices; - fnames[++i] = "RequiredEdges"; fsize[i][0]=this->RequiredEdgesSize[0]; fsize[i][1]=this->RequiredEdgesSize[1]; fpointer[i]=&this->RequiredEdges; - fnames[++i] = "CrackedEdges"; fsize[i][0]=this->CrackedEdgesSize[0]; fsize[i][1]=this->CrackedEdgesSize[1]; fpointer[i]=&this->CrackedEdges; - fnames[++i] = "SubDomains"; fsize[i][0]=this->SubDomainsSize[0]; fsize[i][1]=this->SubDomainsSize[1]; fpointer[i]=&this->SubDomains; - ISSMASSERT(i==numfields-1); + /*Initialize field names*/ + i=0; + fnames[i++] = "Vertices"; + fnames[i++] = "Edges"; + fnames[i++] = "TangentAtEdges"; + fnames[i++] = "RequiredVertices"; + fnames[i++] = "RequiredEdges"; + fnames[i++] = "CrackedEdges"; + fnames[i++] = "SubDomains"; + ISSMASSERT(i==numfields); /*Initialize Matlab structure*/ output=mxCreateStructArray(ndim,dimensions,numfields,fnames); - /*Add every field tu structure*/ - for(i=0;iVerticesSize[0], this->VerticesSize[1], this->Vertices); + i++; SetMatlabStructureField(output,"Edges", this->EdgesSize[0], this->EdgesSize[1], this->Edges); + i++; SetMatlabStructureField(output,"TangentAtEdges", this->TangentAtEdgesSize[0], this->TangentAtEdgesSize[1], this->TangentAtEdges); + i++; SetMatlabStructureField(output,"RequiredVertices",this->RequiredVerticesSize[0],this->RequiredVerticesSize[1],this->RequiredVertices); + i++; SetMatlabStructureField(output,"RequiredEdges", this->RequiredEdgesSize[0], this->RequiredEdgesSize[1], this->RequiredEdges); + i++; SetMatlabStructureField(output,"CrackedEdges", this->CrackedEdgesSize[0], this->CrackedEdgesSize[1], this->CrackedEdges); + i++; SetMatlabStructureField(output,"SubDomains", this->SubDomainsSize[0], this->SubDomainsSize[1], this->SubDomains); + ISSMASSERT(i==numfields); /*Assign output*/ @@ -116,2 +104,32 @@ #endif /*}}}*/ +/*FUNCTION BamgGeom::SetMatlabStructureField{{{1*/ +#ifdef _SERIAL_ +void BamgGeom::SetMatlabStructureField(mxArray* matlab_struct,const char* fieldname,int fieldrows,int fieldcols,double* fieldpointer){ + + /*Intermediary*/ + int i1,i2; + mxArray* pfield=NULL; + mxArray* pfield2=NULL; + + /*Copy field*/ + double* fieldcopy=NULL; + if (fieldrows*fieldcols){ + fieldcopy=(double*)xmalloc(fieldrows*fieldcols*sizeof(double)); + for(i1=0;i1TrianglesSize[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; + this->SubDomainsFromGeomSize[0]=0, this->SubDomainsFromGeomSize[1]=0; this->SubDomainsFromGeom=NULL; + this->CrackedVerticesSize[0]=0, this->CrackedVerticesSize[1]=0; this->CrackedVertices=NULL; + this->CrackedEdgesSize[0]=0, this->CrackedEdgesSize[1]=0; this->CrackedEdges=NULL; + this->VerticesOnGeometricVertexSize[0]=0, this->VerticesOnGeometricVertexSize[1]=0;this->VerticesOnGeometricVertex=NULL; this->VerticesOnGeometricEdgeSize[0]=0, this->VerticesOnGeometricEdgeSize[1]=0; this->VerticesOnGeometricEdge=NULL; this->EdgesOnGeometricEdgeSize[0]=0, this->EdgesOnGeometricEdgeSize[1]=0; this->EdgesOnGeometricEdge=NULL; - this->SubDomainsSize[0]=0, this->SubDomainsSize[1]=0; this->SubDomains=NULL; - this->SubDomainsFromGeomSize[0]=0, this->SubDomainsFromGeomSize[1]=0; this->SubDomainsFromGeom=NULL; + this->hVertices=NULL; + this->IssmEdgesSize[0]=0, this->IssmEdgesSize[1]=0; this->IssmEdges=NULL; this->IssmSegmentsSize[0]=0, this->IssmSegmentsSize[1]=0; this->IssmSegments=NULL; + this->ElementConnectivitySize[0]=0, this->ElementConnectivitySize[1]=0; this->ElementConnectivity=NULL; this->NodalConnectivitySize[0]=0, this->NodalConnectivitySize[1]=0; this->NodalConnectivity=NULL; this->NodalElementConnectivitySize[0]=0, this->NodalElementConnectivitySize[1]=0; this->NodalElementConnectivity=NULL; - this->CrackedVerticesSize[0]=0, this->CrackedVerticesSize[1]=0; this->CrackedVertices=NULL; - this->CrackedEdgesSize[0]=0, this->CrackedEdgesSize[1]=0; this->CrackedEdges=NULL; + } +/*}}}*/ +/*FUNCTION BamgMesh::BamgMesh(mxArray* matlab_struct){{{1*/ +#ifdef _SERIAL_ +BamgMesh::BamgMesh(mxArray* matlab_struct){ + + int lines,cols; + + FetchData(&this->Vertices, &this->VerticesSize[0], &this->VerticesSize[1], mxGetField(matlab_struct,0,"Vertices")); + FetchData(&this->Edges, &this->EdgesSize[0], &this->EdgesSize[1], mxGetField(matlab_struct,0,"Edges")); + FetchData(&this->Triangles, &this->TrianglesSize[0], &this->TrianglesSize[1], mxGetField(matlab_struct,0,"Triangles")); + this->QuadrilateralsSize[0]=0, this->QuadrilateralsSize[1]=0; this->Quadrilaterals=NULL; + + this->SubDomainsSize[0]=0, this->SubDomainsSize[1]=0; this->SubDomains=NULL; + this->SubDomainsFromGeomSize[0]=0, this->SubDomainsFromGeomSize[1]=0; this->SubDomainsFromGeom=NULL; + this->CrackedVerticesSize[0]=0, this->CrackedVerticesSize[1]=0; this->CrackedVertices=NULL; + FetchData(&this->CrackedEdges, &this->CrackedEdgesSize[0], &this->CrackedEdgesSize[1], mxGetField(matlab_struct,0,"CrackedEdges")); + + FetchData(&this->VerticesOnGeometricEdge, &this->VerticesOnGeometricEdgeSize[0], &this->VerticesOnGeometricEdgeSize[1], mxGetField(matlab_struct,0,"VerticesOnGeometricEdge")); + FetchData(&this->VerticesOnGeometricVertex,&this->VerticesOnGeometricVertexSize[0],&this->VerticesOnGeometricVertexSize[1],mxGetField(matlab_struct,0,"VerticesOnGeometricVertex")); + FetchData(&this->EdgesOnGeometricEdge, &this->EdgesOnGeometricEdgeSize[0], &this->EdgesOnGeometricEdgeSize[1], mxGetField(matlab_struct,0,"EdgesOnGeometricEdge")); + + FetchData(&this->hVertices, &lines, &cols, mxGetField(matlab_struct,0,"hVertices")); + + this->IssmEdgesSize[0]=0, this->IssmEdgesSize[1]=0; this->IssmEdges=NULL; + FetchData(&this->IssmSegments, &this->IssmSegmentsSize[0], &this->IssmSegmentsSize[1], mxGetField(matlab_struct,0,"IssmSegments")); + + this->ElementConnectivitySize[0]=0, this->ElementConnectivitySize[1]=0; this->ElementConnectivity=NULL; + this->NodalConnectivitySize[0]=0, this->NodalConnectivitySize[1]=0; this->NodalConnectivity=NULL; + this->NodalElementConnectivitySize[0]=0, this->NodalElementConnectivitySize[1]=0; this->NodalElementConnectivity=NULL; + + /*Checks*/ + if (this->hVertices && (cols!=1 || lines!=this->VerticesSize[0])){ISSMERROR("the size of 'hVertices' should be [%i %i]",this->VerticesSize[0],1);} + +} +#endif /*}}}*/ /*FUNCTION BamgMesh::~BamgMesh(){{{1*/ @@ -34,17 +75,23 @@ xfree((void**)&this->Triangles); xfree((void**)&this->Quadrilaterals); + + xfree((void**)&this->SubDomains); + xfree((void**)&this->SubDomainsFromGeom); + xfree((void**)&this->CrackedVertices); + xfree((void**)&this->CrackedEdges); + xfree((void**)&this->VerticesOnGeometricVertex); xfree((void**)&this->VerticesOnGeometricEdge); xfree((void**)&this->EdgesOnGeometricEdge); - xfree((void**)&this->SubDomains); - xfree((void**)&this->SubDomainsFromGeom); + xfree((void**)&this->hVertices); + xfree((void**)&this->IssmEdges); xfree((void**)&this->IssmSegments); + xfree((void**)&this->ElementConnectivity); xfree((void**)&this->NodalConnectivity); xfree((void**)&this->NodalElementConnectivity); - xfree((void**)&this->CrackedVertices); - xfree((void**)&this->CrackedEdges); + } @@ -52,24 +99,4 @@ /*Methods*/ -/*FUNCTION BamgMesh::GetMatlabStructureFields{{{1*/ -#ifdef _SERIAL_ -void BamgMesh::GetMatlabStructureFields(mxArray* matlab_struct){ - - int lines,cols; - - FetchData(&this->Triangles, &this->TrianglesSize[0], &this->TrianglesSize[1], mxGetField(matlab_struct,0,"Triangles")); - FetchData(&this->Vertices, &this->VerticesSize[0], &this->VerticesSize[1], mxGetField(matlab_struct,0,"Vertices")); - FetchData(&this->Edges, &this->EdgesSize[0], &this->EdgesSize[1], mxGetField(matlab_struct,0,"Edges")); - FetchData(&this->IssmSegments, &this->IssmSegmentsSize[0], &this->IssmSegmentsSize[1], mxGetField(matlab_struct,0,"IssmSegments")); - FetchData(&this->CrackedEdges, &this->CrackedEdgesSize[0], &this->CrackedEdgesSize[1], mxGetField(matlab_struct,0,"CrackedEdges")); - FetchData(&this->EdgesOnGeometricEdge, &this->EdgesOnGeometricEdgeSize[0], &this->EdgesOnGeometricEdgeSize[1], mxGetField(matlab_struct,0,"EdgesOnGeometricEdge")); - FetchData(&this->VerticesOnGeometricEdge, &this->VerticesOnGeometricEdgeSize[0], &this->VerticesOnGeometricEdgeSize[1], mxGetField(matlab_struct,0,"VerticesOnGeometricEdge")); - FetchData(&this->VerticesOnGeometricVertex,&this->VerticesOnGeometricVertexSize[0],&this->VerticesOnGeometricVertexSize[1],mxGetField(matlab_struct,0,"VerticesOnGeometricVertex")); - FetchData(&this->hVertices, &lines, &cols, mxGetField(matlab_struct,0,"hVertices")); - if (this->hVertices && (cols!=1 || lines!=this->VerticesSize[0])){ISSMERROR("the size of 'hVertices' should be [%i %i]",this->VerticesSize[0],1);} - -} -#endif -/*}}}*/ /*FUNCTION BamgMesh::SetMatlabStructureFields{{{1*/ #ifdef _SERIAL_ @@ -77,60 +104,53 @@ /*Intermediary*/ - int i,i1,i2; - mxArray* pfield=NULL; - mxArray* pfield2=NULL; + int i; mxArray* output=NULL; int numfields=16; const char* fnames[numfields]; - int fsize[numfields][2]; - double** fpointer[numfields]; mwSize ndim=2; mwSize dimensions[2]={1,1}; - /*Build fnames and fsize names and sizes of each field*/ - i=-1; - fnames[++i] = "Triangles"; fsize[i][0]=this->TrianglesSize[0]; fsize[i][1]=this->TrianglesSize[1]; fpointer[i]=&this->Triangles; - fnames[++i] = "Vertices"; fsize[i][0]=this->VerticesSize[0]; fsize[i][1]=this->VerticesSize[1]; fpointer[i]=&this->Vertices; - fnames[++i] = "Edges"; fsize[i][0]=this->EdgesSize[0]; fsize[i][1]=this->EdgesSize[1]; fpointer[i]=&this->Edges; - fnames[++i] = "IssmSegments"; fsize[i][0]=this->IssmSegmentsSize[0]; fsize[i][1]=this->IssmSegmentsSize[1]; fpointer[i]=&this->IssmSegments; - fnames[++i] = "IssmEdges"; fsize[i][0]=this->IssmEdgesSize[0]; fsize[i][1]=this->IssmEdgesSize[1]; fpointer[i]=&this->IssmEdges; - fnames[++i] = "Quadrilaterals"; fsize[i][0]=this->QuadrilateralsSize[0]; fsize[i][1]=this->QuadrilateralsSize[1]; fpointer[i]=&this->Quadrilaterals; - fnames[++i] = "VerticesOnGeometricVertex";fsize[i][0]=this->VerticesOnGeometricVertexSize[0]; fsize[i][1]=this->VerticesOnGeometricVertexSize[1];fpointer[i]=&this->VerticesOnGeometricVertex; - fnames[++i] = "VerticesOnGeometricEdge"; fsize[i][0]=this->VerticesOnGeometricEdgeSize[0]; fsize[i][1]=this->VerticesOnGeometricEdgeSize[1]; fpointer[i]=&this->VerticesOnGeometricEdge; - fnames[++i] = "EdgesOnGeometricEdge"; fsize[i][0]=this->EdgesOnGeometricEdgeSize[0]; fsize[i][1]=this->EdgesOnGeometricEdgeSize[1]; fpointer[i]=&this->EdgesOnGeometricEdge; - fnames[++i] = "SubDomains"; fsize[i][0]=this->SubDomainsSize[0]; fsize[i][1]=this->SubDomainsSize[1]; fpointer[i]=&this->SubDomains; - fnames[++i] = "SubDomainsFromGeom"; fsize[i][0]=this->SubDomainsFromGeomSize[0]; fsize[i][1]=this->SubDomainsFromGeomSize[1]; fpointer[i]=&this->SubDomainsFromGeom; - fnames[++i] = "ElementConnectivity"; fsize[i][0]=this->ElementConnectivitySize[0]; fsize[i][1]=this->ElementConnectivitySize[1]; fpointer[i]=&this->ElementConnectivity; - fnames[++i] = "NodalConnectivity"; fsize[i][0]=this->NodalConnectivitySize[0]; fsize[i][1]=this->NodalConnectivitySize[1]; fpointer[i]=&this->NodalConnectivity; - fnames[++i] = "NodalElementConnectivity"; fsize[i][0]=this->NodalElementConnectivitySize[0]; fsize[i][1]=this->NodalElementConnectivitySize[1]; fpointer[i]=&this->NodalElementConnectivity; - fnames[++i] = "CrackedVertices"; fsize[i][0]=this->CrackedVerticesSize[0]; fsize[i][1]=this->CrackedVerticesSize[1]; fpointer[i]=&this->CrackedVertices; - fnames[++i] = "CrackedEdges"; fsize[i][0]=this->CrackedEdgesSize[0]; fsize[i][1]=this->CrackedEdgesSize[1]; fpointer[i]=&this->CrackedEdges; - ISSMASSERT(i==numfields-1); + /*Initialize field names*/ + i=0; + fnames[i++] = "Triangles"; + fnames[i++] = "Vertices"; + fnames[i++] = "Edges"; + fnames[i++] = "IssmSegments"; + fnames[i++] = "IssmEdges"; + fnames[i++] = "Quadrilaterals"; + fnames[i++] = "VerticesOnGeometricVertex"; + fnames[i++] = "VerticesOnGeometricEdge"; + fnames[i++] = "EdgesOnGeometricEdge"; + fnames[i++] = "SubDomains"; + fnames[i++] = "SubDomainsFromGeom"; + fnames[i++] = "ElementConnectivity"; + fnames[i++] = "NodalConnectivity"; + fnames[i++] = "NodalElementConnectivity"; + fnames[i++] = "CrackedVertices"; + fnames[i++] = "CrackedEdges"; + ISSMASSERT(i==numfields); /*Initialize Matlab structure*/ output=mxCreateStructArray(ndim,dimensions,numfields,fnames); - /*Add every field tu structure*/ - for(i=0;iTrianglesSize[0], this->TrianglesSize[1], this->Triangles); + i++; SetMatlabStructureField(output,"Vertices", this->VerticesSize[0], this->VerticesSize[1], this->Vertices); + i++; SetMatlabStructureField(output,"Edges", this->EdgesSize[0], this->EdgesSize[1], this->Edges); + i++; SetMatlabStructureField(output,"IssmSegments", this->IssmSegmentsSize[0], this->IssmSegmentsSize[1], this->IssmSegments); + i++; SetMatlabStructureField(output,"IssmEdges", this->IssmEdgesSize[0], this->IssmEdgesSize[1], this->IssmEdges); + i++; SetMatlabStructureField(output,"Quadrilaterals", this->QuadrilateralsSize[0], this->QuadrilateralsSize[1], this->Quadrilaterals); + i++; SetMatlabStructureField(output,"VerticesOnGeometricVertex",this->VerticesOnGeometricVertexSize[0],this->VerticesOnGeometricVertexSize[1], this->VerticesOnGeometricVertex); + i++; SetMatlabStructureField(output,"VerticesOnGeometricEdge", this->VerticesOnGeometricEdgeSize[0], this->VerticesOnGeometricEdgeSize[1], this->VerticesOnGeometricEdge); + i++; SetMatlabStructureField(output,"EdgesOnGeometricEdge", this->EdgesOnGeometricEdgeSize[0], this->EdgesOnGeometricEdgeSize[1], this->EdgesOnGeometricEdge); + i++; SetMatlabStructureField(output,"SubDomains", this->SubDomainsSize[0], this->SubDomainsSize[1], this->SubDomains); + i++; SetMatlabStructureField(output,"SubDomainsFromGeom", this->SubDomainsFromGeomSize[0], this->SubDomainsFromGeomSize[1], this->SubDomainsFromGeom); + i++; SetMatlabStructureField(output,"ElementConnectivity", this->ElementConnectivitySize[0], this->ElementConnectivitySize[1], this->ElementConnectivity); + i++; SetMatlabStructureField(output,"NodalConnectivity", this->NodalConnectivitySize[0], this->NodalConnectivitySize[1], this->NodalConnectivity); + i++; SetMatlabStructureField(output,"NodalElementConnectivity", this->NodalElementConnectivitySize[0], this->NodalElementConnectivitySize[1], this->NodalElementConnectivity); + i++; SetMatlabStructureField(output,"CrackedVertices", this->CrackedVerticesSize[0], this->CrackedVerticesSize[1], this->CrackedVertices); + i++; SetMatlabStructureField(output,"CrackedEdges", this->CrackedEdgesSize[0], this->CrackedEdgesSize[1], this->CrackedEdges); + ISSMASSERT(i==numfields); /*Assign output*/ @@ -140,2 +160,32 @@ #endif /*}}}*/ +/*FUNCTION BamgMesh::SetMatlabStructureField{{{1*/ +#ifdef _SERIAL_ +void BamgMesh::SetMatlabStructureField(mxArray* matlab_struct,const char* fieldname,int fieldrows,int fieldcols,double* fieldpointer){ + + /*Intermediary*/ + int i1,i2; + mxArray* pfield=NULL; + mxArray* pfield2=NULL; + + /*Copy field*/ + double* fieldcopy=NULL; + if (fieldrows*fieldcols){ + fieldcopy=(double*)xmalloc(fieldrows*fieldcols*sizeof(double)); + for(i1=0;i1iso=0; this->maxnbv=0; this->MaxCornerAngle=0; + this->Crack=0; this->Hessiantype=0; this->Metrictype=0; this->KeepVertices=0; - this->Crack=0; this->maxsubdiv=0; this->power=0; @@ -38,4 +38,44 @@ } +/*}}}*/ +/*FUNCTION BamgOpts::BamgOpts(mxArray* matlab_struct){{{1*/ +#ifdef _SERIAL_ +BamgOpts::BamgOpts(mxArray* matlab_struct){ + + int lines,cols; + + FetchData(&this->maxnbv,mxGetField(matlab_struct,0,"maxnbv")); + FetchData(&this->MaxCornerAngle,mxGetField(matlab_struct,0,"MaxCornerAngle")); + FetchData(&this->Crack,mxGetField(matlab_struct,0,"Crack")); + FetchData(&this->Hessiantype,mxGetField(matlab_struct,0,"Hessiantype")); + FetchData(&this->Metrictype,mxGetField(matlab_struct,0,"Metrictype")); + FetchData(&this->KeepVertices,mxGetField(matlab_struct,0,"KeepVertices")); + FetchData(&this->maxsubdiv,mxGetField(matlab_struct,0,"maxsubdiv")); + FetchData(&this->power,mxGetField(matlab_struct,0,"power")); + FetchData(&this->anisomax,mxGetField(matlab_struct,0,"anisomax")); + FetchData(&this->nbsmooth,mxGetField(matlab_struct,0,"nbsmooth")); + FetchData(&this->nbjacobi,mxGetField(matlab_struct,0,"nbjacobi")); + FetchData(&this->omega,mxGetField(matlab_struct,0,"omega")); + FetchData(&this->hmin,mxGetField(matlab_struct,0,"hmin")); + FetchData(&this->hmax,mxGetField(matlab_struct,0,"hmax")); + FetchData(&this->hminVertices,&lines,&cols,mxGetField(matlab_struct,0,"hminVertices")); + FetchData(&this->hmaxVertices,&lines,&cols,mxGetField(matlab_struct,0,"hmaxVertices")); + FetchData(&this->gradation,mxGetField(matlab_struct,0,"gradation")); + FetchData(&this->cutoff,mxGetField(matlab_struct,0,"cutoff")); + FetchData(&this->splitcorners,mxGetField(matlab_struct,0,"splitcorners")); + FetchData(&this->geometricalmetric,mxGetField(matlab_struct,0,"geometricalmetric")); + FetchData(&this->verbose,mxGetField(matlab_struct,0,"verbose")); + FetchData(&this->field,&lines,&this->numfields,mxGetField(matlab_struct,0,"field")); + FetchData(&this->err,NULL,&cols,mxGetField(matlab_struct,0,"err")); + if (this->numfields!=0 && cols!=this->numfields){ISSMERROR("the size of 'err' should be the same as 'field'");} + FetchData(&this->errg,mxGetField(matlab_struct,0,"errg")); + FetchData(&this->coeff,mxGetField(matlab_struct,0,"coeff")); + FetchData(&this->metric,&lines,&cols,mxGetField(matlab_struct,0,"metric")); + + /*Additional checks*/ + this->Check(); + +} +#endif /*}}}*/ /*FUNCTION BamgOpts::~BamgOpts() {{{1*/ Index: /issm/trunk/src/c/objects/Bamg/BamgOpts.h =================================================================== --- /issm/trunk/src/c/objects/Bamg/BamgOpts.h (revision 5176) +++ /issm/trunk/src/c/objects/Bamg/BamgOpts.h (revision 5177) @@ -6,9 +6,12 @@ #define _BAMGOPTS_H_ +#ifdef _SERIAL_ +#include +#endif + class BamgOpts{ public: - int iso; int maxnbv; double MaxCornerAngle; @@ -40,4 +43,7 @@ BamgOpts(); + #ifdef _SERIAL_ + BamgOpts(mxArray* matlab_struct); + #endif ~BamgOpts(); Index: /issm/trunk/src/c/objects/Bamg/Geometry.cpp =================================================================== --- /issm/trunk/src/c/objects/Bamg/Geometry.cpp (revision 5176) +++ /issm/trunk/src/c/objects/Bamg/Geometry.cpp (revision 5177) @@ -111,4 +111,6 @@ * where D is the longest side of the domain (direction x or y) * so that (x-pmin.x)/D is in ]0 1[ + * + * coefIcoor = (2^30 -1)/D */ coefIcoor=(MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y)); @@ -398,19 +400,17 @@ /*Methods*/ - /*FUNCTION Geometry::AfterRead(){{{1*/ + /*FUNCTION Geometry::AfterRead{{{1*/ void Geometry::AfterRead(){ /*Original code from Frederic Hecht (BAMG v1.01, MeshGeom.cpp/AfterRead)*/ - long int verbose=0; - - long i,j,k; - int jj; - long* head_v=new long[nbv]; - long* next_p=new long[2*nbe]; - float* eangle=new float[nbe]; - double eps=1e-20; - QuadTree quadtree; // build quadtree to find duplicates - BamgVertex* v0=vertices; - GeometricalVertex* v0g=(GeometricalVertex*) (void*)v0; + long i,j,k; + int jj; + long *head_v = new long[nbv]; + long *next_p = new long[2*nbe]; + float *eangle = new float[nbe]; + double eps = 1e-20; + QuadTree quadtree; // build quadtree to find duplicates + BamgVertex *v0 = vertices; + GeometricalVertex *v0g = (GeometricalVertex*) (void*)v0; k=0; @@ -424,24 +424,17 @@ All the coordinates are transformed to ]0,1[^2 then, the integer coordinates are computed using - the transformation ]0,1[^2 -> [0 2^(30-1)[^2 for a quadtree of depth 30*/ + the transformation ]0,1[^2 -> [0 2^30-1[^2 for a quadtree of depth 30*/ vertices[i].i=toI2(vertices[i].r); - //find nearest vertex already present in the quadtree (NULL if empty) + /*find nearest vertex already present in the quadtree (NULL if empty)*/ 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 - if( v && Norme1(v->r - vertices[i]) < eps ){ - printf("WARNING: two points of the geometry are very closed to each other\n"); - } - - //Add vertices[i] to the quadtree + /*if there is a vertex found that is to close to vertices[i] -> error*/ + if( v && Norme1(v->r - vertices[i].r) < eps ){ + ISSMERROR("two points of the geometry are very closed to each other"); + } + + /*Add vertices[i] to the quadtree*/ quadtree.Add(vertices[i]); - } - - //if k>0, there are some duplicate vertices -> error - if (k) { - printf("number of distinct vertices= %i, over %i\n",nbv - k,nbv); - printf("List of duplicate vertices:\n"); - ISSMERROR("See above"); } @@ -767,11 +760,11 @@ nbv=0; nbe=0; - quadtree=0; - curves=0; - edges=0; - vertices=0; + quadtree=NULL; + curves=NULL; + edges=NULL; + vertices=NULL; nbsubdomains=0; nbcurves=0; - subdomains=0; + subdomains=NULL; MaxCornerAngle = 10*Pi/180; //default is 10 degres } @@ -941,6 +934,14 @@ /*FUNCTION Geometry::toI2{{{1*/ I2 Geometry::toI2(const R2 & P) const { - return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) - ,(Icoor1) (coefIcoor*(P.y-pmin.y)) ); + /*coefIcoor is the coefficient used for integer coordinates: + * (x-pmin.x) + * Icoor x = (2^30 -1) ------------ + * D + * where D is the longest side of the domain (direction x or y) + * so that (x-pmin.x)/D is in ]0 1[ + * + * coefIcoor = (2^30 -1)/D + */ + return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) ,(Icoor1) (coefIcoor*(P.y-pmin.y)) ); }/*}}}*/ /*FUNCTION Geometry::UnMarkEdges{{{1*/ Index: /issm/trunk/src/c/objects/Bamg/QuadTree.cpp =================================================================== --- /issm/trunk/src/c/objects/Bamg/QuadTree.cpp (revision 5176) +++ /issm/trunk/src/c/objects/Bamg/QuadTree.cpp (revision 5177) @@ -72,4 +72,23 @@ /*Constructors/Destructors*/ + /*FUNCTION QuadTree::QuadTree(){{{1*/ + QuadTree::QuadTree() : + /*Original code from Frederic Hecht (BAMG v1.01, QuadTree.cpp/QuadTree)*/ + + lenStorageQuadTreeBox(100), // by default 100 vertices by box + th(NULL), // initial mesh = NULL + NbQuadTreeBox(0), // initial number of quadtree boxes = 0 + NbVertices(0), // initial number of vertices = 0 + NbQuadTreeBoxSearch(0), // initial ?? = 0 + NbVerticesSearch(0){ // initial ?? = 0 + + //create lenStorageQuadTreeBox (100) StorageQuadTreeBox elements + sb =new StorageQuadTreeBox(lenStorageQuadTreeBox); + + //root=QuadTreeBox* pointer toward next quadtree box + root=NewQuadTreeBox(); + + } + /*}}}1*/ /*FUNCTION QuadTree::QuadTree(Mesh * t,long nbv){{{1*/ QuadTree::QuadTree(Mesh * t,long nbv) : @@ -84,5 +103,5 @@ { if (nbv == -1) nbv = t->nbv; - sb =new StorageQuadTreeBox(lenStorageQuadTreeBox); + sb = new StorageQuadTreeBox(lenStorageQuadTreeBox); root=NewQuadTreeBox(); if ( MaxISize <= MaxICoor){ @@ -93,20 +112,4 @@ } /*}}}1*/ - /*FUNCTION QuadTree::QuadTree(){{{1*/ - QuadTree::QuadTree() : - /*Original code from Frederic Hecht (BAMG v1.01, QuadTree.cpp/QuadTree)*/ - - lenStorageQuadTreeBox(100), // by default 100 vertices by box - th(0), // initial mesh = NULL - NbQuadTreeBox(0), // initial number of quadtree boxes = 0 - NbVertices(0), // initial number of vertices = 0 - NbQuadTreeBoxSearch(0), // initial ?? = 0 - NbVerticesSearch(0){ // initial ?? = 0 - //create lenStorageQuadTreeBox (100) StorageQuadTreeBox elements - sb =new StorageQuadTreeBox(lenStorageQuadTreeBox); - //root=QuadTreeBox* pointer roward ?? - root=NewQuadTreeBox(); - } - /*}}}1*/ /*FUNCTION QuadTree::~QuadTree(){{{1*/ QuadTree::~QuadTree() { @@ -115,4 +118,25 @@ delete sb; root=0; + } + /*}}}1*/ + /*FUNCTION QuadTree::StorageQuadTreeBox::StorageQuadTreeBox{{{1*/ + QuadTree::StorageQuadTreeBox::StorageQuadTreeBox(long ll,StorageQuadTreeBox *nn) { + /*Original code from Frederic Hecht (BAMG v1.01, QuadTree.cpp/StorageQuadTreeBox)*/ + + /*Initilalize variables*/ + len = ll; // number of quadtree boxes + n = nn; // next StorageQuadTreeBox pointer + b = new QuadTreeBox[ll]; // quadtree boxes + ISSMASSERT(b); // check allocation + + /*Initialize all quadtree boxes (empty)*/ + for (int i = 0; i (BAMG v1.01, QuadTree.cpp/Add)*/ - QuadTreeBox** pb; - QuadTreeBox* b; + QuadTreeBox** pb=NULL; + QuadTreeBox* b=NULL; + + /*Get integer coodinate of current point w*/ register long i=w.i.x, j=w.i.y; + + /*Initialize level*/ register long level=MaxISize; - //Get inital box (the largest) + /*Get inital box (the largest)*/ pb = &root; @@ -460,20 +488,4 @@ } /*}}}1*/ - /*FUNCTION QuadTree::StorageQuadTreeBox::StorageQuadTreeBox{{{1*/ - QuadTree::StorageQuadTreeBox::StorageQuadTreeBox(long ll,StorageQuadTreeBox *nn) { - /*Original code from Frederic Hecht (BAMG v1.01, QuadTree.cpp/StorageQuadTreeBox)*/ - - len = ll; - n = nn; - b = new QuadTreeBox[ll]; - for (int i = 0; i BamgVertex else => QuadTreeBox; + + int n; // number of current vertices in the box + union{ QuadTreeBox* b[4]; - BamgVertex* v[4]; + BamgVertex* v[4]; }; }; + class StorageQuadTreeBox{ + public: + + /*Fields*/ QuadTreeBox *b,*bc,*be; long len; StorageQuadTreeBox* n; // next StorageQuadTreeBox + + /*Methods*/ StorageQuadTreeBox(long ,StorageQuadTreeBox* =NULL); ~StorageQuadTreeBox() { @@ -39,4 +51,6 @@ long SizeOf() const {return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0);} }; + + /*QuadTree private Fields*/ StorageQuadTreeBox* sb; long lenStorageQuadTreeBox; @@ -44,14 +58,13 @@ public: - //fields + /*QuadTree public Fields*/ QuadTreeBox* root; - Mesh* th; + Mesh* th; + long NbQuadTreeBox,NbVertices; + long NbQuadTreeBoxSearch,NbVerticesSearch; - //functions - ~QuadTree(); QuadTree(Mesh *t,long nbv=-1); QuadTree(); - long NbQuadTreeBox,NbVertices; - long NbQuadTreeBoxSearch,NbVerticesSearch; + ~QuadTree(); BamgVertex* NearestVertex(Icoor1 i,Icoor1 j); BamgVertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j); @@ -63,8 +76,18 @@ * a private class and is declared before QuadTree::*/ QuadTreeBox* NewQuadTreeBox(void){ - if(! (sb->bcbe)) sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb); - if (!sb || (sb->bc->n != 0)){ISSMERROR("!sb || (sb->bc->n != 0)");} + + /*if bc==be or bc>be (we have reach the end of the StorageQuadTreeBox) + * create a new StorageQuadTreeBox)*/ + if(!(sb->bcbe)){ + sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb); + } + ISSMASSERT(sb && sb->bc->n==0); + + /*Increase counter*/ NbQuadTreeBox++; + + /*bc now points toward next quadtree box*/ return sb->bc++; + } }; Index: /issm/trunk/src/mex/Bamg/Bamg.cpp =================================================================== --- /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 5176) +++ /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 5177) @@ -11,6 +11,4 @@ /*diverse: */ - int i; - int lines,cols; BamgOpts *bamgopts=NULL; BamgMesh *bamgmesh_in=NULL; @@ -26,47 +24,11 @@ /*Initialize variables*/ - bamgopts=new BamgOpts; - bamggeom_in=new BamgGeom; - bamgmesh_in=new BamgMesh; - bamggeom_out=new BamgGeom; - bamgmesh_out=new BamgMesh; + bamgopts = new BamgOpts(BAMGOPTIONS); + bamggeom_in= new BamgGeom(BAMGGEOMETRY); + bamgmesh_in= new BamgMesh(BAMGMESH); - /*Build objects from matlab structures*/ - bamggeom_in->GetMatlabStructureFields(BAMGGEOMETRY); - bamgmesh_in->GetMatlabStructureFields(BAMGMESH); - - /*create bamg options input*/ - FetchData(&bamgopts->coeff,mxGetField(BAMGOPTIONS,0,"coeff")); - FetchData(&bamgopts->maxsubdiv,mxGetField(BAMGOPTIONS,0,"maxsubdiv")); - FetchData(&bamgopts->Crack,mxGetField(BAMGOPTIONS,0,"Crack")); - FetchData(&bamgopts->Hessiantype,mxGetField(BAMGOPTIONS,0,"Hessiantype")); - FetchData(&bamgopts->Metrictype,mxGetField(BAMGOPTIONS,0,"Metrictype")); - FetchData(&bamgopts->KeepVertices,mxGetField(BAMGOPTIONS,0,"KeepVertices")); - FetchData(&bamgopts->power,mxGetField(BAMGOPTIONS,0,"power")); - FetchData(&bamgopts->errg,mxGetField(BAMGOPTIONS,0,"errg")); - FetchData(&bamgopts->nbjacobi,mxGetField(BAMGOPTIONS,0,"nbjacobi")); - FetchData(&bamgopts->nbsmooth,mxGetField(BAMGOPTIONS,0,"nbsmooth")); - FetchData(&bamgopts->omega,mxGetField(BAMGOPTIONS,0,"omega")); - FetchData(&bamgopts->maxnbv,mxGetField(BAMGOPTIONS,0,"maxnbv")); - FetchData(&bamgopts->hmin,mxGetField(BAMGOPTIONS,0,"hmin")); - FetchData(&bamgopts->hmax,mxGetField(BAMGOPTIONS,0,"hmax")); - FetchData(&bamgopts->anisomax,mxGetField(BAMGOPTIONS,0,"anisomax")); - FetchData(&bamgopts->gradation,mxGetField(BAMGOPTIONS,0,"gradation")); - FetchData(&bamgopts->cutoff,mxGetField(BAMGOPTIONS,0,"cutoff")); - FetchData(&bamgopts->verbose,mxGetField(BAMGOPTIONS,0,"verbose")); - FetchData(&bamgopts->splitcorners,mxGetField(BAMGOPTIONS,0,"splitcorners")); - FetchData(&bamgopts->geometricalmetric,mxGetField(BAMGOPTIONS,0,"geometricalmetric")); - FetchData(&bamgopts->MaxCornerAngle,mxGetField(BAMGOPTIONS,0,"MaxCornerAngle")); - FetchData(&bamgopts->hminVertices,&lines,&cols,mxGetField(BAMGOPTIONS,0,"hminVertices")); - if (bamgopts->hminVertices && (cols!=1 || lines!=bamgmesh_in->VerticesSize[0])){ISSMERROR("the size of 'hminVertices' should be [%i %i]",bamgmesh_in->VerticesSize[0],1);} - FetchData(&bamgopts->hmaxVertices,&lines,&cols,mxGetField(BAMGOPTIONS,0,"hmaxVertices")); - if (bamgopts->hmaxVertices && (cols!=1 || lines!=bamgmesh_in->VerticesSize[0])){ISSMERROR("the size of 'hmaxVertices' should be [%i %i]",bamgmesh_in->VerticesSize[0],1);} - FetchData(&bamgopts->metric,&lines,&cols,mxGetField(BAMGOPTIONS,0,"metric")); - if (bamgopts->metric && (cols!=3 || lines!=bamgmesh_in->VerticesSize[0])){ISSMERROR("the size of 'metric' should be [%i %i]",bamgmesh_in->VerticesSize[0],3);} - FetchData(&bamgopts->field,&lines,&bamgopts->numfields,mxGetField(BAMGOPTIONS,0,"field")); - if (bamgopts->field && lines!=bamgmesh_in->VerticesSize[0]){ISSMERROR("the size of 'field' should be [%i %i]",bamgmesh_in->VerticesSize[0],bamgopts->numfields);} - FetchData(&bamgopts->err,NULL,&cols,mxGetField(BAMGOPTIONS,0,"err")); - if (bamgopts->numfields!=0 && cols!=bamgopts->numfields){ISSMERROR("the size of 'err' should be the same as 'field'");} - bamgopts->Check(); + /*Initialize outputs*/ + bamggeom_out=new BamgGeom(); + bamgmesh_out=new BamgMesh(); /*!Generate internal degree of freedom numbers: */ Index: /issm/trunk/src/mex/BamgConvertMesh/BamgConvertMesh.cpp =================================================================== --- /issm/trunk/src/mex/BamgConvertMesh/BamgConvertMesh.cpp (revision 5176) +++ /issm/trunk/src/mex/BamgConvertMesh/BamgConvertMesh.cpp (revision 5177) @@ -67,4 +67,8 @@ bamgmesh->SetMatlabStructureFields(BAMGMESHOUT); + /*Clean up*/ + delete bamggeom; + delete bamgmesh; + /*end module: */ MODULEEND();