Changeset 3322

Timestamp:

03/23/10 08:26:51 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added nodal connectivity

Location:

issm/trunk/src/c/Bamgx/objects

Files:

• Geometry.cpp (modified) (8 diffs)
• Triangles.cpp (modified) (15 diffs)

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

@@ -355 +355 @@
                 BamgGeomInit(bamggeom);
 
-                //Vertices
+                /*Vertices*/
                 if(verbose>5) printf("      writing Vertices\n");
                 bamggeom->VerticesSize[0]=nbv;
@@ -370 +370 @@
                         }
                 }
-                else{
-                        bamggeom->Vertices=NULL;
-                }
-
-                //Edges
+
+                /*Edges*/
                 if(verbose>5) printf("      writing Edges\n");
                 bamggeom->EdgesSize[0]=nbe;
@@ -394 +391 @@
                         }
                 }
-                else{
-                        bamggeom->Edges=NULL;
-                }
-
-                //CrackedEdges
+
+                /*CrackedEdges*/
                 if(verbose>5) printf("      writing CrackedEdges\n");
                 bamggeom->CrackedEdgesSize[0]=nbcracked;
@@ -415 +409 @@
                         }
                 }
-                else{
-                        bamggeom->CrackedEdges=NULL;
-                }
-
-                //RequiredEdges
+
+                /*RequiredEdges*/
                 if(verbose>5) printf("      writing %i RequiredEdges\n",nbreq);
                 bamggeom->RequiredEdgesSize[0]=nbreq;
@@ -433 +424 @@
                         }
                 }
-                else{
-                        bamggeom->RequiredEdges=NULL;
-                }
 
                 //No corners
 
-                //RequiredVertices
+                /*RequiredVertices*/
                 if(verbose>5) printf("      writing %i RequiredVertices\n",nbreqv);
                 bamggeom->RequiredVerticesSize[0]=nbreqv;
@@ -453 +441 @@
                         }
                 }
-                else{
-                        bamggeom->RequiredVertices=NULL;
-                }
-
-                //SubDomains
+
+                /*SubDomains*/
                 if(verbose>5) printf("      writing SubDomains\n");
                 bamggeom->SubDomainsSize[0]=NbSubDomains;
@@ -470 +455 @@
                         }
                 }
-                else{
-                        bamggeom->SubDomains=NULL;
-                }
-
-                //TangentAtEdges
+
+                /*TangentAtEdges*/
                 if(verbose>5) printf("      writing TangentAtEdges\n");
                 bamggeom->TangentAtEdgesSize[0]=nbtan;
@@ -496 +478 @@
                                 count=count+1;
                         }
-                }
-                else{
-                        bamggeom->TangentAtEdges=NULL;
                 }
         }

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

@@ -554 +554 @@
         void Triangles::WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts){
 
+                /*Intermediary*/
                 int i,j,k,num,i1,i2;
                 long n;
+                int* head_1=NULL;
+                int* next_1=NULL;
+                int* connectivitysize_1=NULL;
+                int  connectivitymax_1=0;
+
                 /*Get options*/
                 int verbose=bamgopts->verbose;
@@ -562 +568 @@
                 BamgMeshInit(bamgmesh);
 
-                //Build reft that holds the number the subdomain number of each triangle
+                /*Build reft that holds the number the subdomain number of each triangle*/
                 long* reft = new long[nbt];
                 long nbInT = TriangleReferenceList(reft);
 
-                //Vertices
+                /*Chaining algorithm used to generate connectivity tables and other outputs*/
+
+                //Memory Allocation
+                head_1=(int*)xmalloc(nbv*sizeof(int));
+                next_1=(int*)xmalloc(3*nbt*sizeof(int));
+                connectivitysize_1=(int*)xmalloc(nbv*sizeof(int));
+
+                //Initialization
+                for (i=0;i<nbv;i++) head_1[i]=-1;
+                for (i=0;i<nbv;i++) connectivitysize_1[i]=0;
+                k=0;
+                //Chains generation
+                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++){
+                                        int v=Number(triangles[i][j]); //jth vertex of the ith triangle
+                                        if (k>3*nbt-1 || k<0) throw ErrorException(__FUNCT__,exprintf("k = %i, nbt = %i",k,nbt));
+                                        next_1[k]=head_1[v];
+                                        if (v>nbv-1 || v<0)   throw ErrorException(__FUNCT__,exprintf("v = %i, nbv = %i",v,nbv));
+                                        head_1[v]=k++;
+                                        connectivitysize_1[v]+=1;
+                                }
+                        }
+                }
+                //Get maximum connectivity
+                connectivitymax_1=0;
+                for (i=0;i<nbv;i++){
+                        if (connectivitysize_1[i]>connectivitymax_1) connectivitymax_1=connectivitysize_1[i];
+                }
+
+                /*OK, now build outputs*/
+
+                /*Vertices*/
                 if(verbose>5) printf("      writing Vertices\n");
                 bamgmesh->VerticesSize[0]=nbv;
@@ -578 +617 @@
                         }
                 }
-                else{
-                        bamgmesh->Vertices=NULL;
-                }
-
-                //Edges
+
+                /*Edges*/
                 if(verbose>5) printf("      writing Edges\n");
                 bamgmesh->EdgesSize[0]=nbe;
@@ -598 +634 @@
                         }
                 }
-                else{
-                        bamgmesh->Edges=NULL;
-                }
-
-                //Element Connectivity
-                if(verbose>5) printf("      writing Element connectivity\n");
-                bamgmesh->ElementConnectivitySize[0]=nbt-NbOutT;
-                bamgmesh->ElementConnectivitySize[1]=3;
-                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++) {
-                        //Do not take into account outside triangles (reft<0)
-                        if (reft[i]>=0){
-                                for (j=0;j<3;j++){
-                                        int v=Number(triangles[i][j]); //jth vertex of the ith triangle
-                                        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));
-                                        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->NodalElementConnectivitySize[0]=nbv;
-                bamgmesh->NodalElementConnectivitySize[1]=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
+
+                /*Element edges*/
                 if(verbose>5) printf("      writing element edges\n");
                 SetOfEdges4* edge4=new SetOfEdges4(nbt*3,nbv);
@@ -707 +675 @@
                 xfree((void**)&elemedge);
 
-                //Segments
+                /*Segments*/
                 if(verbose>5) printf("      writing Segments\n");
                 bamgmesh->SegmentsSize[0]=NumSegments;
@@ -719 +687 @@
                                 int i2=Number(edges[i][1]);
                                 bool stop=false;
-                                for(j=head_v[i1];j!=-1;j=next_p[j]){
+                                for(j=head_1[i1];j!=-1;j=next_1[j]){
                                         for(k=0;k<3;k++){
                                                 if (Number(triangles[(int)j/3][k])==i1){
@@ -749 +717 @@
                         }
                 }
-                //clean up
-                xfree((void**)&head_v);
-                xfree((void**)&next_p);
-
-                //CrackedEdges
+
+                /*CrackedEdges*/
                 if(verbose>5) printf("      writing CrackedEdges\n");
                 bamgmesh->CrackedEdgesSize[0]=NbCrackedEdges;
@@ -764 +729 @@
                         }
                 }
-                else{
-                        bamgmesh->CrackedEdges=NULL;
-                }
-
-                //Triangles
+
+                /*Triangles*/
                 if(verbose>5) printf("      writing Triangles\n");
                 k=nbInT-NbOfQuad*2;
@@ -788 +750 @@
                         }
                 }
-                else{
-                        bamgmesh->Triangles=NULL;
-                }
-
-                //Quadrilaterals
+
+                /*Quadrilaterals*/
                 if(verbose>5) printf("      writing Quadrilaterals\n");
                 bamgmesh->QuadrilateralsSize[0]=NbOfQuad;
@@ -812 +771 @@
                         }
                 }
-                else{
-                        bamgmesh->Quadrilaterals=NULL;
-                }
-
-                //SubDomains
+
+                /*SubDomains*/
                 if(verbose>5) printf("      writing SubDomains\n");
                 bamgmesh->SubDomainsSize[0]=NbSubDomains;
@@ -829 +785 @@
                         }
                 }
-                else{
-                        bamgmesh->SubDomains=NULL;
-                }
-
-                //SubDomainsFromGeom
+
+                /*SubDomainsFromGeom*/
                 if(verbose>5) printf("      writing SubDomainsFromGeom\n");
                 bamgmesh->SubDomainsFromGeomSize[0]=Gh.NbSubDomains;
@@ -846 +799 @@
                         }
                 }
-                else{
-                        bamgmesh->SubDomainsFromGeom=NULL;
-                }
-
-                //VerticesOnGeomVertex
+
+                /*VerticesOnGeomVertex*/
                 if(verbose>5) printf("      writing VerticesOnGeometricVertex\n");
                 bamgmesh->VerticesOnGeometricVertexSize[0]=NbVerticesOnGeomVertex;
@@ -863 +813 @@
                         }
                 }
-                else{
-                        bamgmesh->VerticesOnGeometricVertex=NULL;
-                }
-
-                //VertexOnGeometricEdge
+
+                /*VertexOnGeometricEdge*/
                 if(verbose>5) printf("      writing VerticesOnGeometricEdge\n");
                 bamgmesh->VerticesOnGeometricEdgeSize[0]=NbVerticesOnGeomEdge;
@@ -883 +830 @@
                         }
                 }
-                else{
-                        bamgmesh->VerticesOnGeometricEdge=NULL;
-                }
-
-                //EdgesOnGeometricEdge
+
+                /*EdgesOnGeometricEdge*/
                 if(verbose>5) printf("      writing EdgesOnGeometricEdge\n");
                 k=0;
@@ -906 +850 @@
                         }
                 }
-                else{
-                        bamgmesh->EdgesOnGeometricEdge=NULL;
-                }
-
+
+                /*Element Connectivity*/
+                if(verbose>5) printf("      writing Element connectivity\n");
+                bamgmesh->ElementConnectivitySize[0]=nbt-NbOutT;
+                bamgmesh->ElementConnectivitySize[1]=3;
+                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");
+                bamgmesh->NodalElementConnectivitySize[0]=nbv;
+                bamgmesh->NodalElementConnectivitySize[1]=connectivitymax_1;
+                bamgmesh->NodalElementConnectivity=(double*)xmalloc(connectivitymax_1*nbv*sizeof(double));
+                for (i=0;i<connectivitymax_1*nbv;i++) bamgmesh->NodalElementConnectivity[i]=NAN;
+                for (i=0;i<nbv;i++){
+                        k=0;
+                        for(j=head_1[i];j!=-1;j=next_1[j]){
+                                assert(connectivitymax_1*i+k < connectivitymax_1*nbv);
+                                bamgmesh->NodalElementConnectivity[connectivitymax_1*i+k]=floor(j/3)+1;
+                                k++;
+                        }
+                }
+
+                /*Nodal Connectivity*/
+                if(verbose>5) printf("      writing Nodal connectivity\n");
+                //chaining algorithm (again...)
+                int* head_2=NULL;
+                int* next_2=NULL;
+                int* connectivitysize_2=NULL;
+                int  connectivitymax_2=0;
+                i1=bamgmesh->ElementEdgesSize[0];
+                i2=bamgmesh->ElementEdgesSize[1];
+                head_2=(int*)xmalloc(nbv*sizeof(int));
+                next_2=(int*)xmalloc(2*i1*sizeof(int));
+                connectivitysize_2=(int*)xmalloc(nbv*sizeof(int));
+                //Initialization
+                for (i=0;i<nbv;i++) head_2[i]=-1;
+                for (i=0;i<nbv;i++) connectivitysize_2[i]=0;
+                k=0;
+                //Chains generation
+                for (i=0;i<i1;i++) {
+                        for (j=0;j<2;j++){
+                                int v=(int)bamgmesh->ElementEdges[i*i2+j]-1; //back to C indexing
+                                if (k>2*i1-1 || k<0) throw ErrorException(__FUNCT__,exprintf("Index exceed matrix dimensions (k=%i not in [0 %i]",k,2*i1-1));
+                                next_2[k]=head_2[v];
+                                if (v>nbv-1 || v<0)   throw ErrorException(__FUNCT__,exprintf("Index exceed matrix dimensions (v=%i not in [0 %i])",v,nbv-1));
+                                head_2[v]=k++;
+                                connectivitysize_2[v]+=1;
+                        }
+                }
+                //Get maximum connectivity
+                for (i=0;i<nbv;i++){
+                        if (connectivitysize_2[i]>connectivitymax_2) connectivitymax_2=connectivitysize_2[i];
+                }
+                //Build output
+                bamgmesh->NodalConnectivitySize[0]=nbv;
+                bamgmesh->NodalConnectivitySize[1]=connectivitymax_2;
+                bamgmesh->NodalConnectivity=(double*)xmalloc(connectivitymax_2*nbv*sizeof(double));
+                for (i=0;i<connectivitymax_2*nbv;i++) bamgmesh->NodalConnectivity[i]=NAN;
+                for (i=0;i<nbv;i++){
+                        k=0;
+                        for(j=head_2[i];j!=-1;j=next_2[j]){
+                                assert(connectivitymax_2*i+k < connectivitymax_2*nbv);
+                                num=(int)bamgmesh->ElementEdges[int(j/2)*i2+0];
+                                if (i+1==num){ //carefull, ElementEdge is in M indexing
+                                        //i is the first vertex of the edge, it is therefore connected to the second vertex
+                                        bamgmesh->NodalConnectivity[connectivitymax_2*i+k]=bamgmesh->ElementEdges[int(j/2)*i2+1];
+                                }
+                                else{
+                                        bamgmesh->NodalConnectivity[connectivitymax_2*i+k]=num;
+                                }
+                                k++;
+                        }
+                }
 
                 //clean up
+                xfree((void**)&connectivitysize_1);
+                xfree((void**)&head_1);
+                xfree((void**)&next_1);
+                xfree((void**)&connectivitysize_2);
+                xfree((void**)&head_2);
+                xfree((void**)&next_2);
                 delete [] reft;
         }