[cig-commits] r4044 - in long/3D/Gale/trunk/src/StGermain: .
Discretisation/Geometry/src
walter at geodynamics.org
walter at geodynamics.org
Thu Jul 20 20:03:39 PDT 2006
Author: walter
Date: 2006-07-20 20:03:38 -0700 (Thu, 20 Jul 2006)
New Revision: 4044
Modified:
long/3D/Gale/trunk/src/StGermain/
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.c
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.h
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/QuadEdge.c
Log:
r2530 at earth: boo | 2006-07-20 20:01:36 -0700
r2506 at earth (orig r3674): RaquibulHassan | 2006-07-11 20:33:14 -0700
* Implemented a function that gathers 'n' sub-triangulations on the master processor.
* Implemented mapping tables to translate local node numbers to that of global.
* Implemented a function to generate the global triangle indices.
Property changes on: long/3D/Gale/trunk/src/StGermain
___________________________________________________________________
Name: svk:merge
- 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2498
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3673
+ 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2530
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3674
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.c 2006-07-21 03:03:26 UTC (rev 4043)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.c 2006-07-21 03:03:38 UTC (rev 4044)
@@ -61,6 +61,7 @@
const Type ParallelDelaunay_Type="ParallelDelaunay";
#define PI 3.1415926535897932384626
+#define MERGE_FACTOR 1
/*--------------------------------------------------------------------------------------------------------------------------
** Constructors
@@ -158,16 +159,12 @@
self->_execute = _ParallelDelaunay_Execute;
self->_destroy = _ParallelDelaunay_Destroy;
- self->attributes = attr;
+ self->attributes = Memory_Alloc_Unnamed( DelaunayAttributes );
+ memcpy( self->attributes, attr, sizeof( DelaunayAttributes ) );
+ self->attributes->BuildBoundingTriangle = 0;
+
self->dictionary = dictionary;
-
- if( self->attributes->BuildBoundingTriangle ){
- self->numSites = numSites + 3;
- }
- else{
- self->numSites = numSites;
- }
-
+ self->numSites = numSites;
self->numInputSites = numSites;
self->points = sites;
self->leftProc = 0;
@@ -209,10 +206,15 @@
DelaunayAttributes *attr )
{
ParallelDelaunay *self = NULL;
+ DelaunayAttributes *myAttr = NULL;
+ myAttr = Memory_Alloc_Unnamed( DelaunayAttributes );
+ memcpy( myAttr, attr, sizeof( DelaunayAttributes ) );
+ myAttr->BuildBoundingTriangle = 0;
+
assert( _sizeOfSelf >= sizeof(ParallelDelaunay) );
self = (ParallelDelaunay*)_Delaunay_New( _sizeOfSelf, type, _delete, _print, _copy, _defaultConstructor,
- _construct, _build, _initialise, _execute, _destroy, name, True, dictionary, sites, numSites, 0, attr );
+ _construct, _build, _initialise, _execute, _destroy, name, True, dictionary, sites, numSites, 0, myAttr );
self->points = sites;
self->leftProc = 0;
@@ -279,7 +281,16 @@
assert( self );
- _Stg_Component_Delete( self );
+ Memory_Free( self->localPoints );
+ if( self->mappingTable[0] ) Memory_Free( self->mappingTable[0] );
+ if( self->mappingTable[1] ) Memory_Free( self->mappingTable[1] );
+ Memory_Free( self->mapGlobalToLocal );
+ if( self->processor ) Memory_Free( self->processor );
+ Memory_Free( self->processorLoad );
+ Memory_Free( self->attributes );
+
+ Stg_Class_Delete( self->localTriangulation );
+ _Delaunay_Delete( self );
}
/** Stg_Class_Print() implementation */
@@ -413,7 +424,7 @@
}
self->localTriangulation = Delaunay_New( "delaunay", self->dictionary, self->localPoints, self->numLocalSites, offset, self->attributes );
- self->localTriangulation->qp = QuadEdgePool_New( self->localTriangulation->numSites * 4 );
+ self->localTriangulation->qp = QuadEdgePool_New( self->localTriangulation->numSites * (3 + MERGE_FACTOR) );
Delaunay_SortSites(self->localTriangulation->sites, self->localTriangulation->numSites);
Delaunay_Recurse(self->localTriangulation, 0, self->localTriangulation->numSites,
&self->localTriangulation->leftMost, &self->localTriangulation->rightMost);
@@ -528,68 +539,322 @@
self->localTriangulation = Delaunay_New( "delaunay", self->dictionary, self->localPoints, self->numTotalLocalSites, offset, self->attributes );
Stg_Component_Build( self->localTriangulation, NULL, True );
- for( i=0; i<self->numTotalLocalSites; i++ ){
- if( (self->localTriangulation->sites[i].id) >= (self->localTriangulation->idOffset+self->numLocalSites) ){
+ if( self->attributes->BuildTriangleIndices ){
+ ParallelDelaunay_BuildTriangleIndices( self );
+ }
+}
+
+void _ParallelDelaunay_Initialise( void* pd, void* data )
+{
- if( self->localTriangulation->sites[i].id <
- (self->localTriangulation->idOffset+self->numLocalSites + self->numHaloSites[0]) ){
+}
+
+void _ParallelDelaunay_Execute( void* pd, void* data )
+{
+
+}
+
+void _ParallelDelaunay_Destroy( void* pd, void* data )
+{
+
+}
+
+#define onCurrentProc( pd, id ) ( id < (pd->processorLoad[pd->rank]+pd->localTriangulation->idOffset) )
+#define onLeftProc( pd, id ) ( ( pd->numHaloSites[0]>0 ) && ( id >= (pd->processorLoad[pd->rank]+pd->localTriangulation->idOffset) ) && \
+ ( id < (pd->processorLoad[pd->rank]+pd->numHaloSites[0]+pd->localTriangulation->idOffset) ) )
+#define onRightProc( pd, id ) ( ( pd->numHaloSites[1]>0 ) && ( id >= (pd->processorLoad[pd->rank]+pd->localTriangulation->idOffset+pd->numHaloSites[0]) ) && \
+ ( id < (pd->numTotalLocalSites+pd->localTriangulation->idOffset) ) )
+
+void ParallelDelaunay_BuildTriangleIndices( ParallelDelaunay *pd )
+{
+ int i = 0, triCount;
+ QuadEdgeRef e = 0, eStart = 0, eOnext = 0, eLnext = 0;
+ QuadEdge *edges = NULL;
+ Site *sites = NULL;
+ int maxEdges = 0;
+ int rank = 0;
+ int pass = 0;
+ int **triIndices = NULL;
+ Delaunay *delaunay = NULL;
+
+ delaunay = pd->localTriangulation;
+
+ assert( delaunay );
+
+ delaunay->triangleIndices = Memory_Alloc_Array_Unnamed( int*, delaunay->numFaces );
+ delaunay->triangleIndices[0] = Memory_Alloc_Array_Unnamed( int, delaunay->numFaces * 3 );
+ memset( delaunay->triangleIndices[0] , 0, sizeof(int) * delaunay->numFaces * 3 );
+
+ for( i=0; i<delaunay->numFaces; i++ ){
+ delaunay->triangleIndices[i] = delaunay->triangleIndices[0]+i*3;
+ }
+
+ triIndices = delaunay->triangleIndices;
+
+ edges = delaunay->qp->quadEdges;
+ sites = delaunay->sites;
+ rank = pd->rank;
+ maxEdges = delaunay->qp->numQuadEdges;
+
+ for (i = 0; i < maxEdges; i++) {
+ edges[i].count = 0;
+ }
+
+ triCount = 0;
+ for (i = 0; i < maxEdges; i++) {
+
+ e = eStart = (QuadEdgeRef)((void*)&(edges[i]));
+
+ if( IS_FREE(e) )continue;
+
+ do{
+ eOnext = ONEXT(e);
+ eLnext = LNEXT(e);
- self->localTriangulation->sites[i].id =
- self->mappingTable[0][(self->localTriangulation->sites[i].id -
- (self->localTriangulation->idOffset+self->numLocalSites))];
+ if( (COUNT(e)<2) && (COUNT(LNEXT(e))<2) && (COUNT(eOnext)<2) ){
+ if( ((((Site*)ORG(eLnext)) == ((Site*)DEST(e)))) &&
+ ((((Site*)DEST(eLnext)) == ((Site*)DEST(eOnext)))) ){
+
+ if( onCurrentProc( pd, ((Site*)ORG(e))->id ) ||
+ onCurrentProc( pd, ((Site*)DEST(e))->id ) ||
+ onCurrentProc( pd, ((Site*)DEST(eOnext))->id ) ){
+
+ pass = 0;
+ if( !onCurrentProc( pd, ((Site*)ORG(e))->id ) ){
+ if( !onRightProc( pd, ((Site*)ORG(e))->id ) ) pass = 1;
+ }
+
+ if( !onCurrentProc( pd, ((Site*)DEST(e))->id ) ){
+ if( !onRightProc( pd, ((Site*)DEST(e))->id ) ) pass = 1;
+ }
+
+ if( !onCurrentProc( pd, ((Site*)DEST(eOnext))->id ) ){
+ if( !onRightProc( pd, ((Site*)DEST(eOnext))->id ) ) pass = 1;
+ }
+
+ if( !pass ){
+ triIndices[triCount][0] = ParallelDelaunay_TranslateLocalToGlobal(pd, ((Site*)ORG(e))->id);
+ triIndices[triCount][1] = ParallelDelaunay_TranslateLocalToGlobal(pd, ((Site*)DEST(e))->id);
+ triIndices[triCount][2] = ParallelDelaunay_TranslateLocalToGlobal(pd, ((Site*)DEST(eOnext))->id);
+ triCount++;
+ }
+ }
+
+ COUNT(e)++;
+ COUNT(LNEXT(e))++;
+ COUNT(eOnext)++;
+ }
}
- else{
- self->localTriangulation->sites[i].id =
- self->mappingTable[1][(self->localTriangulation->sites[i].id -
- (self->localTriangulation->idOffset+
- self->numLocalSites+
- self->numHaloSites[0]))];
- }
+ e = eOnext;
+ }while( e != eStart );
+ }
+
+ delaunay->numTriangles = triCount;
+}
+
+#define NEIGHBOURS_TAG 1<<4
+#define VORONOI_SIDES_TAG 1<<5
+#define VORONOI_AREA_TAG 1<<6
+#define NUM_NEIGHBOUR_TAG 1<<7
+#define MAX_NEIGHBOURS 100
+void ParallelDelaunay_GatherTriangulation( ParallelDelaunay *pd )
+{
+ int i, j, count, count1;
+ MPI_Status st;
+ int numNeighboursSum = 0;
+ int stride = 0;
+ int translationArray[MAX_NEIGHBOURS] = { 0 };
+
+ assert( pd );
+
+ if( pd->rank == MASTER_PROC ){
+
+ if( pd->attributes->CalculateVoronoiSurfaceArea ){
+ pd->voronoiArea = Memory_Alloc_Array_Unnamed( float, pd->numInputSites * sizeof( float ) );
}
+
+ if( pd->attributes->FindNeighbours ){
+ pd->numNeighbours = Memory_Alloc_Array_Unnamed( int, pd->numInputSites * sizeof( int ) );
+ }
- for( count=0; count<self->localTriangulation->numNeighbours[i]; count++ ){
- if( (self->localTriangulation->neighbours[i][count]) >= (self->localTriangulation->idOffset+self->numLocalSites) ){
- if( self->localTriangulation->neighbours[i][count] <
- (self->localTriangulation->idOffset+self->numLocalSites + self->numHaloSites[0]) ){
+ count = 0;
+ count1 = 0;
+ numNeighboursSum = 0;
+ for( i=0; i<pd->numInputSites; i++ ){
+ if( pd->attributes->CalculateVoronoiSurfaceArea ){
+ if( pd->processor[i] == MASTER_PROC ){
+ memcpy( &(pd->voronoiArea[i]), &(pd->localTriangulation->voronoiArea[count++]), sizeof( float ) );
+ }
+ else{
+ MPI_Recv( &(pd->voronoiArea[i]), 1, MPI_FLOAT, pd->processor[i], VORONOI_AREA_TAG, (*pd->comm), &st );
+ }
+ }
- self->localTriangulation->neighbours[i][count] =
- self->mappingTable[0][(self->localTriangulation->sites[i].id -
- (self->localTriangulation->idOffset+self->numLocalSites))];
+ if( pd->attributes->FindNeighbours ){
+ if( pd->processor[i] == MASTER_PROC ){
+ memcpy( &(pd->numNeighbours[i]), &(pd->localTriangulation->numNeighbours[count1++]), sizeof( int ) );
}
else{
- self->localTriangulation->neighbours[i][count] =
- self->mappingTable[1][(self->localTriangulation->sites[i].id -
- (self->localTriangulation->idOffset+
- self->numLocalSites+
- self->numHaloSites[0]))];
+ MPI_Recv( &(pd->numNeighbours[i]), 1, MPI_INT, pd->processor[i], NUM_NEIGHBOUR_TAG, (*pd->comm), &st );
}
}
+ numNeighboursSum += pd->numNeighbours[i];
}
+
+ 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 );
+ }
+
+ 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 );
+ }
+
+ stride = 0;
+ for( j=0; j<pd->numInputSites; j++ ){
+ if( pd->attributes->FindNeighbours ){
+ pd->neighbours[j] = pd->neighbours[0]+stride;
+ }
+
+ if( pd->attributes->CalculateVoronoiSides ){
+ pd->voronoiSides[j] = pd->voronoiSides[0]+stride;
+ }
+ stride += pd->numNeighbours[j];
+ }
+
+ count = 0;
+ for(j=0; j<pd->numInputSites; j++){
+ if( pd->processor[j] == MASTER_PROC ){
+ if( pd->attributes->CalculateVoronoiSides ){
+ memcpy( (pd->voronoiSides[j]), (pd->localTriangulation->voronoiSides[count]), sizeof(float)*pd->numNeighbours[j] );
+ }
+ if( pd->attributes->FindNeighbours ){
+ memcpy( (pd->neighbours[j]), (pd->localTriangulation->neighbours[count]), sizeof(int)*pd->numNeighbours[j] );
+ }
+ count++;
+ }
+ else{
+ if( pd->attributes->CalculateVoronoiSides ){
+ MPI_Recv( (pd->voronoiSides[j]), pd->numNeighbours[j], MPI_FLOAT, pd->processor[j], VORONOI_SIDES_TAG, *(pd->comm), &st );
+ }
+ if( pd->attributes->FindNeighbours ){
+ MPI_Recv( (pd->neighbours[j]), pd->numNeighbours[j], MPI_INT, pd->processor[j], NEIGHBOURS_TAG, *(pd->comm), &st );
+ }
+ }
+ }
}
-}
+ else{
+ for( i=0; i<pd->localTriangulation->numInputSites; i++ ){
+ if( pd->attributes->CalculateVoronoiSurfaceArea ){
+ MPI_Send( &(pd->localTriangulation->voronoiArea[i]), 1, MPI_FLOAT, MASTER_PROC, VORONOI_AREA_TAG, (*pd->comm) );
+ }
+
+ if( pd->attributes->FindNeighbours ){
+ MPI_Send( &(pd->localTriangulation->numNeighbours[i]), 1, MPI_INT, MASTER_PROC, NUM_NEIGHBOUR_TAG, (*pd->comm) );
+ }
+ }
+
+ for( i=0; i<pd->localTriangulation->numInputSites; i++ ){
+ if( pd->attributes->CalculateVoronoiSides ){
+ MPI_Send( (pd->localTriangulation->voronoiSides[i]), pd->localTriangulation->numNeighbours[i], MPI_FLOAT, MASTER_PROC, VORONOI_SIDES_TAG, *(pd->comm));
+ }
+ if( pd->attributes->FindNeighbours ){
+ memset( translationArray, 0, sizeof( translationArray ) );
+
+ for( j=0; j<pd->localTriangulation->numNeighbours[i]; j++ ){
+ translationArray[j] = ParallelDelaunay_TranslateLocalToGlobal( pd, pd->localTriangulation->neighbours[i][j] );
+ }
+
+ MPI_Send( translationArray, pd->localTriangulation->numNeighbours[i], MPI_INT, MASTER_PROC, NEIGHBOURS_TAG, *(pd->comm) );
+ }
+ }
+ }
-void _ParallelDelaunay_Initialise( void* pd, void* data )
-{
+ if( pd->attributes->BuildTriangleIndices ){
+ int **triIndices = NULL;
+ int globalNumTriangles = 0;
+ Delaunay *delaunay = NULL;
+ int *triCountArray = NULL;
+ int triCount = 0;
+ int rank = 0;
+
+ delaunay = pd->localTriangulation;
+
+ assert( delaunay );
-}
+ triIndices = delaunay->triangleIndices;
+ triCount = delaunay->numTriangles;
-void _ParallelDelaunay_Execute( void* pd, void* data )
-{
+ MPI_Allreduce( &triCount, &globalNumTriangles, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
+
+ triCountArray = Memory_Alloc_Array_Unnamed( int, sizeof(int)*pd->numProcs );
+ rank = pd->rank;
-}
+ if( rank == MASTER_PROC ){
-void _ParallelDelaunay_Destroy( void* pd, void* data )
-{
-
+ pd->triangleIndices = Memory_Alloc_Array_Unnamed( int*, sizeof(int*) * globalNumTriangles );
+ pd->triangleIndices[0] = Memory_Alloc_Array_Unnamed( int, sizeof(int) * globalNumTriangles*3 );
+
+ triCountArray[0] = triCount;
+ for( i=1; i<pd->numProcs; i++ ){
+ MPI_Recv( &(triCountArray[i]), 1, MPI_INT, i, DATA_TAG, MPI_COMM_WORLD, &st );
+ }
+
+ for( i=0; i<globalNumTriangles; i++ ){
+ pd->triangleIndices[i] = pd->triangleIndices[0] + i*3;
+ }
+
+ for(j=0; j<triCountArray[0]; j++){
+ memcpy( pd->triangleIndices[j], triIndices[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 );
+ }
+ }
+
+ /*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] );
+ }*/
+ }
+ 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 );
+ }
+ }
+
+ Memory_Free( triCountArray );
+ }
}
/*--------------------------------------------------------------------------------------------------------------------------
** Private Member functions
*/
-void ParallelDelaunay_BuildTriangleIndices( Delaunay *delaunay )
+int ParallelDelaunay_TranslateLocalToGlobal( ParallelDelaunay *self, int id )
{
+ if( id >= (self->localTriangulation->idOffset+self->numLocalSites) ){
+
+ if( id < (self->localTriangulation->idOffset+self->numLocalSites + self->numHaloSites[0]) ){
+
+ id = self->mappingTable[0][(id -
+ (self->localTriangulation->idOffset+self->numLocalSites))];
+ }
+ else{
+ id = self->mappingTable[1][(id -
+ (self->localTriangulation->idOffset+
+ self->numLocalSites+
+ self->numHaloSites[0]))];
+ }
+ return id;
+ }
+ else{
+ return id;
+ }
}
#define ORG_TAG 101
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.h 2006-07-21 03:03:26 UTC (rev 4043)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.h 2006-07-21 03:03:38 UTC (rev 4044)
@@ -147,10 +147,16 @@
/*--------------------------------------------------------------------------------------------------------------------------
** Private Member functions
*/
- void ParallelDelaunay_BuildTriangleIndices( Delaunay *delaunay );
+ void ParallelDelaunay_BuildTriangleIndices( ParallelDelaunay *pd );
void ParallelDelaunaySendEdge( QuadEdgeRef edge, int rank, MPI_Comm *comm, MPI_Request *req );
QuadEdgeRef ParallelDelaunayRecvEdge( Delaunay *d, int rank, MPI_Comm *comm );
QuadEdgeRef ParallelDelaunayFindLowestQuadEdge( ParallelDelaunay *pd, MPI_Comm *comm, int rank );
void ParallelDelaunayMerge( ParallelDelaunay *pd, MPI_Comm *comm, int rank );
+ int ParallelDelaunay_TranslateLocalToGlobal( ParallelDelaunay *self, int id );
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Public Member functions
+ */
+ void ParallelDelaunay_GatherTriangulation( ParallelDelaunay *pd );
+
#endif /* __Discretisation_Geometry_ParallelDelaunay_h__ */
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/QuadEdge.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/QuadEdge.c 2006-07-21 03:03:26 UTC (rev 4043)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/QuadEdge.c 2006-07-21 03:03:38 UTC (rev 4044)
@@ -112,7 +112,14 @@
QuadEdge *QuadEdge_New( QuadEdgePool *qp )
{
QuadEdge *e = NULL;
-
+ int index = 0;
+
+ index = qp->numQuadEdgesFree - 1;
+
+ if( index < 0 ){
+ return NULL;
+ }
+
e = qp->pool[--(qp->numQuadEdgesFree)];
SET_IN_USE( (QuadEdgeRef)e );
@@ -172,9 +179,15 @@
QuadEdgeRef MakeQuadEdge( QuadEdgePool *qp )
{
- QuadEdgeRef e;
+ QuadEdgeRef e = 0;
e = (QuadEdgeRef) QuadEdge_New( qp );
+
+ if( e == 0 ){
+ fprintf( stderr, "Out of memory..!\n Aborting..!\n" );
+ assert( 0 );
+ }
+
ONEXT(e) = e;
SYMDNEXT(e) = SYM(e);
ROTRNEXT(e) = TOR(e);
More information about the cig-commits
mailing list