[cig-commits] r4313 - in long/3D/Gale/trunk/src/StGermain: . Discretisation/Geometry/src

walter at geodynamics.org walter at geodynamics.org
Thu Aug 17 17:16:33 PDT 2006


Author: walter
Date: 2006-08-17 17:16:32 -0700 (Thu, 17 Aug 2006)
New Revision: 4313

Modified:
   long/3D/Gale/trunk/src/StGermain/
   long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Delaunay.c
   long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.c
Log:
 r2691 at earth:  boo | 2006-08-17 17:14:11 -0700
  r2645 at earth (orig r3726):  RaquibulHassan | 2006-07-28 02:36:42 -0700
  Fixed a bug in the routine for finding the neighbouring elements of each triangle.
  
  
 



Property changes on: long/3D/Gale/trunk/src/StGermain
___________________________________________________________________
Name: svk:merge
   - 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2690
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3725
   + 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2691
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3726

Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Delaunay.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Delaunay.c	2006-08-18 00:16:29 UTC (rev 4312)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Delaunay.c	2006-08-18 00:16:32 UTC (rev 4313)
@@ -166,7 +166,7 @@
 	self->idOffset = idOffset;
 	
 	if( sites != NULL ){
-		self->sites = Memory_Alloc_Array_Unnamed( Site, sizeof( Site ) * self->numSites );
+		self->sites = Memory_Alloc_Array_Unnamed( Site, self->numSites );
 		memset( self->boundingTriangle, 0, sizeof( self->boundingTriangle ) );
 
 		for( i=0; i<self->numSites; i++ ){
@@ -228,7 +228,7 @@
 	self->idOffset = idOffset;
 	
 	if( sites != NULL ){
-		self->sites = Memory_Alloc_Array_Unnamed( Site, sizeof( Site ) * self->numSites );
+		self->sites = Memory_Alloc_Array_Unnamed( Site, self->numSites );
 		memset( self->boundingTriangle, 0, sizeof( self->boundingTriangle ) );
 
 		for( i=0; i<self->numSites; i++ ){
@@ -333,6 +333,11 @@
 		Memory_Free( self->triangleIndices );
 	}
 
+	if( self->triangleNeighbours ){
+		Memory_Free( self->triangleNeighbours[0] );
+		Memory_Free( self->triangleNeighbours );
+	}
+
 	if( self->numNeighbours ){
 		Memory_Free( self->numNeighbours );
 	}
@@ -590,7 +595,7 @@
 	
 	start = le = delaunay->leftMost;
 	
-	delaunay->hull = Memory_Alloc_Array_Unnamed( int, sizeof( int ) * delaunay->numSites );
+	delaunay->hull = Memory_Alloc_Array_Unnamed( int, delaunay->numSites );
 	memset( delaunay->hull, 0, sizeof( int ) * delaunay->numSites );
 	
 	do{

Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.c	2006-08-18 00:16:29 UTC (rev 4312)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.c	2006-08-18 00:16:32 UTC (rev 4313)
@@ -351,17 +351,17 @@
 		}
 	}
 	
-	self->mapGlobalToLocal = Memory_Alloc_Array_Unnamed( int, sizeof(int)*numSites );
-	self->processorLoad = Memory_Alloc_Array_Unnamed( int, sizeof(int)*numProcs );
+	self->mapGlobalToLocal = Memory_Alloc_Array_Unnamed( int, numSites );
+	self->processorLoad = Memory_Alloc_Array_Unnamed( int, numProcs );
 	memset( self->processorLoad, 0, sizeof( int )*numProcs );
 
 	if( self->rank == MASTER_PROC ){
-		self->processor = Memory_Alloc_Array_Unnamed( int, sizeof(int)*numSites );
+		self->processor = Memory_Alloc_Array_Unnamed( int, numSites );
 
-		alloced = Memory_Alloc_Array_Unnamed( int, sizeof(int)*numSites );
+		alloced = Memory_Alloc_Array_Unnamed( int, numSites );
 		memset( alloced, 0, sizeof( int )*numSites );
 	
-		self->initialOrder = Memory_Alloc_Array_Unnamed( int, sizeof(int)*numSites );
+		self->initialOrder = Memory_Alloc_Array_Unnamed( int, numSites );
 		memset( self->initialOrder, 0, sizeof( int )*numSites );
 		
 		Delaunay_FindMinMax( self->sites, self->numSites, &_minX, &_minY, &_maxX, &_maxY );
@@ -397,7 +397,7 @@
 		}
 
 		self->numLocalSites = self->processorLoad[MASTER_PROC];
-		self->localPoints = Memory_Alloc_Array_Unnamed( CoordF, sizeof( CoordF )*self->numLocalSites );
+		self->localPoints = Memory_Alloc_Array_Unnamed( CoordF, self->numLocalSites );
 		
 		count = 0;
 		for( i=0; i<numSites; i++ ){
@@ -416,7 +416,7 @@
 		
 		MPI_Recv( &self->numLocalSites, 1, MPI_INT, MASTER_PROC, LOAD_TAG, *self->comm, &status );
 
-		self->localPoints = Memory_Alloc_Array_Unnamed( CoordF, sizeof( CoordF )*self->numLocalSites );
+		self->localPoints = Memory_Alloc_Array_Unnamed( CoordF, self->numLocalSites );
 
 		for( i=0; i<self->numLocalSites; i++ ){
 			MPI_Recv( &(self->localPoints[i]), sizeof(CoordF), MPI_BYTE, MASTER_PROC, DATA_TAG, *self->comm, &status );
@@ -466,7 +466,7 @@
 		self->localPoints = Memory_Realloc_Array
 			( self->localPoints, CoordF, sizeof(CoordF) * (self->numTotalLocalSites + self->haloSites[0]->nodeCount) );
 		
-		self->mappingTable[0] = Memory_Alloc_Array_Unnamed( int, sizeof( int ) * self->haloSites[0]->nodeCount );
+		self->mappingTable[0] = Memory_Alloc_Array_Unnamed( int, self->haloSites[0]->nodeCount );
 		memset( self->mappingTable[0], 0, sizeof( int ) * self->haloSites[0]->nodeCount );
 		
 		count=0;
@@ -515,7 +515,7 @@
 		self->localPoints = Memory_Realloc_Array
 			( self->localPoints, CoordF, sizeof(CoordF) * (self->numTotalLocalSites + self->haloSites[1]->nodeCount) );
 		
-		self->mappingTable[1] = Memory_Alloc_Array_Unnamed( int, sizeof( int ) * self->haloSites[1]->nodeCount );
+		self->mappingTable[1] = Memory_Alloc_Array_Unnamed( int, self->haloSites[1]->nodeCount );
 		memset( self->mappingTable[1], 0, sizeof( int ) * self->haloSites[1]->nodeCount );
 		
 		count=0;
@@ -587,8 +587,6 @@
 	int pass = 0;
 	int **triIndices = NULL;
 	Delaunay *delaunay = NULL;
-	int **edgeToTriangle = NULL;
-	int index = 0;
 
 	delaunay = pd->localTriangulation;
 
@@ -598,35 +596,10 @@
 	delaunay->triangleIndices[0] = Memory_Alloc_Array_Unnamed( int, delaunay->numFaces * 3 );
 	memset( delaunay->triangleIndices[0] , 0, sizeof(int) * delaunay->numFaces * 3 );
 
-	if( delaunay->attributes->BuildTriangleNeighbours ){
-		delaunay->triangleNeighbours = Memory_Alloc_Array_Unnamed( int*, delaunay->numFaces );
-		delaunay->triangleNeighbours[0] = Memory_Alloc_Array_Unnamed( int, delaunay->numFaces*3 );
-	
-		edgeToTriangle = Memory_Alloc_Array_Unnamed( int*, delaunay->qp->numQuadEdges );
-		edgeToTriangle[0] = Memory_Alloc_Array_Unnamed( int, delaunay->qp->numQuadEdges * 2 );
-	}
-	
 	for( i=0; i<delaunay->numFaces; i++ ){
 		delaunay->triangleIndices[i] = delaunay->triangleIndices[0]+i*3;
-
-		if ( delaunay->attributes->BuildTriangleNeighbours ){
-			delaunay->triangleNeighbours[i] = delaunay->triangleNeighbours[0]+i*3;
-			
-			delaunay->triangleNeighbours[i][0] = delaunay->numFaces;
-			delaunay->triangleNeighbours[i][1] = delaunay->numFaces;
-			delaunay->triangleNeighbours[i][2] = delaunay->numFaces;
-		}
 	}
 	
-	if( delaunay->attributes->BuildTriangleNeighbours ){
-		for( i=0; i<delaunay->qp->numQuadEdges; i++ ){
-			edgeToTriangle[i] = edgeToTriangle[0]+i*2;
-			
-			edgeToTriangle[i][0] = delaunay->numFaces;
-			edgeToTriangle[i][1] = delaunay->numFaces;
-		}
-	}
-	
 	triIndices = delaunay->triangleIndices;
 	
 	edges = delaunay->qp->quadEdges;
@@ -675,18 +648,6 @@
 							triIndices[triCount][1] = ParallelDelaunay_TranslateLocalToGlobal(pd, ((Site*)DEST(e))->id);
 							triIndices[triCount][2] = ParallelDelaunay_TranslateLocalToGlobal(pd, ((Site*)DEST(eOnext))->id);
 
-							if ( delaunay->attributes->BuildTriangleNeighbours ){
-								
-								index = (int)(((QuadEdge*)((void*)e)) - delaunay->qp->quadEdges);
-								edgeToTriangle[index][COUNT(e)] = triCount;
-							
-								index = (int)(((QuadEdge*)((void*)eOnext)) - delaunay->qp->quadEdges);
-								edgeToTriangle[index][COUNT(eOnext)] = triCount;
-								
-								index = (int)(((QuadEdge*)((void*)eLnext)) - delaunay->qp->quadEdges);
-								edgeToTriangle[index][COUNT(eLnext)] = triCount;
-							}
-
 							triCount++;
 						}
 					}
@@ -701,27 +662,6 @@
 	}
 	
 	delaunay->numTriangles = triCount;
-	
-	if( delaunay->attributes->BuildTriangleNeighbours ){
-		int *triangleNeighbourCount = NULL;
-
-		triangleNeighbourCount = Memory_Alloc_Array_Unnamed( int, delaunay->numFaces );
-		memset( triangleNeighbourCount, 0, sizeof( int ) * delaunay->numFaces );
-		
-		for( i=0; i<delaunay->qp->numQuadEdges; i++ ){
-			if( IS_FREE( (QuadEdgeRef)(&(delaunay->qp->quadEdges[i])) ) ) continue;
-
-		
-			if( edgeToTriangle[i][0] != delaunay->numFaces )
-				delaunay->triangleNeighbours[edgeToTriangle[i][0]][triangleNeighbourCount[edgeToTriangle[i][0]]++] = edgeToTriangle[i][1];
-		
-			if( edgeToTriangle[i][1] > delaunay->numFaces )
-				delaunay->triangleNeighbours[edgeToTriangle[i][1]][triangleNeighbourCount[edgeToTriangle[i][1]]++] = edgeToTriangle[i][0];
-		}
-		
-		Memory_Free( edgeToTriangle );
-		Memory_Free( triangleNeighbourCount );
-	}
 }
 
 #define NEIGHBOURS_TAG 1<<4
@@ -732,7 +672,7 @@
 
 void ParallelDelaunay_GatherTriangulation( ParallelDelaunay *pd )
 {
-	int i, j, count, count1;
+	int i, j, k, l, count, count1;
 	MPI_Status st;
 	int numNeighboursSum = 0;
 	int stride = 0;
@@ -743,11 +683,11 @@
 	if( pd->rank == MASTER_PROC ){
 		
 		if( pd->attributes->CalculateVoronoiSurfaceArea ){
-			pd->voronoiArea = Memory_Alloc_Array_Unnamed( float, pd->numInputSites * sizeof( float ) );
+			pd->voronoiArea = Memory_Alloc_Array_Unnamed( float, pd->numInputSites );
 		}
 
 		if( pd->attributes->FindNeighbours ){
-			pd->numNeighbours = Memory_Alloc_Array_Unnamed( int, pd->numInputSites * sizeof( int ) );
+			pd->numNeighbours = Memory_Alloc_Array_Unnamed( int, pd->numInputSites );
 		}
 		
 		count = 0;
@@ -775,13 +715,13 @@
 		}
 
 		if( pd->attributes->FindNeighbours ){
-			pd->neighbours = Memory_Alloc_Array_Unnamed( int*, sizeof(int*) * pd->numInputSites );
-			pd->neighbours[0] = Memory_Alloc_Array_Unnamed( int, sizeof(int) * numNeighboursSum );
+			pd->neighbours = Memory_Alloc_Array_Unnamed( int*, pd->numInputSites );
+			pd->neighbours[0] = Memory_Alloc_Array_Unnamed( int, numNeighboursSum );
 		}
 	
 		if( pd->attributes->CalculateVoronoiSides ){
-			pd->voronoiSides = Memory_Alloc_Array_Unnamed( float*, sizeof(float*) * pd->numInputSites );
-			pd->voronoiSides[0] = Memory_Alloc_Array_Unnamed( float, sizeof(float) * numNeighboursSum );
+			pd->voronoiSides = Memory_Alloc_Array_Unnamed( float*, pd->numInputSites );
+			pd->voronoiSides[0] = Memory_Alloc_Array_Unnamed( float, numNeighboursSum );
 		}
 		
 		stride = 0;
@@ -872,17 +812,17 @@
 	
 		pd->numTriangles = globalNumTriangles;
 		
-		triCountArray = Memory_Alloc_Array_Unnamed( int, sizeof(int)*pd->numProcs );
+		triCountArray = Memory_Alloc_Array_Unnamed( int, pd->numProcs );
 		rank = pd->rank;
 	
 		if( rank == MASTER_PROC ){
 
-			pd->triangleIndices = Memory_Alloc_Array_Unnamed( int*,  sizeof(int*) * globalNumTriangles );
-			pd->triangleIndices[0] = Memory_Alloc_Array_Unnamed( int, sizeof(int) * globalNumTriangles*3 );
+			pd->triangleIndices = Memory_Alloc_Array_Unnamed( int*,  globalNumTriangles );
+			pd->triangleIndices[0] = Memory_Alloc_Array_Unnamed( int, globalNumTriangles*3 );
 
 			if ( pd->attributes->BuildTriangleNeighbours ){
-				pd->triangleNeighbours = Memory_Alloc_Array_Unnamed( int*,  sizeof(int*) * globalNumTriangles );
-				pd->triangleNeighbours[0] = Memory_Alloc_Array_Unnamed( int, sizeof(int) * globalNumTriangles*3 );
+				pd->triangleNeighbours = Memory_Alloc_Array_Unnamed( int*,  globalNumTriangles );
+				pd->triangleNeighbours[0] = Memory_Alloc_Array_Unnamed( int, globalNumTriangles*3 );
 			}
 
 			triCountArray[0] = triCount;
@@ -895,25 +835,20 @@
 	
 				if ( pd->attributes->BuildTriangleNeighbours ){
 					pd->triangleNeighbours[i] = pd->triangleNeighbours[0] + i*3;
+
+					pd->triangleNeighbours[i][0] = pd->numTriangles;
+					pd->triangleNeighbours[i][1] = pd->numTriangles;
+					pd->triangleNeighbours[i][2] = pd->numTriangles;
 				}
 			}
 
 			for(j=0; j<triCountArray[0]; j++){
 				memcpy( pd->triangleIndices[j], triIndices[j], sizeof(int)*3 );
-				
-				if ( pd->attributes->BuildTriangleNeighbours ){
-					memcpy( pd->triangleNeighbours[j], pd->localTriangulation->triangleNeighbours[j], sizeof(int)*3 );
-				}
 			}
 		
 			for( i=1; i<pd->numProcs; i++ ){
 				for( count=0; count < triCountArray[i]; count++ ){
 					MPI_Recv( pd->triangleIndices[j], 3, MPI_INT, i, DATA_TAG, MPI_COMM_WORLD, &st );					
-					
-					if ( pd->attributes->BuildTriangleNeighbours ){
-						MPI_Recv( pd->triangleNeighbours[j], 3, MPI_INT, i, DATA_TAG, MPI_COMM_WORLD, &st );
-					}
-
 					j++;
 				}
 			}
@@ -927,15 +862,84 @@
 			/*for( j=0; j<globalNumTriangles; j++ ){
 				printf( "indices[%d] = [%d, %d, %d]\n", j, pd->triangleIndices[j][0], pd->triangleIndices[j][1], pd->triangleIndices[j][2] );
 			}*/
+			
+			stride = 0;
+			if ( pd->attributes->BuildTriangleNeighbours ){
+				int **nodeToTriangle = NULL;
+				int *nodeToTriangleCount = NULL;
+				int *triangleNeighboursCount = NULL;
+				int triangles[2];
+				int counter = 0;
+				int found = 0;
+				int m = 0;
+
+				nodeToTriangle = Memory_Alloc_Array_Unnamed( int*, pd->numInputSites );
+				nodeToTriangle[0] = Memory_Alloc_Array_Unnamed( int, numNeighboursSum );
+
+				nodeToTriangleCount = Memory_Alloc_Array_Unnamed( int, pd->numInputSites );
+				memset( nodeToTriangleCount, 0, sizeof( int ) * pd->numInputSites );
+				
+				triangleNeighboursCount = Memory_Alloc_Array_Unnamed( int, pd->numTriangles );
+				memset( triangleNeighboursCount, 0, sizeof( int ) * pd->numTriangles );
+
+				for( i=0; i<pd->numSites; i++ ){
+					nodeToTriangle[i] = nodeToTriangle[0] + stride;
+
+					stride += pd->numNeighbours[i];
+				}
+
+				for( i=0; i<pd->numTriangles; i++ ){
+					nodeToTriangle[pd->triangleIndices[i][0]][nodeToTriangleCount[pd->triangleIndices[i][0]]++] = i;
+					nodeToTriangle[pd->triangleIndices[i][1]][nodeToTriangleCount[pd->triangleIndices[i][1]]++] = i;
+					nodeToTriangle[pd->triangleIndices[i][2]][nodeToTriangleCount[pd->triangleIndices[i][2]]++] = i;
+				}
+				
+				for( i=0; i<pd->numInputSites; i++ ){
+					for( j=0; j<pd->numNeighbours[i]; j++ ){
+						triangles[0] = -1;
+						triangles[1] = -1;
+						counter = 0;
+						for( k=0; k<nodeToTriangleCount[i]; k++ ){
+							for( l=0; l<3; l++ ){
+								if( pd->neighbours[i][j] == pd->triangleIndices[nodeToTriangle[i][k]][l] ){
+									triangles[counter++] = nodeToTriangle[i][k];
+								}
+							}
+						}
+						
+						if( (triangles[0] > -1) ){
+							
+							found = 0;
+							for( m=0; m<3; m++ ){
+								if( pd->triangleNeighbours[triangles[0]][m] == triangles[1] ) found = 1;
+							}
+							if( !found ){
+								pd->triangleNeighbours[triangles[0]][triangleNeighboursCount[triangles[0]]++] = triangles[1];
+							}
+						}
+						if( (triangles[1] > -1) ){
+						
+							found = 0;
+							for( m=0; m<3; m++ ){
+								if( pd->triangleNeighbours[triangles[1]][m] == triangles[0] ) found = 1;
+							}
+							if( !found ){
+								pd->triangleNeighbours[triangles[1]][triangleNeighboursCount[triangles[1]]++] = triangles[0];
+							}
+						}
+					}
+				}
+				
+				Memory_Free( nodeToTriangle[0] );
+				Memory_Free( nodeToTriangle );
+				Memory_Free( nodeToTriangleCount );
+				Memory_Free( triangleNeighboursCount );
+			}
 		}
 		else{
 			MPI_Send( &triCount, 1, MPI_INT, MASTER_PROC, DATA_TAG, MPI_COMM_WORLD );
 			for( j=0; j<triCount; j++ ){
 				MPI_Send( triIndices[j], 3, MPI_INT, MASTER_PROC, DATA_TAG, MPI_COMM_WORLD );
-				
-				if ( pd->attributes->BuildTriangleNeighbours ){
-					MPI_Send( pd->localTriangulation->triangleNeighbours[j], 3, MPI_INT, MASTER_PROC, DATA_TAG, MPI_COMM_WORLD );
-				}
 			}
 		}
 
@@ -1006,10 +1010,10 @@
 	
 	assert(edge);
 
-	s = Memory_Alloc_Array_Unnamed( Site, sizeof( Site ) * 2 );
+	s = Memory_Alloc_Array_Unnamed( Site, 2 );
 	memset( s, 0, sizeof( Site ) * 2 );
 	
-	c = Memory_Alloc_Array_Unnamed( CoordF, sizeof( CoordF ) * 2 );
+	c = Memory_Alloc_Array_Unnamed( CoordF, 2 );
 	memset( c, 0, sizeof( CoordF ) * 2 );
 	
 	MPI_Recv( &(sp[0]), sizeof(SitePacket), MPI_BYTE, rank, ORG_TAG, *comm, &st );



More information about the cig-commits mailing list