Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 3321) +++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 3322) @@ -355,5 +355,5 @@ BamgGeomInit(bamggeom); - //Vertices + /*Vertices*/ if(verbose>5) printf(" writing Vertices\n"); bamggeom->VerticesSize[0]=nbv; @@ -370,9 +370,6 @@ } } - else{ - bamggeom->Vertices=NULL; - } - - //Edges + + /*Edges*/ if(verbose>5) printf(" writing Edges\n"); bamggeom->EdgesSize[0]=nbe; @@ -394,9 +391,6 @@ } } - else{ - bamggeom->Edges=NULL; - } - - //CrackedEdges + + /*CrackedEdges*/ if(verbose>5) printf(" writing CrackedEdges\n"); bamggeom->CrackedEdgesSize[0]=nbcracked; @@ -415,9 +409,6 @@ } } - else{ - bamggeom->CrackedEdges=NULL; - } - - //RequiredEdges + + /*RequiredEdges*/ if(verbose>5) printf(" writing %i RequiredEdges\n",nbreq); bamggeom->RequiredEdgesSize[0]=nbreq; @@ -433,11 +424,8 @@ } } - else{ - bamggeom->RequiredEdges=NULL; - } //No corners - //RequiredVertices + /*RequiredVertices*/ if(verbose>5) printf(" writing %i RequiredVertices\n",nbreqv); bamggeom->RequiredVerticesSize[0]=nbreqv; @@ -453,9 +441,6 @@ } } - else{ - bamggeom->RequiredVertices=NULL; - } - - //SubDomains + + /*SubDomains*/ if(verbose>5) printf(" writing SubDomains\n"); bamggeom->SubDomainsSize[0]=NbSubDomains; @@ -470,9 +455,6 @@ } } - else{ - bamggeom->SubDomains=NULL; - } - - //TangentAtEdges + + /*TangentAtEdges*/ if(verbose>5) printf(" writing TangentAtEdges\n"); bamggeom->TangentAtEdgesSize[0]=nbtan; @@ -496,7 +478,4 @@ count=count+1; } - } - else{ - bamggeom->TangentAtEdges=NULL; } } Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3321) +++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3322) @@ -554,6 +554,12 @@ 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,9 +568,42 @@ 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=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;iconnectivitymax_1) connectivitymax_1=connectivitysize_1[i]; + } + + /*OK, now build outputs*/ + + /*Vertices*/ if(verbose>5) printf(" writing Vertices\n"); bamgmesh->VerticesSize[0]=nbv; @@ -578,9 +617,6 @@ } } - else{ - bamgmesh->Vertices=NULL; - } - - //Edges + + /*Edges*/ if(verbose>5) printf(" writing Edges\n"); bamgmesh->EdgesSize[0]=nbe; @@ -598,74 +634,6 @@ } } - 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=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=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;imax) max=connectivitysize[i]; - } - //Build output - bamgmesh->NodalElementConnectivitySize[0]=nbv; - bamgmesh->NodalElementConnectivitySize[1]=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 + + /*Element edges*/ if(verbose>5) printf(" writing element edges\n"); SetOfEdges4* edge4=new SetOfEdges4(nbt*3,nbv); @@ -707,5 +675,5 @@ xfree((void**)&elemedge); - //Segments + /*Segments*/ if(verbose>5) printf(" writing Segments\n"); bamgmesh->SegmentsSize[0]=NumSegments; @@ -719,5 +687,5 @@ 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,9 +717,6 @@ } } - //clean up - xfree((void**)&head_v); - xfree((void**)&next_p); - - //CrackedEdges + + /*CrackedEdges*/ if(verbose>5) printf(" writing CrackedEdges\n"); bamgmesh->CrackedEdgesSize[0]=NbCrackedEdges; @@ -764,9 +729,6 @@ } } - else{ - bamgmesh->CrackedEdges=NULL; - } - - //Triangles + + /*Triangles*/ if(verbose>5) printf(" writing Triangles\n"); k=nbInT-NbOfQuad*2; @@ -788,9 +750,6 @@ } } - else{ - bamgmesh->Triangles=NULL; - } - - //Quadrilaterals + + /*Quadrilaterals*/ if(verbose>5) printf(" writing Quadrilaterals\n"); bamgmesh->QuadrilateralsSize[0]=NbOfQuad; @@ -812,9 +771,6 @@ } } - else{ - bamgmesh->Quadrilaterals=NULL; - } - - //SubDomains + + /*SubDomains*/ if(verbose>5) printf(" writing SubDomains\n"); bamgmesh->SubDomainsSize[0]=NbSubDomains; @@ -829,9 +785,6 @@ } } - else{ - bamgmesh->SubDomains=NULL; - } - - //SubDomainsFromGeom + + /*SubDomainsFromGeom*/ if(verbose>5) printf(" writing SubDomainsFromGeom\n"); bamgmesh->SubDomainsFromGeomSize[0]=Gh.NbSubDomains; @@ -846,9 +799,6 @@ } } - else{ - bamgmesh->SubDomainsFromGeom=NULL; - } - - //VerticesOnGeomVertex + + /*VerticesOnGeomVertex*/ if(verbose>5) printf(" writing VerticesOnGeometricVertex\n"); bamgmesh->VerticesOnGeometricVertexSize[0]=NbVerticesOnGeomVertex; @@ -863,9 +813,6 @@ } } - else{ - bamgmesh->VerticesOnGeometricVertex=NULL; - } - - //VertexOnGeometricEdge + + /*VertexOnGeometricEdge*/ if(verbose>5) printf(" writing VerticesOnGeometricEdge\n"); bamgmesh->VerticesOnGeometricEdgeSize[0]=NbVerticesOnGeomEdge; @@ -883,9 +830,6 @@ } } - else{ - bamgmesh->VerticesOnGeometricEdge=NULL; - } - - //EdgesOnGeometricEdge + + /*EdgesOnGeometricEdge*/ if(verbose>5) printf(" writing EdgesOnGeometricEdge\n"); k=0; @@ -906,10 +850,99 @@ } } - 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=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;iNodalElementConnectivity[i]=NAN; + for (i=0;iNodalElementConnectivity[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;iElementEdges[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;iconnectivitymax_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;iNodalConnectivity[i]=NAN; + for (i=0;iElementEdges[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; }