Changeset 3309

Timestamp:

03/19/10 15:13:26 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added very useful connectivity fields output from Bamg

Location:

issm/trunk/src

Files:

• c/Bamgx/objects/SetOfE4.cpp (modified) (2 diffs)
• c/Bamgx/objects/Triangles.cpp (modified) (4 diffs)
• c/objects/BamgGeom.cpp (modified) (2 diffs)
• c/objects/BamgGeom.h (modified) (1 diff)
• c/objects/BamgMesh.cpp (modified) (6 diffs)
• c/objects/BamgMesh.h (modified) (3 diffs)
• m/classes/@bamgmesh/bamgmesh.m (modified) (2 diffs)
• mex/Bamg/Bamg.cpp (modified) (1 diff)

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

@@ -1 @@
-
 #include <assert.h>
 #include "BamgObjects.h"
@@ -39 +38 @@
 
                 //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)

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

@@ -482 +482 @@
                         if(verbose>5) printf("      no Edges found\n");
                 }
+
                 //CrackedEdges
                 if(bamgmesh->CrackedEdges){
@@ -554 +555 @@
         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 +601 @@
                 }
 
-                //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<nbt;i++){
+                        if (reft[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;i<nbv;i++) head_v[i]=-1;
+                for (i=0;i<nbv;i++) connectivitysize[i]=0;
                 k=0;
                 for (i=0;i<nbt;i++) {
@@ -618 +639 @@
                                         if (k>3*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;i<nbv;i++){
+                        if (connectivitysize[i]>max) max=connectivitysize[i];
+                }
+                //Build output
+                bamgmesh->NumNodalElementConnectivity=max;
+                bamgmesh->NodalElementConnectivity=(double*)xmalloc(max*nbv*sizeof(double));
+                for (i=0;i<max*nbv;i++) bamgmesh->NodalElementConnectivity[i]=NAN;
+                for (i=0;i<nbv;i++){
+                        k=0;
+                        for(j=head_v[i];j!=-1;j=next_p[j]){
+                                assert(max*i+k < max*nbv);
+                                bamgmesh->NodalElementConnectivity[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<nbt;i++){
+                        //Do not take into account outside triangles (reft<0)
+                        if (reft[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;i<edge4->nb();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));

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

@@ -1 @@
-#ifdef _SERIAL_
-#include "mex.h"
 #include "stdio.h"
 #include "./BamgGeom.h"
@@ -28 +26 @@
 }
 
+#ifdef _SERIAL_
 void BamgGeomWrite(mxArray** pbamggeom_mat, BamgGeom* bamggeom){
 

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

@@ -39 +39 @@
 
 #ifdef _SERIAL_
+#include "mex.h"
 void BamgGeomWrite(mxArray** bamggeom_mat,BamgGeom* bamggeom);
 #endif

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

@@ -1 @@
-#ifdef _SERIAL_
 #include "stdio.h"
-#include "mex.h"
 #include "./BamgMesh.h"
 
@@ -12 +10 @@
         bamgmesh->NumEdges=0;
         bamgmesh->Edges=NULL;
+        bamgmesh->NumAllEdges=0;
+        bamgmesh->AllEdges=NULL;
         bamgmesh->NumSegments=0;
         bamgmesh->Segments=NULL;
@@ -30 +30 @@
         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 +49 @@
         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 +123 @@
         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 +202 @@
         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*/

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

@@ -15 +15 @@
         int     NumEdges;
         double* Edges;
+
+        int     NumAllEdges;
+        double* AllEdges;
 
         int     NumSegments;
@@ -42 +45 @@
 
         double* hVertices;
+
+        int     NumElementConnectivity;
+        double* ElementConnectivity;
+
+        int     NumNodalConnectivity;
+        double* NodalConnectivity;
+
+        int     NumNodalElementConnectivity;
+        double* NodalElementConnectivity;
 };
 
@@ -47 +59 @@
 
 #ifdef _SERIAL_
+#include "mex.h"
 void BamgMeshWrite(mxArray** bamgmesh_mat,BamgMesh* bamgmesh);
 #endif

issm/trunk/src/m/classes/@bamgmesh/bamgmesh.m

@@ -19 +19 @@
         bm.NumQuadrilaterals=0;
         bm.Quadrilaterals=[];
+
+        bm.NumAllEdges=0;
+        bm.AllEdges=zeros(0,2);
 
         bm.NumSegments=0;
@@ -42 +45 @@
         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=[];
 

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

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