Index: /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp (revision 3308) +++ /issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp (revision 3309) @@ -1,3 +1,2 @@ - #include #include "BamgObjects.h" @@ -39,7 +38,5 @@ //check that head is not empty - if (head==NULL) { - throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add no head is NULL (How come?)")); - } + assert(head); //get n from h (usually h=ii) Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3308) +++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3309) @@ -482,4 +482,5 @@ if(verbose>5) printf(" no Edges found\n"); } + //CrackedEdges if(bamgmesh->CrackedEdges){ @@ -554,5 +555,6 @@ void Triangles::WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){ - int i,j,k,num; + int i,j,k,num,i1,i2; + long n; /*Get options*/ int verbose=bamgopts->verbose; @@ -599,15 +601,34 @@ } - //Segments - bamgmesh->NumSegments=0; - if(verbose>5) printf(" writing Segments\n"); - - //chaining algorithm + //Element Connectivity + if(verbose>5) printf(" writing Element connectivity\n"); + bamgmesh->NumElementConnectivity=nbt-NbOutT; + bamgmesh->ElementConnectivity=(double*)xmalloc(3*(nbt-NbOutT)*sizeof(double)); + for (i=0;i<3*(nbt-NbOutT);i++) bamgmesh->ElementConnectivity[i]=NAN; + num=0; + for (i=0;i=0){ + for (j=0;j<3;j++){ + k=Number(triangles[i].TriangleAdj(j)); + if (reft[k]>=0){ + assert(3*num+j<3*(nbt-NbOutT)); + bamgmesh->ElementConnectivity[3*num+j]=k+1; // back to Matlab indexing + } + } + num+=1; + } + } + + //ElementNodal Connectivity + if(verbose>5) printf(" writing Nodal element connectivity\n"); + //chaining algorithm to get nodal connectivity int* head_v=NULL; head_v=(int*)xmalloc(nbv*sizeof(int)); int* next_p=NULL; next_p=(int*)xmalloc(3*nbt*sizeof(int)); - + int* connectivitysize=NULL; + connectivitysize=(int*)xmalloc(nbv*sizeof(int)); for (i=0;i3*nbt-1 || k<0) throw ErrorException(__FUNCT__,exprintf("k = %i, nbt = %i",k,nbt)); next_p[k]=head_v[v]; - if (v>nbv-1 || v<0) throw ErrorException(__FUNCT__,exprintf("v = %i, nbv = %i",v,nbv)); + if (v>nbv-1 || v<0) throw ErrorException(__FUNCT__,exprintf("v = %i, nbv = %i",v,nbv)); head_v[v]=k++; - } - } - } + connectivitysize[v]+=1; + } + } + } + int max=0; + for (i=0;imax) max=connectivitysize[i]; + } + //Build output + bamgmesh->NumNodalElementConnectivity=max; + bamgmesh->NodalElementConnectivity=(double*)xmalloc(max*nbv*sizeof(double)); + for (i=0;iNodalElementConnectivity[i]=NAN; + for (i=0;iNodalElementConnectivity[max*i+k]=floor(j/3)+1; + k++; + } + } + //clean up (do not erase head_v and next_p-> used and destroyed by segments) + xfree((void**)&connectivitysize); + + //Build Unique edges + if(verbose>5) printf(" writing all edges\n"); + SetOfEdges4* edge4=new SetOfEdges4(nbt*3,nbv); + double* elemedge=new double[3*nbt]; + for (i=0;i<3*nbt;i++) elemedge[i]=NAN; + k=0; + for (i=0;i=0){ + for (j=0;j<3;j++) { + i1=Number(triangles[i][VerticesOfTriangularEdge[j][0]]); + i2=Number(triangles[i][VerticesOfTriangularEdge[j][1]]); + n =edge4->SortAndFind(i1,i2); + if (n==-1){ + //first time + n=edge4->SortAndAdd(i1,i2); + elemedge[n*2+0]=double(k); + } + else{ + //second time + elemedge[n*2+1]=double(k); + } + } + k++; + } + } + bamgmesh->NumAllEdges=edge4->nb(); + bamgmesh->AllEdges=(double*)xmalloc(4*edge4->nb()*sizeof(double)); + for (i=0;inb();i++){ + bamgmesh->AllEdges[i*4+0]=edge4->i(i)+1;// back to M indexing + bamgmesh->AllEdges[i*4+1]=edge4->j(i)+1;// back to M indexing + bamgmesh->AllEdges[i*4+2]=elemedge[2*i+0]+1; // back to M indexing + bamgmesh->AllEdges[i*4+3]=elemedge[2*i+1]+1; // back to M indexing + } + //clean up + delete edge4; + delete elemedge; + + //Segments + bamgmesh->NumSegments=0; + if(verbose>5) printf(" writing Segments\n"); bamgmesh->NumSegments=NumSegments; bamgmesh->Segments=(double*)xmalloc(3*NumSegments*sizeof(double)); Index: /issm/trunk/src/c/objects/BamgGeom.cpp =================================================================== --- /issm/trunk/src/c/objects/BamgGeom.cpp (revision 3308) +++ /issm/trunk/src/c/objects/BamgGeom.cpp (revision 3309) @@ -1,4 +1,2 @@ -#ifdef _SERIAL_ -#include "mex.h" #include "stdio.h" #include "./BamgGeom.h" @@ -28,4 +26,5 @@ } +#ifdef _SERIAL_ void BamgGeomWrite(mxArray** pbamggeom_mat, BamgGeom* bamggeom){ Index: /issm/trunk/src/c/objects/BamgGeom.h =================================================================== --- /issm/trunk/src/c/objects/BamgGeom.h (revision 3308) +++ /issm/trunk/src/c/objects/BamgGeom.h (revision 3309) @@ -39,4 +39,5 @@ #ifdef _SERIAL_ +#include "mex.h" void BamgGeomWrite(mxArray** bamggeom_mat,BamgGeom* bamggeom); #endif Index: /issm/trunk/src/c/objects/BamgMesh.cpp =================================================================== --- /issm/trunk/src/c/objects/BamgMesh.cpp (revision 3308) +++ /issm/trunk/src/c/objects/BamgMesh.cpp (revision 3309) @@ -1,5 +1,3 @@ -#ifdef _SERIAL_ #include "stdio.h" -#include "mex.h" #include "./BamgMesh.h" @@ -12,4 +10,6 @@ bamgmesh->NumEdges=0; bamgmesh->Edges=NULL; + bamgmesh->NumAllEdges=0; + bamgmesh->AllEdges=NULL; bamgmesh->NumSegments=0; bamgmesh->Segments=NULL; @@ -30,7 +30,14 @@ bamgmesh->SubDomainsFromGeom=NULL; bamgmesh->hVertices=NULL; + bamgmesh->NumElementConnectivity=0; + bamgmesh->ElementConnectivity=NULL; + bamgmesh->NumNodalConnectivity=0; + bamgmesh->NodalConnectivity=NULL; + bamgmesh->NumNodalElementConnectivity=0; + bamgmesh->NodalElementConnectivity=NULL; } +#ifdef _SERIAL_ void BamgMeshWrite(mxArray** pbamgmesh_mat, BamgMesh* bamgmesh){ @@ -42,32 +49,40 @@ mxArray* pfield=NULL; mxArray* pfield2=NULL; - int numfields=23; + int numfields=31; const char* fnames[numfields]; mwSize ndim=2; mwSize dimensions[2]={1,1}; - fnames[0] = "NumTriangles"; - fnames[1] = "Triangles"; - fnames[2] = "NumVertices"; - fnames[3] = "Vertices"; - fnames[4] = "NumEdges"; - fnames[5] = "Edges"; - fnames[6] = "NumSegments"; - fnames[7] = "Segments"; - fnames[8] = "SegmentsMarkers"; - fnames[9] = "NumCrackedEdges"; - fnames[10] = "CrackedEdges"; - fnames[11] = "NumQuadrilaterals"; - fnames[12] = "Quadrilaterals"; - fnames[13] = "NumVerticesOnGeometricVertex"; - fnames[14] = "VerticesOnGeometricVertex"; - fnames[15] = "NumVerticesOnGeometricEdge"; - fnames[16] = "VerticesOnGeometricEdge"; - fnames[17] = "NumEdgesOnGeometricEdge"; - fnames[18] = "EdgesOnGeometricEdge"; - fnames[19] = "NumSubDomains"; - fnames[20] = "SubDomains"; - fnames[21] = "NumSubDomainsFromGeom"; - fnames[22] = "SubDomainsFromGeom"; + fnames[0] = "NumTriangles"; + fnames[1] = "Triangles"; + fnames[2] = "NumVertices"; + fnames[3] = "Vertices"; + fnames[4] = "NumEdges"; + fnames[5] = "Edges"; + fnames[6] = "NumSegments"; + fnames[7] = "Segments"; + fnames[8] = "NumAllEdges"; + fnames[9] = "AllEdges"; + fnames[10] = "SegmentsMarkers"; + fnames[11] = "NumCrackedEdges"; + fnames[12] = "CrackedEdges"; + fnames[13] = "NumQuadrilaterals"; + fnames[14] = "Quadrilaterals"; + fnames[15] = "NumVerticesOnGeometricVertex"; + fnames[16] = "VerticesOnGeometricVertex"; + fnames[17] = "NumVerticesOnGeometricEdge"; + fnames[18] = "VerticesOnGeometricEdge"; + fnames[19] = "NumEdgesOnGeometricEdge"; + fnames[20] = "EdgesOnGeometricEdge"; + fnames[21] = "NumSubDomains"; + fnames[22] = "SubDomains"; + fnames[23] = "NumSubDomainsFromGeom"; + fnames[24] = "SubDomainsFromGeom"; + fnames[25] = "NumElementConnectivity"; + fnames[26] = "ElementConnectivity"; + fnames[27] = "NumNodalConnectivity"; + fnames[28] = "NodalConnectivity"; + fnames[29] = "NumNodalElementConnectivity"; + fnames[30] = "NodalElementConnectivity"; bamgmesh_mat=mxCreateStructArray(ndim,dimensions,numfields,fnames); @@ -108,4 +123,13 @@ mexCallMATLAB(1,&pfield2,1,&pfield,"transpose");//transpose mxSetField(bamgmesh_mat,0,"Segments",pfield2); + + mxSetField(bamgmesh_mat,0,"NumAllEdges",mxCreateDoubleScalar(bamgmesh->NumAllEdges)); + + pfield=mxCreateDoubleMatrix(0,0,mxREAL); + mxSetM(pfield,4); + mxSetN(pfield,bamgmesh->NumAllEdges); + mxSetPr(pfield,bamgmesh->AllEdges); + mexCallMATLAB(1,&pfield2,1,&pfield,"transpose");//transpose + mxSetField(bamgmesh_mat,0,"AllEdges",pfield2); pfield=mxCreateDoubleMatrix(0,0,mxREAL); @@ -178,4 +202,31 @@ mexCallMATLAB(1,&pfield2,1,&pfield,"transpose");//transpose mxSetField(bamgmesh_mat,0,"SubDomainsFromGeom",pfield2); + + mxSetField(bamgmesh_mat,0,"NumElementConnectivity",mxCreateDoubleScalar(bamgmesh->NumElementConnectivity)); + + pfield=mxCreateDoubleMatrix(0,0,mxREAL); + mxSetM(pfield,3); + mxSetN(pfield,bamgmesh->NumElementConnectivity); + mxSetPr(pfield,bamgmesh->ElementConnectivity); + mexCallMATLAB(1,&pfield2,1,&pfield,"transpose");//transpose + mxSetField(bamgmesh_mat,0,"ElementConnectivity",pfield2); + + mxSetField(bamgmesh_mat,0,"NumNodalConnectivity",mxCreateDoubleScalar(bamgmesh->NumNodalConnectivity)); + + pfield=mxCreateDoubleMatrix(0,0,mxREAL); + mxSetM(pfield,2); + mxSetN(pfield,bamgmesh->NumNodalConnectivity); + mxSetPr(pfield,bamgmesh->NodalConnectivity); + mexCallMATLAB(1,&pfield2,1,&pfield,"transpose");//transpose + mxSetField(bamgmesh_mat,0,"NodalConnectivity",pfield2); + + mxSetField(bamgmesh_mat,0,"NumNodalElementConnectivity",mxCreateDoubleScalar(bamgmesh->NumNodalElementConnectivity)); + + pfield=mxCreateDoubleMatrix(0,0,mxREAL); + mxSetM(pfield,bamgmesh->NumNodalElementConnectivity); + mxSetN(pfield,bamgmesh->NumVertices); + mxSetPr(pfield,bamgmesh->NodalElementConnectivity); + mexCallMATLAB(1,&pfield2,1,&pfield,"transpose");//transpose + mxSetField(bamgmesh_mat,0,"NodalElementConnectivity",pfield2); /*Assign output pointer*/ Index: /issm/trunk/src/c/objects/BamgMesh.h =================================================================== --- /issm/trunk/src/c/objects/BamgMesh.h (revision 3308) +++ /issm/trunk/src/c/objects/BamgMesh.h (revision 3309) @@ -15,4 +15,7 @@ int NumEdges; double* Edges; + + int NumAllEdges; + double* AllEdges; int NumSegments; @@ -42,4 +45,13 @@ double* hVertices; + + int NumElementConnectivity; + double* ElementConnectivity; + + int NumNodalConnectivity; + double* NodalConnectivity; + + int NumNodalElementConnectivity; + double* NodalElementConnectivity; }; @@ -47,4 +59,5 @@ #ifdef _SERIAL_ +#include "mex.h" void BamgMeshWrite(mxArray** bamgmesh_mat,BamgMesh* bamgmesh); #endif Index: /issm/trunk/src/m/classes/@bamgmesh/bamgmesh.m =================================================================== --- /issm/trunk/src/m/classes/@bamgmesh/bamgmesh.m (revision 3308) +++ /issm/trunk/src/m/classes/@bamgmesh/bamgmesh.m (revision 3309) @@ -19,4 +19,7 @@ bm.NumQuadrilaterals=0; bm.Quadrilaterals=[]; + + bm.NumAllEdges=0; + bm.AllEdges=zeros(0,2); bm.NumSegments=0; @@ -42,4 +45,13 @@ bm.SubDomainsFromGeom=zeros(0,4); + bm.NumElementConnectivity=0; + bm.ElementConnectivity=zeros(0,3); + + bm.NumNodalConnectivity=0; + bm.NodalConnectivity=zeros(0,2); + + bm.NumNodalElementConnectivity=0; + bm.NodalElementConnectivity=[]; + bm.hVertices=[]; Index: /issm/trunk/src/mex/Bamg/Bamg.cpp =================================================================== --- /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 3308) +++ /issm/trunk/src/mex/Bamg/Bamg.cpp (revision 3309) @@ -50,4 +50,7 @@ FetchData(&bamgmesh_in.Triangles,NULL,NULL,mxGetField(BAMGMESH,0,"Triangles")); FetchData(&bamgmesh_in.hVertices,&lines,&cols,mxGetField(BAMGMESH,0,"hVertices")); + FetchData(&bamgmesh_in.NumSegments,mxGetField(BAMGMESH,0,"NumSegments")); + FetchData(&bamgmesh_in.Segments,NULL,NULL,mxGetField(BAMGMESH,0,"Segments")); + FetchData(&bamgmesh_in.SegmentsMarkers,NULL,NULL,mxGetField(BAMGMESH,0,"SegmentsMarkers")); FetchData(&bamgmesh_in.NumCrackedEdges,mxGetField(BAMGMESH,0,"NumCrackedEdges")); FetchData(&bamgmesh_in.CrackedEdges,NULL,NULL,mxGetField(BAMGMESH,0,"CrackedEdges"));