[cig-commits] r3992 - in long/3D/Gale/trunk/src/StGermain: .
Discretisation/Geometry/src
walter at geodynamics.org
walter at geodynamics.org
Sun Jul 9 14:36:25 PDT 2006
Author: walter
Date: 2006-07-09 14:36:24 -0700 (Sun, 09 Jul 2006)
New Revision: 3992
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/Geometry.h
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.c
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.h
Log:
r2496 at earth: boo | 2006-07-09 14:30:56 -0700
r2491 at earth (orig r3671): RaquibulHassan | 2006-07-07 00:42:27 -0700
Committing a working and tested version of the parallelDelaunay triangulator class.
Yet to write a formal test for the class.
Property changes on: long/3D/Gale/trunk/src/StGermain
___________________________________________________________________
Name: svk:merge
- 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2495
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3670
+ 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2496
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3671
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Delaunay.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Delaunay.c 2006-07-09 21:36:21 UTC (rev 3991)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Delaunay.c 2006-07-09 21:36:24 UTC (rev 3992)
@@ -165,17 +165,19 @@
self->numInputSites = numSites;
self->idOffset = idOffset;
- self->sites = Memory_Alloc_Array_Unnamed( Site, sizeof( Site ) * self->numSites );
- memset( self->boundingTriangle, 0, sizeof( self->boundingTriangle ) );
+ if( sites != NULL ){
+ self->sites = Memory_Alloc_Array_Unnamed( Site, sizeof( Site ) * self->numSites );
+ memset( self->boundingTriangle, 0, sizeof( self->boundingTriangle ) );
- for( i=0; i<self->numSites; i++ ){
- if( i < self->numInputSites ){
- self->sites[i].coord = &(sites[i]);
+ for( i=0; i<self->numSites; i++ ){
+ if( i < self->numInputSites ){
+ self->sites[i].coord = &(sites[i]);
+ }
+ else{
+ self->sites[i].coord = &(self->boundingTriangle[i%3]);
+ }
+ self->sites[i].id = i + self->idOffset;
}
- else{
- self->sites[i].coord = &(self->boundingTriangle[i%3]);
- }
- self->sites[i].id = i + self->idOffset;
}
_Stg_Class_Init( (Stg_Class*)self );
@@ -225,20 +227,21 @@
self->numInputSites = numSites;
self->idOffset = idOffset;
- self->sites = Memory_Alloc_Array_Unnamed( Site, sizeof( Site ) * self->numSites );
- memset( self->boundingTriangle, 0, sizeof( self->boundingTriangle ) );
+ if( sites != NULL ){
+ self->sites = Memory_Alloc_Array_Unnamed( Site, sizeof( Site ) * self->numSites );
+ memset( self->boundingTriangle, 0, sizeof( self->boundingTriangle ) );
- for( i=0; i<self->numSites; i++ ){
- if( i < self->numInputSites ){
- self->sites[i].coord = &(sites[i]);
+ for( i=0; i<self->numSites; i++ ){
+ if( i < self->numInputSites ){
+ self->sites[i].coord = &(sites[i]);
+ }
+ else{
+ self->sites[i].coord = &(self->boundingTriangle[i%3]);
+ }
+ self->sites[i].id = i + idOffset;
}
- else{
- self->sites[i].coord = &(self->boundingTriangle[i%3]);
- }
- self->sites[i].id = i + idOffset;
}
-
if( initFlag ){
_Delaunay_Init( self );
}
@@ -256,6 +259,9 @@
*minX = INFINITY;
*minY = INFINITY;
+
+ if (sites == NULL ) return;
+
for( i=0; i<count; i++ ){
if( *maxX < (*(sites[i].coord))[0] ){
*maxX = (*(sites[i].coord))[0];
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Geometry.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Geometry.h 2006-07-09 21:36:21 UTC (rev 3991)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Geometry.h 2006-07-09 21:36:24 UTC (rev 3992)
@@ -66,6 +66,7 @@
#include "IrregGeometry.h"
#include "QuadEdge.h"
#include "Delaunay.h"
+ #include "ParallelDelaunay.h"
#include "Init.h"
#include "Finalise.h"
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-09 21:36:21 UTC (rev 3991)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.c 2006-07-09 21:36:24 UTC (rev 3992)
@@ -58,6 +58,7 @@
#include <assert.h>
#include <string.h>
#include <math.h>
+#include "visualizer.c"
const Type ParallelDelaunay_Type="ParallelDelaunay";
#define PI 3.1415926535897932384626
@@ -86,6 +87,9 @@
NULL,
NULL,
0,
+ 0,
+ 0,
+ NULL,
NULL );
return d;
@@ -96,6 +100,9 @@
Dictionary* dictionary,
CoordF *sites,
int numSites,
+ int rank,
+ int numProcs,
+ MPI_Comm *comm,
DelaunayAttributes *attr )
{
ParallelDelaunay *d = _ParallelDelaunay_New(
@@ -115,11 +122,16 @@
dictionary,
sites,
numSites,
+ rank,
+ numProcs,
+ comm,
attr );
return d;
}
+#define MASTER_PROC 0
+
/** Initialise a ParallelDelaunay */
void ParallelDelaunay_Init(
ParallelDelaunay* self,
@@ -127,6 +139,9 @@
Dictionary* dictionary,
CoordF *sites,
int numSites,
+ int rank,
+ int numProcs,
+ MPI_Comm *comm,
DelaunayAttributes *attr )
{
self->type = ParallelDelaunay_Type;
@@ -161,6 +176,9 @@
self->haloSites[0] = NULL;
self->haloSites[1] = NULL;
self->localTriangulation = NULL;
+ self->rank = rank;
+ self->numProcs = numProcs;
+ self->comm = comm;
_Stg_Class_Init( (Stg_Class*)self );
_Stg_Object_Init( (Stg_Object*)self, name, NON_GLOBAL );
@@ -186,6 +204,9 @@
Dictionary *dictionary,
CoordF *sites,
int numSites,
+ int rank,
+ int numProcs,
+ MPI_Comm *comm,
DelaunayAttributes *attr )
{
ParallelDelaunay *self = NULL;
@@ -200,6 +221,9 @@
self->haloSites[0] = NULL;
self->haloSites[1] = NULL;
self->localTriangulation = NULL;
+ self->rank = rank;
+ self->numProcs = numProcs;
+ self->comm = comm;
if( initFlag ){
_ParallelDelaunay_Init( self );
@@ -208,31 +232,305 @@
return self;
}
+int ParallelDelaunayBtreeCompareFunction( void *a, void *b )
+{
+ Site *s1, *s2;
+
+ s1 = (Site*)a;
+ s2 = (Site*)b;
+
+ if( (*(s1->coord))[0] > (*(s2->coord))[0] ){
+ return 1;
+ }
+ else if( (*(s1->coord))[0] == (*(s2->coord))[0] ){
+ if( (*(s1->coord))[1] > (*(s2->coord))[1] ){
+ return 1;
+ }
+ else if( (*(s1->coord))[1] < (*(s2->coord))[1] ){
+ return -1;
+ }
+ else{
+ return 0;
+ }
+ }
+ else{
+ return -1;
+ }
+}
+
+#define epsilon 0.0001
+#define LOAD_TAG 1
+#define DATA_TAG 1<<1
void _ParallelDelaunay_Init( ParallelDelaunay* self )
{
- float maxX, minX, maxY, minY;
- float centreX, centreY;
- float radius;
+ float _minX, _maxX;
+ float _minY, _maxY;
+ int numProcs, numSites, i, j, count;
+ float stride, start;
+ int *alloced = NULL;
+ int offset;
assert( self );
+ numProcs = self->numProcs;
+ numSites = self->numSites;
+
+ self->numHaloSites[0] = 0;
+ self->numHaloSites[1] = 0;
+
+ if( numProcs == 1 ){
+ self->leftProc = numProcs;
+ self->rightProc = numProcs;
+ }
+ else{
+ if( self->rank == MASTER_PROC ){
+ self->leftProc = numProcs;
+ self->rightProc = self->rank + 1;
+ }
+ else if( self->rank == (numProcs-1) ){
+ self->leftProc = self->rank - 1;
+ self->rightProc = numProcs;
+ }
+ else{
+ self->leftProc = self->rank - 1;
+ self->rightProc = self->rank + 1;
+ }
+ }
- centreX = 0; centreY = 0;
+ self->mapGlobalToLocal = Memory_Alloc_Array_Unnamed( int, sizeof(int)*numSites );
+ self->processorLoad = Memory_Alloc_Array_Unnamed( int, sizeof(int)*numProcs );
+ memset( self->processorLoad, 0, sizeof( int )*numProcs );
+
+ if( self->rank == MASTER_PROC ){
+ self->processor = Memory_Alloc_Array_Unnamed( int, sizeof(int)*numSites );
+
+ alloced = Memory_Alloc_Array_Unnamed( int, sizeof(int)*numSites );
+ memset( alloced, 0, sizeof( int )*numSites );
- Delaunay_FindMinMax( ((Delaunay*)self)->sites, self->numSites, &minX, &minY, &maxX, &maxY );
+
+ Delaunay_FindMinMax( self->sites, self->numSites, &_minX, &_minY, &_maxX, &_maxY );
+ Delaunay_SortSites( self->sites, self->numSites );
- radius = (sqrt((maxX - minX) * (maxX - minX) + (maxY - minY) * (maxY - minY)));
+ stride = (_maxX - _minX)/((float)numProcs);
+
+ start = _minX;
+ for( i=0; i<numProcs; i++ ){
+ for( j=0; j<numSites; j++ ){
+ if( ((*(self->sites[j].coord))[0] >= start-epsilon) &&
+ ((*(self->sites[j].coord))[0] <= (start+stride+epsilon)) &&
+ (!alloced[j]) ){
+
+ alloced[j] = 1;
+ self->processorLoad[i]++;
+ self->processor[j] = i;
+ }
+ }
+ start+=stride;
+ }
+
+ for( i=0; i<numProcs; i++ ){
+ printf( "processorLoad[%d] = %d\n", i, self->processorLoad[i] );
+ }
+
+ for( i=MASTER_PROC+1; i<numProcs; i++ ){
+ MPI_Send( &(self->processorLoad[i]), 1, MPI_INT, i, LOAD_TAG, *self->comm );
+ }
+
+ self->numLocalSites = self->processorLoad[MASTER_PROC];
+ self->localPoints = Memory_Alloc_Array_Unnamed( CoordF, sizeof( CoordF )*self->numLocalSites );
+
+ count = 0;
+ for( i=0; i<numSites; i++ ){
+ if( self->processor[i] == MASTER_PROC ){
+ memcpy( &(self->localPoints[count++]), self->sites[i].coord, sizeof(CoordF) );
+ }
+ else{
+ MPI_Send( self->sites[i].coord, sizeof( CoordF ), MPI_BYTE, self->processor[i], DATA_TAG, *self->comm );
+ }
+ }
+
+ Memory_Free( alloced );
+
+ {/*TODO remove this block*/
- centreX = minX + (maxX - minX) / 2.0f;
- centreY = minY + (maxY - minY) / 2.0f;
+ maxX = _maxX;
+ maxY = _maxY;
+ minX = _minX;
+ minY = _minY;
+ for( i=MASTER_PROC+1; i<numProcs; i++ ){
+ MPI_Send( &minX, 1, MPI_FLOAT, i, LOAD_TAG, *self->comm );
+ MPI_Send( &minY, 1, MPI_FLOAT, i, LOAD_TAG, *self->comm );
+ MPI_Send( &maxX, 1, MPI_FLOAT, i, LOAD_TAG, *self->comm );
+ MPI_Send( &maxY, 1, MPI_FLOAT, i, LOAD_TAG, *self->comm );
+ }
+ }
+ }
+ else{
+ MPI_Status status;
+
+ MPI_Recv( &self->numLocalSites, 1, MPI_INT, MASTER_PROC, LOAD_TAG, *self->comm, &status );
- self->boundingTriangle[0][0] = centreX - tan(PI/3.0f)*radius;
- self->boundingTriangle[0][1] = centreY - radius;
+ self->localPoints = Memory_Alloc_Array_Unnamed( CoordF, sizeof( CoordF )*self->numLocalSites );
- self->boundingTriangle[1][0] = centreX + tan(PI/3.0f)*radius;
- self->boundingTriangle[1][1] = centreY - radius;
+ for( i=0; i<self->numLocalSites; i++ ){
+ MPI_Recv( &(self->localPoints[i]), sizeof(CoordF), MPI_BYTE, MASTER_PROC, DATA_TAG, *self->comm, &status );
+ }
+
+ {/*TODO remove this block*/
+ MPI_Recv( &minX, 1, MPI_FLOAT, MASTER_PROC, LOAD_TAG, *self->comm, &status );
+ MPI_Recv( &minY, 1, MPI_FLOAT, MASTER_PROC, LOAD_TAG, *self->comm, &status );
+ MPI_Recv( &maxX, 1, MPI_FLOAT, MASTER_PROC, LOAD_TAG, *self->comm, &status );
+ MPI_Recv( &maxY, 1, MPI_FLOAT, MASTER_PROC, LOAD_TAG, *self->comm, &status );
+ }
+ }
+
+ MPI_Bcast( self->processorLoad, numProcs, MPI_INT, MASTER_PROC, *self->comm );
+ self->numTotalLocalSites = self->numLocalSites;
+
+ offset = 0;
+ for( i=0; i<self->rank; i++ ){
+ offset += self->processorLoad[i];
+ }
+
+ self->localTriangulation = Delaunay_New( "delaunay", self->dictionary, self->localPoints, self->numLocalSites, offset, self->attributes );
+ self->localTriangulation->qp = QuadEdgePool_New( self->localTriangulation->numSites * 4 );
+ Delaunay_SortSites(self->localTriangulation->sites, self->localTriangulation->numSites);
+ Delaunay_Recurse(self->localTriangulation, 0, self->localTriangulation->numSites,
+ &self->localTriangulation->leftMost, &self->localTriangulation->rightMost);
+
+ for( i=0; i<numSites; i++ ){
+ self->mapGlobalToLocal[i] = numSites;
+ }
+
+ for( i=0; i<self->numLocalSites; i++ ){
+ self->mapGlobalToLocal[self->localTriangulation->sites[i].id] = self->localTriangulation->sites[i].id;
+ }
+
+ if( self->leftProc != self->numProcs ){
+ BTreeIterator *iter = NULL;
+ Site *result = NULL;
+
+ ParallelDelaunayMerge( self, self->comm, self->leftProc );
+
+ self->haloSites[0] = BTree_New( ParallelDelaunayBtreeCompareFunction, NULL, NULL, NULL, BTREE_NO_DUPLICATES );
+
+ for( i=0; i<self->localTriangulation->qp->numQuadEdges; i++ ){
+ if( IS_FREE((QuadEdgeRef)&(self->localTriangulation->qp->quadEdges[i])) ) continue;
+
+ if( self->mapGlobalToLocal[((Site*)self->localTriangulation->qp->quadEdges[i].data[0])->id] == numSites )
+ BTree_InsertNode( self->haloSites[0], ((Site*)self->localTriangulation->qp->quadEdges[i].data[0]), sizeof( Site* ) );
- self->boundingTriangle[2][0] = centreX;
- self->boundingTriangle[2][1] = centreY + radius/cos(PI/3.0f);
+ if( self->mapGlobalToLocal[((Site*)self->localTriangulation->qp->quadEdges[i].data[2])->id] == numSites )
+ BTree_InsertNode( self->haloSites[0], ((Site*)self->localTriangulation->qp->quadEdges[i].data[2]), sizeof( Site* ) );
+ }
+
+ 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 );
+ memset( self->mappingTable[0], 0, sizeof( int ) * self->haloSites[0]->nodeCount );
+
+ count=0;
+ i = self->numTotalLocalSites;
+ iter = BTreeIterator_New( self->haloSites[0] );
+ for( result=(Site*)BTreeIterator_First(iter);
+ result;
+ result=(Site*)BTreeIterator_Next(iter)){
+
+ self->mappingTable[0][count++] = result->id;
+ memcpy( &(self->localPoints[i++]), result->coord, sizeof( CoordF ) );
+ }
+ Stg_Class_Delete( self->localTriangulation );
+
+ self->localTriangulation = Delaunay_New( "delaunay", self->dictionary, self->localPoints, self->numLocalSites, offset, self->attributes );
+ self->localTriangulation->qp = QuadEdgePool_New( self->localTriangulation->numSites * 4 );
+ Delaunay_SortSites(self->localTriangulation->sites, self->localTriangulation->numSites);
+ Delaunay_Recurse(self->localTriangulation, 0, self->localTriangulation->numSites,
+ &self->localTriangulation->leftMost, &self->localTriangulation->rightMost);
+
+ self->numHaloSites[0] = self->haloSites[0]->nodeCount;
+ self->numTotalLocalSites += self->haloSites[0]->nodeCount;
+ b[0] = 0;
+ }
+
+ if( self->rightProc != self->numProcs ){
+ BTreeIterator *iter = NULL;
+ Site *result = NULL;
+
+ ParallelDelaunayMerge( self, self->comm, self->rightProc );
+
+ self->haloSites[1] = BTree_New( ParallelDelaunayBtreeCompareFunction, NULL, NULL, NULL, BTREE_NO_DUPLICATES );
+
+ for( i=0; i<self->localTriangulation->qp->numQuadEdges; i++ ){
+ if( IS_FREE((QuadEdgeRef)&(self->localTriangulation->qp->quadEdges[i])) ) continue;
+
+ if( self->mapGlobalToLocal[((Site*)self->localTriangulation->qp->quadEdges[i].data[0])->id] == numSites )
+ BTree_InsertNode( self->haloSites[1], ((Site*)self->localTriangulation->qp->quadEdges[i].data[0]), sizeof( Site* ) );
+
+ if( self->mapGlobalToLocal[((Site*)self->localTriangulation->qp->quadEdges[i].data[2])->id] == numSites )
+ BTree_InsertNode( self->haloSites[1], ((Site*)self->localTriangulation->qp->quadEdges[i].data[2]), sizeof( Site* ) );
+ }
+
+ 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 );
+ memset( self->mappingTable[1], 0, sizeof( int ) * self->haloSites[1]->nodeCount );
+
+ count=0;
+ i = self->numTotalLocalSites;
+ iter = BTreeIterator_New( self->haloSites[1] );
+ for( result=(Site*)BTreeIterator_First(iter);
+ result;
+ result=(Site*)BTreeIterator_Next(iter)){
+
+ self->mappingTable[1][count++] = result->id;
+ memcpy( &(self->localPoints[i++]), result->coord, sizeof( CoordF ) );
+ }
+ Stg_Class_Delete( self->localTriangulation );
+
+ self->localTriangulation = Delaunay_New( "delaunay", self->dictionary, self->localPoints, self->numLocalSites, offset, self->attributes );
+ self->localTriangulation->qp = QuadEdgePool_New( self->localTriangulation->numSites * 4 );
+ Delaunay_SortSites(self->localTriangulation->sites, self->localTriangulation->numSites);
+ Delaunay_Recurse(self->localTriangulation, 0, self->localTriangulation->numSites,
+ &self->localTriangulation->leftMost, &self->localTriangulation->rightMost);
+
+ self->numHaloSites[1] = self->haloSites[1]->nodeCount;
+ self->numTotalLocalSites += self->haloSites[1]->nodeCount;
+ b[1] = 0;
+ }
+
+ Stg_Class_Delete( self->localTriangulation );
+ 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->localTriangulation->sites[i].id <
+ (self->localTriangulation->idOffset+self->numLocalSites + self->numHaloSites[0]) ){
+
+ self->localTriangulation->sites[i].id =
+ self->mappingTable[0][(self->localTriangulation->sites[i].id -
+ (self->localTriangulation->idOffset+self->numLocalSites))];
+ }
+ else{
+ self->localTriangulation->sites[i].id =
+ self->mappingTable[1][(self->localTriangulation->sites[i].id -
+ (self->localTriangulation->idOffset+
+ self->numLocalSites+
+ self->numHaloSites[0]))];
+ }
+ }
+ }
+
+ {/*TODO remove this block*/
+
+ createWindow();
+ LoadGLTexture();
+ BuildFont();
+
+ windowEvents( self->localTriangulation, self->localPoints );
+ }
}
/*--------------------------------------------------------------------------------------------------------------------------
@@ -300,3 +598,339 @@
}
+#define ORG_TAG 101
+#define DEST_TAG 102
+#define BREAK_TAG 103
+typedef struct SitePacket_t{
+ int id;
+ float xyz[3];
+}SitePacket;
+void ParallelDelaunaySendEdge( QuadEdgeRef edge, int rank, MPI_Comm *comm, MPI_Request *req )
+{
+ SitePacket sp;
+ assert(edge);
+
+ memcpy( sp.xyz, ((Site*)ORG(edge))->coord, sizeof( CoordF ) );
+ sp.id = ((Site*)ORG(edge))->id;
+ MPI_Isend( &sp, sizeof(SitePacket), MPI_BYTE, rank, ORG_TAG, *comm, &(req[0]) );
+
+ memcpy( sp.xyz, ((Site*)DEST(edge))->coord, sizeof( CoordF ) );
+ sp.id = ((Site*)DEST(edge))->id;
+ MPI_Isend( &sp, sizeof(SitePacket), MPI_BYTE, rank, DEST_TAG, *comm, &(req[1]) );
+}
+
+QuadEdgeRef ParallelDelaunayRecvEdge( Delaunay *d, int rank, MPI_Comm *comm )
+{
+ QuadEdgeRef edge = 0;
+ CoordF *c;
+ Site *s;
+ MPI_Status st;
+ SitePacket sp[2];
+
+ assert( d );
+
+ edge = MakeQuadEdge( d->qp );
+
+ assert(edge);
+
+ s = Memory_Alloc_Array_Unnamed( Site, sizeof( Site ) * 2 );
+ memset( s, 0, sizeof( Site ) * 2 );
+
+ c = Memory_Alloc_Array_Unnamed( CoordF, sizeof( CoordF ) * 2 );
+ memset( c, 0, sizeof( CoordF ) * 2 );
+
+ MPI_Recv( &(sp[0]), sizeof(SitePacket), MPI_BYTE, rank, ORG_TAG, *comm, &st );
+ MPI_Recv( &(sp[1]), sizeof(SitePacket), MPI_BYTE, rank, DEST_TAG, *comm, &st );
+
+ memcpy( &(c[0]), sp[0].xyz, sizeof( CoordF ) );
+ memcpy( &(c[1]), sp[1].xyz, sizeof( CoordF ) );
+
+ s[0].id = sp[0].id;
+ s[0].coord = &(c[0]);
+ s[1].id = sp[1].id;
+ s[1].coord = &(c[1]);
+
+ ORG(edge)=&(s[0]);
+ DEST(edge)=&(s[1]);
+
+ return edge;
+}
+
+int ParallelDelaunayListCompareFunction( void *a, void *b )
+{
+ QuadEdgeRef e1, e2;
+
+ e1 = (QuadEdgeRef)a;
+ e2 = (QuadEdgeRef)b;
+
+ if( e1 > e2 ){
+ return 1;
+ }
+ else if( e1 < e2 ){
+ return -1;
+ }
+ else{
+ return 0;
+ }
+}
+
+void ParallelDelaunayListDeleteFunction( void *a )
+{
+
+}
+
+QuadEdgeRef ParallelDelaunayFindLowestQuadEdge( ParallelDelaunay *pd, MPI_Comm *comm, int rank )
+{
+ QuadEdgeRef ldi=0, rdi=0;
+ MPI_Request r[2];
+ MPI_Status s[2];
+ int globalBreak, localBreak;
+ LinkedList *list = NULL;
+ LinkedListIterator *iter;
+ QuadEdgeRef result = 0;
+ Delaunay *d = NULL;
+
+ assert( pd );
+ d = pd->localTriangulation;
+
+ list = LinkedList_New( ParallelDelaunayListCompareFunction, NULL, NULL, ParallelDelaunayListDeleteFunction, LINKEDLIST_UNSORTED );
+ iter = LinkedListIterator_New( list );
+
+ if( rank == pd->leftProc ){
+ rdi = d->leftMost;
+ }
+ else if( rank == pd->rightProc ){
+ ldi = d->rightMost;
+ }
+ else{
+ fprintf( stderr, "Failed to find lowest edge on rank %d..!\n Aborting..!\n", rank );
+ }
+
+ if( rank == pd->numProcs ) return 0;
+
+ localBreak = 0;
+ globalBreak = 0;
+ while (1)
+ {
+ localBreak = 0;
+ globalBreak = 0;
+
+ if( rank == pd->leftProc ){
+ ParallelDelaunaySendEdge( rdi, rank, comm, &(r[0]) );
+ ldi = ParallelDelaunayRecvEdge( d, rank, comm );
+ MPI_Waitall( 2, r, s );
+ LinkedList_InsertNode( list, (void*)ldi, sizeof( QuadEdgeRef* ) );
+
+ if (RightOf(ORG(ldi), rdi)){
+ rdi = ONEXT(SYM(rdi));
+ localBreak = 1;
+ }
+ else{
+ }
+ }
+
+ if( rank == pd->rightProc ){
+ ParallelDelaunaySendEdge( ldi, rank, comm, &(r[0]) );
+ rdi = ParallelDelaunayRecvEdge( d, rank, comm );
+ MPI_Waitall( 2, r, s );
+ LinkedList_InsertNode( list, (void*)rdi, sizeof( QuadEdgeRef* ) );
+
+ if (LeftOf(ORG(rdi), ldi)){
+ ldi = LNEXT(ldi);
+ localBreak = 1;
+ }
+ else{
+ }
+ }
+
+ MPI_Isend( &localBreak, 1, MPI_INT, rank, BREAK_TAG, *comm, &(r[0]) );
+ MPI_Recv( &(globalBreak), 1, MPI_INT, rank, BREAK_TAG, *comm, &(s[0]) );
+ MPI_Waitall( 1, r, s );
+ globalBreak |= localBreak;
+
+ if( globalBreak == 0 ){
+ break;
+ }
+ }
+
+ for( result=(QuadEdgeRef)((void*)LinkedListIterator_First(iter));
+ result != 0;
+ result=(QuadEdgeRef)((void*)LinkedListIterator_Next(iter)) ){
+ Site *s;
+
+ s = (Site*)ORG(result);
+
+ Memory_Free( s->coord );
+ Memory_Free( s );
+
+ DeleteQuadEdge( d->qp, result );
+ }
+
+ Stg_Class_Delete( list );
+ Stg_Class_Delete( iter );
+
+ if( rank == pd->leftProc ){
+ return rdi;
+ }
+ else if( rank == pd->rightProc ){
+ return ldi;
+ }
+
+ return 0;
+}
+
+void ParallelDelaunayMerge( ParallelDelaunay *pd, MPI_Comm *comm, int rank )
+{
+ QuadEdgeRef lowest = 0, lcand = 0, rcand = 0, basel = 0, baselPrev;
+ MPI_Request r[2];
+ MPI_Status s[2];
+ int localBreak=0, globalBreak=0;
+ double result = 0.0f;
+ Delaunay *d = NULL;
+
+ assert( pd );
+
+ d = pd->localTriangulation;
+
+ lowest = ParallelDelaunayFindLowestQuadEdge( pd, comm, rank );
+
+ localBreak = 0;
+ globalBreak = 0;
+
+ if( rank == pd->numProcs ){
+ fprintf( stderr, "Failed to merge rank %d with rank %d..!\nAborting..!\n", pd->rank, rank );
+ assert( 0 );
+ }
+
+ if( rank == pd->leftProc ){
+ b[0] = rcand = lowest;
+
+ ParallelDelaunaySendEdge( lowest, rank, comm, &(r[0]) );
+ lcand = ParallelDelaunayRecvEdge( d, rank, comm );
+
+ MPI_Waitall( 2, r, s );
+
+ basel = MakeQuadEdge( d->qp );
+ ORG(basel) = DEST(SYM(rcand));
+ DEST(basel) = ORG(lcand);
+
+ DeleteQuadEdge( d->qp, lcand );
+
+ CCW( ((Site*)DEST(rcand))->coord, ((Site*)DEST(basel))->coord, ((Site*)ORG(basel))->coord, &result);
+ if( result == 0.0f ){
+ rcand = OPREV(rcand);
+ }
+
+ while(1){
+ localBreak = 0;
+ if (RightOf(DEST(rcand), basel)){
+ while (InCircle(DEST(basel), ORG(basel), DEST(rcand), DEST(OPREV(rcand)))){
+ QuadEdgeRef t = OPREV(rcand);
+
+ DeleteQuadEdge(d->qp, rcand);
+ rcand = t;
+ }
+ }
+
+ if (!RightOf(DEST(rcand), basel)) localBreak = 1;
+
+ MPI_Isend( &localBreak, 1, MPI_INT, rank, BREAK_TAG, *comm, &(r[0]) );
+ MPI_Recv( &(globalBreak), 1, MPI_INT, rank, BREAK_TAG, *comm, &(s[0]) );
+ MPI_Waitall( 1, r, s );
+ globalBreak &= localBreak;
+
+ if( globalBreak ){
+ break;
+ }
+
+ ParallelDelaunaySendEdge( rcand, rank, comm, &(r[0]) );
+ lcand = ParallelDelaunayRecvEdge( d, rank, comm );
+ MPI_Waitall( 2, r, s );
+
+ if ( !RightOf(DEST(lcand), basel) ||
+ ( RightOf(DEST(rcand), basel) &&
+ InCircle(DEST(lcand), ORG(lcand), ORG(rcand), DEST(rcand)))){
+
+ baselPrev = basel;
+ basel = MakeQuadEdge(d->qp);
+ ORG(basel) = DEST(rcand);
+ DEST(basel) = ORG(SYM(baselPrev));
+
+ rcand=LNEXT(rcand);
+ DeleteQuadEdge( d->qp, lcand );
+ }
+ else{
+ baselPrev = basel;
+ basel = MakeQuadEdge(d->qp);
+ ORG(basel) = DEST(SYM(baselPrev));
+ DEST(basel) = ORG(SYM(lcand));
+ }
+ }
+ }
+ else if( rank == pd->rightProc ){
+ b[1] = lcand = lowest;
+
+ ParallelDelaunaySendEdge( lowest, rank, comm, &(r[0]) );
+ rcand = ParallelDelaunayRecvEdge( d, rank, comm );
+
+ MPI_Waitall( 2, r, s );
+
+ basel = MakeQuadEdge(d->qp);
+ ORG(basel) = DEST(SYM(rcand));
+ DEST(basel) = ORG(lcand);
+
+ DeleteQuadEdge( d->qp, rcand );
+
+ CCW( ((Site*)DEST(lcand))->coord, ((Site*)DEST(basel))->coord, ((Site*)ORG(basel))->coord, &result);
+ if( result == 0.0f ){
+ lcand = ONEXT(lcand);
+ }
+
+ while(1){
+ localBreak = 0;
+ if (RightOf(DEST(lcand), basel)){
+ while (InCircle(DEST(basel), ORG(basel), DEST(lcand), DEST(ONEXT(lcand)))){
+ QuadEdgeRef t = ONEXT(lcand);
+
+ DeleteQuadEdge(d->qp, lcand);
+ lcand = t;
+ }
+ }
+
+ if (!RightOf(DEST(lcand), basel)) localBreak = 1;
+
+ MPI_Isend( &localBreak, 1, MPI_INT, rank, BREAK_TAG, *comm, &(r[0]) );
+ MPI_Recv( &(globalBreak), 1, MPI_INT, rank, BREAK_TAG, *comm, &(s[0]) );
+ MPI_Waitall( 1, r, s );
+ globalBreak &= localBreak;
+
+ if( globalBreak ){
+ break;
+ }
+
+ ParallelDelaunaySendEdge( lcand, rank, comm, &(r[0]) );
+ rcand = ParallelDelaunayRecvEdge( d, rank, comm );
+ MPI_Waitall( 2, r, s );
+
+ if ( !RightOf(DEST(lcand), basel) ||
+ ( RightOf(DEST(rcand), basel) &&
+ InCircle(DEST(lcand), ORG(lcand), ORG(rcand), DEST(rcand)))){
+
+ baselPrev = basel;
+ basel = MakeQuadEdge(d->qp);
+ ORG(basel) = DEST(rcand);
+ DEST(basel) = ORG(SYM(baselPrev));
+ }
+ else{
+ baselPrev = basel;
+ basel = MakeQuadEdge(d->qp);
+ ORG(basel) = DEST(SYM(baselPrev));
+ DEST(basel) = ORG(SYM(lcand));
+
+ lcand = RPREV(lcand);
+ DeleteQuadEdge( d->qp, rcand );
+ }
+ }
+ }
+}
+
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-09 21:36:21 UTC (rev 3991)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ParallelDelaunay.h 2006-07-09 21:36:24 UTC (rev 3992)
@@ -49,11 +49,21 @@
#define __ParallelDelaunay \
__Delaunay \
CoordF *points; \
+ CoordF *localPoints; \
int leftProc; \
int rightProc; \
int numProcs; \
BTree *haloSites[2]; \
- Delaunay *localTriangulation;
+ int numHaloSites[2]; \
+ Delaunay *localTriangulation; \
+ int *mappingTable[2]; \
+ int *mapGlobalToLocal; \
+ int *processor; \
+ int *processorLoad; \
+ int rank; \
+ int numLocalSites; \
+ int numTotalLocalSites; \
+ MPI_Comm *comm;
struct ParallelDelaunay { __ParallelDelaunay };
@@ -69,6 +79,9 @@
Dictionary* dictionary,
CoordF *sites,
int numSites,
+ int rank,
+ int numProcs,
+ MPI_Comm *comm,
DelaunayAttributes *attr );
/** Initialise a ParallelDelaunay */
@@ -78,6 +91,9 @@
Dictionary* dictionary,
CoordF *sites,
int numSites,
+ int rank,
+ int numProcs,
+ MPI_Comm *comm,
DelaunayAttributes *attr );
/** Creation implementation */
@@ -98,6 +114,9 @@
Dictionary *dictionary,
CoordF *sites,
int numSites,
+ int rank,
+ int numProcs,
+ MPI_Comm *comm,
DelaunayAttributes *attr );
void _ParallelDelaunay_Init( ParallelDelaunay* self );
@@ -129,5 +148,9 @@
** Private Member functions
*/
void ParallelDelaunay_BuildTriangleIndices( Delaunay *delaunay );
+ 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 );
#endif /* __Discretisation_Geometry_ParallelDelaunay_h__ */
More information about the cig-commits
mailing list