SPModel: libSPModel/src/SurfaceMesh.c Source File
VPAC - Computational Software Development
Main | SPModel | StGermain FrameWork |
Main Page | Alphabetical List | Class List | Directories | File List | Class Members | File Members

SurfaceMesh.c

Go to the documentation of this file.
00001 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00002 **
00003 ** Copyright (C), 2004, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street, Melbourne, 3053, Australia.
00004 **
00005 ** Authors:
00006 **  Ogar R. Widjaja, Computational Scientist, VPAC.
00007 **  Raquibul Hassan, Software Engineer, VPAC. (raq@vpac.org)
00008 **  Keith Hsuan, Computational Scientist, VPAC (keith@vpac.org)
00009 **  William F. Appelbe, Director, VPAC. (bill@vpac.org)
00010 **  Stevan M. Quenette, Senior Software Engineer, VPAC. (steve@vpac.org)
00011 **  Patrick D. Sunter, Software Engineer, VPAC. (patrick@vpac.org)
00012 **
00013 ** This file may be distributed under the terms of the VPAC Public License
00014 ** as defined by VPAC of Australia and appearing in the file
00015 ** LICENSE.VPL included in the packaging of this file.
00016 **
00017 ** This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
00018 ** WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00019 **
00020 ** $Id: SurfaceMesh.c 262 2006-02-16 02:51:37Z RaquibulHassan $
00021 **
00022 **~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
00023 
00024 #include <mpi.h>
00025 #include <StGermain/StGermain.h>
00026 #include <Cascade/cascade.h>
00027 
00028 
00029 #include "types.h"
00030 #include "SurfaceMesh.h"
00031 
00032 #include <stdio.h>
00033 #include <stdlib.h>
00034 #include <string.h>
00035 #include <assert.h>
00036 #include <math.h>
00037 #include "SurfaceMesh.h"
00038 #include "SurfaceMeshLoader.h"
00039 #include "Misc.h"
00040 #include "CommHandler.h"
00041 
00042 /* Textual name of this class */
00043 const Type SurfaceMesh_Type = "SurfaceMesh";
00044 SurfaceMesh *globalMesh = NULL;
00045 SurfaceMesh *localMesh = NULL;
00046 
00047 SurfaceMesh* SurfaceMesh_DefaultNew( Name name )
00048 {
00049     return _SurfaceMesh_New(
00050         sizeof(SurfaceMesh),
00051         SurfaceMesh_Type,
00052         _SurfaceMesh_Delete, 
00053         _SurfaceMesh_Print,
00054         NULL,
00055         (void*)SurfaceMesh_DefaultNew,
00056         _SurfaceMesh_Construct,
00057         _SurfaceMesh_Build,
00058         _SurfaceMesh_Initialise,
00059         _SurfaceMesh_Execute, 
00060         _SurfaceMesh_Destroy, 
00061         name,
00062         False,
00063         0,
00064         0,
00065         NULL,
00066         NULL );
00067 }
00068 
00069 SurfaceMesh* SurfaceMesh_New(
00070         Name                    name,
00071         int                     numProcs,
00072         int                     rank,
00073         void*                   extension_Register,
00074         Dictionary*             dictionary )
00075 {
00076     return _SurfaceMesh_New(
00077         sizeof(SurfaceMesh),
00078         SurfaceMesh_Type,
00079         _SurfaceMesh_Delete, 
00080         _SurfaceMesh_Print,
00081         NULL,
00082         (void*)SurfaceMesh_DefaultNew,
00083         _SurfaceMesh_Construct,
00084         _SurfaceMesh_Build,
00085         _SurfaceMesh_Initialise,
00086         _SurfaceMesh_Execute, 
00087         _SurfaceMesh_Destroy, 
00088         name,
00089         True,
00090         numProcs,
00091         rank,
00092         extension_Register,
00093         dictionary );
00094 }
00095 
00096 SurfaceMesh* _SurfaceMesh_New(
00097         SizeT                   _sizeOfSelf,
00098         Type                    type,
00099         Stg_Class_DeleteFunction*               _delete,
00100         Stg_Class_PrintFunction*                _print,
00101         Stg_Class_CopyFunction*             _copy, 
00102         Stg_Component_DefaultConstructorFunction*   _defaultConstructor,
00103         Stg_Component_ConstructFunction*            _construct,
00104         Stg_Component_BuildFunction*        _build,
00105         Stg_Component_InitialiseFunction*       _initialise,
00106         Stg_Component_ExecuteFunction*      _execute,
00107         Stg_Component_DestroyFunction*      _destroy,
00108         Name                    name,
00109         Bool                            initFlag,
00110         int                     numProcs,
00111         int                     rank,
00112         void*                   extension_Register,
00113         Dictionary*             dictionary )
00114 {
00115     SurfaceMesh*            self;
00116     
00117     /* Allocate memory */
00118     assert( _sizeOfSelf >= sizeof(SurfaceMesh) );
00119     /* Construct using parent */
00120     self = (SurfaceMesh*)_Stg_Component_New( _sizeOfSelf, type, _delete, _print, NULL, 
00121             _defaultConstructor, _construct, _build, _initialise, _execute, _destroy,
00122         name, NON_GLOBAL );
00123 
00124     /* General info */
00125     self->dictionary = dictionary;
00126     self->numProcs = numProcs;
00127     self->rank = rank;
00128     
00129     /* Virtual functions */
00130     
00131     /* SurfaceMesh info */
00132     if( initFlag ){
00133         _SurfaceMesh_Init( self );
00134     }
00135     
00136     return self;
00137 }
00138 
00139 void _SurfaceMesh_Init( SurfaceMesh* self ) 
00140 {
00141     Dictionary* meshDict = NULL;
00142     int nx, ny; 
00143     char *inputType;
00144     
00145     /* General and Virtual info should already be set */
00146     
00147     /* SurfaceMesh info */
00148     self->debug = Journal_Register( Info_Type, "debugStream" );
00149     meshDict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( self->dictionary, "mesh" ) ); 
00150 
00151     if ( !meshDict ) {
00152         inputType = (char*)MESH_TYPE_AUTOGENERATED;
00153         meshDict = self->dictionary;
00154     }
00155     else {
00156         inputType = Dictionary_GetString_WithDefault( meshDict, "type", (char*)MESH_TYPE_AUTOGENERATED );
00157     }
00158     
00159     fprintf(stderr,"INPUT TYPE %s \n", inputType);
00160 
00161     if ( 0 == strcmp( inputType, MESH_TYPE_AUTOGENERATED ) ) {
00162         nx = Dictionary_GetUnsignedInt_WithDefault( meshDict, "nx", 6); /* Number of nodes in x-direction */
00163         ny = Dictionary_GetUnsignedInt_WithDefault( meshDict, "ny", 6); /* Number of nodes in y-direction */
00164 
00165         self->numNodes = nx*ny;
00166     }
00167     else if ( 0 == strcmp( inputType, MESH_TYPE_FROM_TEXT_INPUT ) ) {
00168         
00169         self->numNodes = Dictionary_GetUnsignedInt_WithDefault( meshDict, "nodeCount", 2000 );
00170     }
00171     else if ( 0 == strcmp( inputType, MESH_TYPE_FROM_ENVI_INPUT ) ) {
00172 
00173     }
00174     else if ( 0 == strcmp( inputType, MESH_TYPE_FROM_SPM_INPUT ) ) {
00175         self->numNodes = Dictionary_GetUnsignedInt_WithDefault( meshDict, "nodeCount", 2000 );
00176     }
00177     else if ( 0 == strcmp( inputType, MESH_TYPE_FROM_OUTPUT_FILES ) ) {
00178         
00179         self->numNodes = Dictionary_GetUnsignedInt_WithDefault( meshDict, "nodeCount", 2000 );
00180     }
00181 
00182     self->iAdapt = Dictionary_GetBool_WithDefault( meshDict, "iadapt", False );
00183     self->sideX = Dictionary_GetFloat_WithDefault( meshDict, "sidex", 100000 );
00184     self->sideY = Dictionary_GetFloat_WithDefault( meshDict, "sidey", 100000 );
00185     /* Maximum number of neighbours per node */
00186     self->maxNeighboursPerNode = Dictionary_GetUnsignedInt_WithDefault(self->dictionary, "maxNeighboursPerNode", 20);
00187 
00188     
00189     self->surfaceMeshLoader = SurfaceMeshLoader_New( "surfaceMeshLoader", self->dictionary );
00190     SurfaceMeshLoader_ScanMesh( self->surfaceMeshLoader, self );
00191 }
00192 
00193 void _SurfaceMesh_Delete( void* surfaceMesh ) {
00194     SurfaceMesh* self = (SurfaceMesh*)surfaceMesh;
00195     
00196     Journal_DPrintf( self->debug, "In %s - for %s\n", __func__, self->name );
00197     
00198     Stg_Class_Delete( self->surfaceMeshLoader );
00199     SurfaceMesh_ReleaseMemory( self );
00200     _Stg_Component_Delete( self );
00201 }
00202 
00203 void _SurfaceMesh_Print( void* surfaceMesh, Stream* stream ) {
00204     SurfaceMesh* self = (SurfaceMesh*)surfaceMesh;
00205     int i = 0;
00206     void *data;
00207     
00208     /* General info */
00209     Journal_Printf( stream, "SurfaceMesh (ptr): %p\n", self );
00210     
00211     for( i=0; i<self->numNodes; i++ ){
00212         LinkedListIterator *iter = LinkedListIterator_New( self->nodeProviders[i] );
00213         
00214         Journal_Printf( stream, "node id \t%d\n", self->id[i] );
00215         Journal_Printf( stream, "receiver id \t%d\n", self->receiver[i] );
00216         Journal_Printf( stream, "num providers \t%d\n", self->nodeProviders[i]->nodeCount );
00217         Journal_Printf( stream, "provider ids\n" );
00218         for( data=LinkedListIterator_First( iter ); data!=NULL; data=LinkedListIterator_Next(iter) ){
00219             Journal_Printf( stream, "\t %d\n", *(int*)data );
00220         }
00221         Journal_Printf( stream, "size \t\t%d\n", self->size[i] );
00222         
00223         Stg_Class_Delete( iter );
00224     }
00225     /* Print parent */
00226     
00227     /* Virtual info */
00228     
00229     /* SurfaceMesh info */
00230 }
00231 
00232 void _SurfaceMesh_Construct( void* surfaceMesh, Stg_ComponentFactory* cf )
00233 {
00234     
00235 }
00236 
00237 int heightCompareFunc( void *data1, void *data2 )
00238 {
00239     int *node1, *node2;
00240 
00241     node1 = (int*)data1;
00242     node2 = (int*)data2;
00243 
00244     if( data1 == NULL || data2 == NULL ){
00245         return 0;
00246     }
00247     
00248     if(globalMesh->h[*node1-1] < globalMesh->h[*node2-1]){
00249         return 1;
00250     }
00251     else if (globalMesh->h[*node1-1] == globalMesh->h[*node2-1]){
00252         return 0;
00253     }
00254     else{
00255         return -1;
00256     }
00257 }
00258 
00259 int idCompareFunc( const void *data1, const void *data2 )
00260 {
00261     int *id1, *id2;
00262 
00263     id1 = (int*)data1;
00264     id2 = (int*)data2;
00265 
00266     if( data1 == NULL || data2 == NULL ){
00267         return 0;
00268     }
00269     
00270     if(localMesh->h[localMesh->mapGlobalToLocal[*id1-1]] <
00271             localMesh->h[localMesh->mapGlobalToLocal[*id2-1]]){
00272         return 1;
00273     }
00274     else if (localMesh->h[localMesh->mapGlobalToLocal[*id1-1]] == 
00275             localMesh->h[localMesh->mapGlobalToLocal[*id2-1]]){
00276         return 0;
00277     }
00278     else{
00279         return -1;
00280     }
00281 }
00282 
00283 void nodePrintFunction( void *id, Stream *stream )
00284 {
00285     int *nodeId = NULL;
00286 
00287     assert( id );
00288     assert( stream );
00289 
00290     nodeId = (int*) id;
00291     printf( "node id %d, height %lf \n", *nodeId, globalMesh->h[*nodeId] );
00292 }
00293 
00294 void _SurfaceMesh_Build( void* surfaceMesh, void* data ) 
00295 {
00296     SurfaceMesh* self = NULL;
00297     
00298     assert( surfaceMesh );
00299 
00300     self = (SurfaceMesh*) surfaceMesh;
00301     
00302     if( self->rank == MASTER_PROC ){
00303     
00304         Stg_Component_Initialise( self, (void*)(&(self->numNodes)), False );
00305         
00306         globalMesh = self;  
00307     
00308         fprintf(stderr,"IN SURFACEMESH_BUILD \n");
00309         Journal_DPrintf( self->debug, "In %s - for %s\n", __func__, self->name );
00310         
00311         Stg_Component_Build( self->surfaceMeshLoader, NULL, False );
00312         SurfaceMeshLoader_LoadMesh( self->surfaceMeshLoader, self );
00313     
00314         _SurfaceMesh_FindNeighbours( self );
00315 
00316         SurfaceMesh_BuildRiverNetwork( self );
00317     }
00318 }
00319 
00320 void SurfaceMesh_BuildRiverNetwork( SurfaceMesh *self )
00321 {
00322     int i = 0;
00323     assert( self );
00324     
00325     /* determine the water flows (each node gives water to its lowest natural neighbour) */ 
00326     for( i=0; i<self->numNodes; i++ ){
00327         LinkedList_DeleteAllNodes( self->nodeProviders[i] );
00328     }
00329     _SurfaceMesh_DetermineFlows( self, self->boundaryConditions );
00330     printf("Determined flows \n");
00331 
00332     /* find the sill nodes and determine lake nodes */
00333     _SurfaceMesh_FindDonors( self );
00334     printf("Found Donors \n");
00335 
00336     /* determine the depth and size of each node in the list */ 
00337     _SurfaceMesh_DetermineCatchmentSizes( self );   
00338 
00339     printf("Determined CatchmentSizes \n");
00340 }
00341 
00342 void _SurfaceMesh_Initialise( void* surfaceMesh, void* data )
00343 {
00344     int i;
00345     SurfaceMesh *mesh = NULL;
00346     int arrayElements = 0;
00347 
00348     mesh = ( SurfaceMesh* )surfaceMesh;
00349     assert( mesh );
00350     assert( data );
00351     
00352     arrayElements = *(int*)data;
00353 
00354     assert( arrayElements > 0 );
00355     
00356     mesh->arrayElements = arrayElements;
00357     mesh->haloNodes = (SurfaceMeshHaloNodes*)malloc( mesh->numProcs * sizeof( SurfaceMeshHaloNodes ) );
00358     memset( mesh->haloNodes, 0, sizeof( SurfaceMeshHaloNodes ) * mesh->numProcs );
00359     
00360     mesh->foreignHaloNodes = (SurfaceMeshForeignHaloNodes*)malloc( mesh->numProcs * sizeof( SurfaceMeshForeignHaloNodes ) );
00361     memset( mesh->foreignHaloNodes, 0, sizeof( SurfaceMeshForeignHaloNodes ) * mesh->numProcs );
00362 
00363     mesh->id = (int*)malloc( arrayElements*sizeof(int) ); memset( mesh->id, 0, arrayElements*sizeof(int) );
00364     mesh->sortedId = (int*)malloc( arrayElements*sizeof(int) ); memset( mesh->sortedId, 0, arrayElements*sizeof(int) );
00365     mesh->nodeType = (int*)malloc( arrayElements*sizeof(int) ); memset( mesh->nodeType, 0, arrayElements*sizeof(int) );
00366     mesh->x = (float*)malloc( arrayElements*sizeof(float) ); memset( mesh->x, 0, arrayElements*sizeof(float) );
00367     mesh->y = (float*)malloc( arrayElements*sizeof(float) ); memset( mesh->y, 0, arrayElements*sizeof(float) );
00368     mesh->h = (float*)malloc( arrayElements*sizeof(float) ); memset( mesh->h, 0, arrayElements*sizeof(float) );
00369     mesh->hi = (float*)malloc( sizeof(float)*arrayElements ); memset( mesh->hi, 0, arrayElements*sizeof(float) );
00370     mesh->h0 = (float*)malloc( sizeof(float)*arrayElements ); memset( mesh->h0, 0, arrayElements*sizeof(float) );
00371     
00372     mesh->receiver = (int*)malloc( sizeof(int)*arrayElements ); memset( mesh->receiver, 0, sizeof(int)*arrayElements);
00373     mesh->highestNeighbour = (int*)malloc( sizeof(int)*arrayElements ); memset( mesh->highestNeighbour, 0, sizeof(int)*arrayElements );
00374     mesh->processor = (int*)malloc( sizeof(int)*arrayElements ); memset( mesh->processor, 0, arrayElements*sizeof(int) );
00375     mesh->size = (int*)malloc( sizeof(int)*arrayElements ); memset( mesh->size, 0, arrayElements*sizeof(int) );
00376     mesh->nLake = (int*)malloc( sizeof(int)*arrayElements ); memset( mesh->nLake, 0, arrayElements*sizeof(int) );
00377     mesh->nCatch = (int*)malloc( sizeof(int)*arrayElements ); memset( mesh->nCatch, 0, arrayElements*sizeof(int) );
00378     mesh->slope = (float*)malloc( sizeof(float)*arrayElements ); memset( mesh->slope, 0, arrayElements*sizeof(float) );
00379     mesh->length = (float*)malloc( sizeof(float)*arrayElements ); memset( mesh->length, 0, arrayElements*sizeof(float) );
00380     mesh->surface = (float*)malloc( sizeof(float)*arrayElements ); memset( mesh->surface, 0, arrayElements*sizeof(float) );
00381 
00382     mesh->nodeNeighbours = (int**)malloc( arrayElements*sizeof(int*) );
00383     mesh->nodeNeighbours[0] = (int*)malloc( mesh->maxNeighboursPerNode*arrayElements*sizeof(int));
00384     memset( mesh->nodeNeighbours[0], 0, mesh->maxNeighboursPerNode*arrayElements*sizeof(int) );
00385     for( i=0; i<arrayElements; i++ ) {
00386         mesh->nodeNeighbours[i]= mesh->nodeNeighbours[0]+ i * mesh->maxNeighboursPerNode;
00387     }
00388     
00389     mesh->vertices = (int**)malloc( 3*(arrayElements)*sizeof(int*) );
00390     mesh->vertices[0] = (int*)malloc( 3*3*arrayElements*sizeof(int));
00391     memset( mesh->vertices[0], 0, sizeof( int )*arrayElements* 3 * 3 );
00392     for( i=0; i<arrayElements*3; i++ ) {
00393         mesh->vertices[i]= mesh->vertices[0]+i*3;
00394     }
00395 
00396     mesh->sides = (float**)malloc( arrayElements*sizeof(float*) );
00397     mesh->sides[0] = (float*)malloc( mesh->maxNeighboursPerNode * arrayElements * sizeof(float));
00398     memset( mesh->sides[0], 0, sizeof( float ) * arrayElements * mesh->maxNeighboursPerNode );
00399     for( i=0; i<arrayElements; i++ ) {
00400         mesh->sides[i]= mesh->sides[0]+i*mesh->maxNeighboursPerNode;
00401     }
00402     
00403     mesh->nodeProviders = (LinkedList**)malloc( sizeof( LinkedList* ) * arrayElements );
00404     for( i=0; i<arrayElements; i++ ){
00405         mesh->nodeProviders[i] = LinkedList_New(providersCompareFunction,
00406                                                 providersDataCopyFunction,
00407                                                 providersDataPrintFunction,
00408                                                 NULL,
00409                                                 LINKEDLIST_UNSORTED);
00410     }
00411     
00412     mesh->boundaryConditions = (float*)malloc( mesh->numNodes*sizeof(float) );
00413     memset( mesh->boundaryConditions, 0, sizeof( float ) * mesh->numNodes );
00414     
00415     mesh->mapGlobalToLocal = (int*) malloc( mesh->numNodes * sizeof(int) ); /* array that maps from global node number to local node index */
00416     mesh->numNeigh = (int*) malloc( mesh->numNodes * sizeof(int) );         /* array of number of neighbours of each node  */
00417 
00418     mesh->foreignProviderChanges = (LinkedList**) malloc( sizeof( LinkedList* ) * mesh->numProcs );
00419     mesh->localProviderChanges = (LinkedList**) malloc( sizeof( LinkedList* ) * mesh->numProcs );
00420 
00421     for( i=0; i<mesh->numProcs; i++ ){
00422         mesh->foreignProviderChanges[i] = LinkedList_New(providerSyncCompareFunction,
00423                                                 providerSyncDataCopyFunction,
00424                                                 providerSyncDataPrintFunction,
00425                                                 NULL,
00426                                                 LINKEDLIST_UNSORTED);
00427 
00428         mesh->localProviderChanges[i] = LinkedList_New(providerSyncCompareFunction,
00429                                                 providerSyncDataCopyFunction,
00430                                                 providerSyncDataPrintFunction,
00431                                                 NULL,
00432                                                 LINKEDLIST_UNSORTED);
00433     }
00434     
00435     mesh->requestTable = (MPI_Request*)malloc( sizeof( MPI_Request ) * arrayElements );
00436     mesh->statusTable = (MPI_Status*)malloc( sizeof( MPI_Status ) * arrayElements );
00437 }
00438 
00439 void _SurfaceMesh_Execute( void* surfaceMesh, void* data ) {
00440     /* Do nothing */
00441 }
00442 
00443 void _SurfaceMesh_Destroy( void* surfaceMesh, void* data ) {
00444     /* Do nothing */
00445 }
00446 
00447 /* Private functions */
00448 
00449 void _SurfaceMesh_FindNeighbours( SurfaceMesh* mesh ) {
00450     int* vis_tlist;
00451     int* vis_elist;
00452     int* add_tlist;
00453     float *bb, *xy;
00454     float **PP, **aa;
00455     int *nodes, **neighbour;
00456     float *newsurface;
00457     double precision = 0.0, **points;
00458     Index i;
00459     
00460     /* Allocate temporary working arrays */
00461     
00462     points = (double**)malloc( mesh->numNodes*sizeof(double*) );
00463     points[0] = (double*)malloc( (2)*mesh->numNodes*sizeof(double));
00464     for( i=0; i<mesh->numNodes; i++ ) {
00465         points[i]= points[0]+i*2;
00466     }
00467     
00468     neighbour = (int**)malloc( 3*mesh->numNodes*sizeof(int*) );
00469     neighbour[0] = (int*)malloc( 3*3*mesh->numNodes*sizeof(int));
00470     for( i=0; i<mesh->numNodes*3; i++ ) {
00471         neighbour[i] = neighbour[0]+i*3;
00472     }
00473         
00474     vis_tlist = (int*)malloc( mesh->numNodes*sizeof(int) );
00475     vis_elist = (int*)malloc( mesh->numNodes*sizeof(int) );
00476     add_tlist = (int*)malloc( mesh->numNodes*sizeof(int) );
00477     nodes = (int*)malloc( mesh->numNodes*sizeof(int) );
00478     newsurface = (float*)malloc( mesh->numNodes*sizeof(float));
00479     
00480     bb = (float*)malloc( mesh->maxNeighboursPerNode*sizeof(float) );
00481     xy = (float*)malloc( 2*sizeof(float) );
00482     PP = (float**)malloc( mesh->maxNeighboursPerNode*sizeof(float*) );
00483     PP[0] = (float*)malloc( 2*mesh->maxNeighboursPerNode*sizeof(float));
00484     for( i=0; i<mesh->maxNeighboursPerNode; i++ ) {
00485         PP[i]= PP[0]+i*2;
00486     }
00487     
00488     aa = (float**)malloc( 2*sizeof(float*) );
00489     aa[0] = (float*)malloc( 2*mesh->maxNeighboursPerNode*sizeof(float));
00490     for( i=0; i<2; i++ ) {
00491         aa[i]= aa[0]+i*mesh->maxNeighboursPerNode;
00492     }
00493 
00494     find_neighbours( mesh->x, mesh->y, mesh->nodeNeighbours[0], mesh->numNeigh, &mesh->numNodes,
00495         &mesh->maxNeighboursPerNode, points[0], mesh->vertices[0], neighbour[0], nodes,
00496         vis_tlist, vis_elist, add_tlist, &mesh->elementGlobalCount,
00497         mesh->surface, newsurface, &precision, xy, 
00498         PP[0], aa[0], bb, &mesh->surfScale, mesh->sides[0]  );
00499 
00500     printf ("Found neighbours \n");
00501 
00502     /* Free temporary working arrays */
00503     free( vis_tlist );
00504     free( vis_elist );
00505     free( add_tlist );
00506     free( nodes );
00507     free( newsurface );
00508     free( PP[0] );
00509     free( PP );
00510     free( aa[0] );
00511     free( points[0] );
00512     free( points );
00513     free( neighbour[0] );
00514     free( neighbour );
00515     free( aa );
00516     free( bb );
00517     free( xy );
00518 }
00519 
00520 void SurfaceMesh_AllocateMemoryForHaloNodes( SurfaceMesh *mesh )
00521 {
00522     int i;
00523 
00524     assert( mesh );
00525     
00526     for( i=0; i<mesh->numProcs; i++ ){
00527         if( mesh->haloNodes[i].numHaloNodes > 0 ){
00528             
00529             mesh->haloNodes[i].haloId = (int*)malloc( mesh->haloNodes[i].numHaloNodes*sizeof(int) );
00530             memset( mesh->haloNodes[i].haloId, 0, mesh->haloNodes[i].numHaloNodes*sizeof(int) );
00531         
00532             mesh->haloNodes[i].x = (float*)malloc( mesh->haloNodes[i].numHaloNodes*sizeof(float) );
00533             memset( mesh->haloNodes[i].x, 0, mesh->haloNodes[i].numHaloNodes*sizeof(float) );
00534         
00535             mesh->haloNodes[i].y = (float*)malloc( mesh->haloNodes[i].numHaloNodes*sizeof(float) );
00536             memset( mesh->haloNodes[i].y, 0, mesh->haloNodes[i].numHaloNodes*sizeof(float) );
00537         
00538             mesh->haloNodes[i].h = (float*)malloc( mesh->haloNodes[i].numHaloNodes*sizeof(float) );
00539             memset( mesh->haloNodes[i].h, 0, sizeof(float)*mesh->haloNodes[i].numHaloNodes );
00540         
00541             mesh->haloNodes[i].receiver = (int*)malloc( sizeof(int)*mesh->haloNodes[i].numHaloNodes );
00542             memset( mesh->haloNodes[i].receiver, 0, sizeof(int)*mesh->haloNodes[i].numHaloNodes );
00543         
00544             mesh->haloNodes[i].surface = (float*)malloc( sizeof(float)*mesh->haloNodes[i].numHaloNodes );
00545             memset( mesh->haloNodes[i].surface, 0, sizeof(float)*mesh->haloNodes[i].numHaloNodes );
00546         }
00547         
00548         if( mesh->foreignHaloNodes[i].numForeignHaloNodes > 0 ){
00549             mesh->foreignHaloNodes[i].foreignHaloId = (int*)malloc( sizeof( int ) * mesh->foreignHaloNodes[i].numForeignHaloNodes );
00550             memset(mesh->foreignHaloNodes[i].foreignHaloId, 0, sizeof( int ) * mesh->foreignHaloNodes[i].numForeignHaloNodes );
00551         }
00552     }
00553 }
00554 
00555 void SurfaceMesh_ReleaseMemory( SurfaceMesh* mesh )
00556 {
00557     int i = 0;
00558 
00559     assert( mesh );
00560     
00561     if( ((Stg_Component*)mesh)->isInitialised ){
00562 
00563         free( mesh->id );
00564         free( mesh->sortedId );
00565         free( mesh->nodeType );
00566         free( mesh->x );
00567         free( mesh->y );
00568         free( mesh->h );
00569         free( mesh->hi );
00570         free( mesh->h0 );
00571     
00572         free( mesh->receiver );
00573         free( mesh->highestNeighbour );
00574         free( mesh->processor );
00575         free( mesh->size );
00576         free( mesh->nLake );
00577         free( mesh->nCatch );
00578         free( mesh->slope );
00579         free( mesh->length );
00580         free( mesh->surface );
00581 
00582         free( mesh->nodeNeighbours[0] );
00583         free( mesh->nodeNeighbours );
00584 
00585         free( mesh->vertices[0] );
00586         free( mesh->vertices );
00587     
00588         free( mesh->sides[0] );
00589         free( mesh->sides );
00590 
00591         for( i=0; i<mesh->arrayElements; i++ ){
00592             Stg_Class_Delete( mesh->nodeProviders[i] );
00593         }
00594         free( mesh->nodeProviders );
00595     
00596         free( mesh->boundaryConditions );
00597     
00598         free( mesh->mapGlobalToLocal );
00599         free( mesh->numNeigh );
00600     
00601         for( i=0; i<mesh->numProcs; i++ ){
00602             Stg_Class_Delete( mesh->foreignProviderChanges[i] );
00603             Stg_Class_Delete( mesh->localProviderChanges[i] );
00604         }
00605         free( mesh->foreignProviderChanges );
00606         free( mesh->localProviderChanges );
00607 
00608         free( mesh->requestTable );
00609         free( mesh->statusTable );
00610     }
00611     
00612     if( mesh->haloNodes ){
00613         for( i=0; i<mesh->numProcs; i++ ){
00614             if( mesh->haloNodes[i].haloId ){
00615                 free( mesh->haloNodes[i].haloId );
00616             }
00617             if( mesh->haloNodes[i].x ){
00618                 free( mesh->haloNodes[i].x );
00619             }
00620             if( mesh->haloNodes[i].y ){
00621                 free( mesh->haloNodes[i].y );
00622             }
00623             if( mesh->haloNodes[i].h ){
00624                 free( mesh->haloNodes[i].h );
00625             }
00626             if( mesh->haloNodes[i].receiver ){
00627                 free( mesh->haloNodes[i].receiver );
00628             }
00629             if( mesh->haloNodes[i].surface ){
00630                 free( mesh->haloNodes[i].surface );
00631             }
00632         }
00633         free( mesh->haloNodes );
00634     }
00635     
00636     if( mesh->foreignHaloNodes ){
00637         for( i=0; i<mesh->numProcs; i++ ){
00638             if( mesh->foreignHaloNodes[i].foreignHaloId ){
00639                 free( mesh->foreignHaloNodes[i].foreignHaloId );
00640             }
00641         }
00642         free( mesh->foreignHaloNodes );
00643     }
00644 }
00645 
00646 void _SurfaceMesh_DetermineFlows( SurfaceMesh* mesh, float* boundaryConditions )
00647 {
00648     int i;
00649     float dh, dhmin, dd, dhmax;
00650     int *currNode, neighbourNode;
00651     BTreeIterator *iterator = NULL;
00652     Stream *stream;
00653     BTree *startPtr = NULL;
00654 
00655     stream = Journal_Register( InfoStream_Type, "myStream" );
00656     /*Print( sPtr, stream );*/
00657     
00658     if( mesh->rank == 0 ){
00659         startPtr = BTree_New(
00660                         heightCompareFunc,
00661                         NULL,
00662                         NULL,
00663                         nodePrintFunction,
00664                         BTREE_ALLOW_DUPLICATES);
00665     
00666         for (i=0; i<mesh->numNodes; i++) {
00667             mesh->size[i] = 0;
00668             BTree_InsertNode( startPtr, &mesh->id[i], sizeof( int ) );
00669         }   
00670 
00671         iterator = BTreeIterator_New( startPtr );
00672     
00673         for (currNode = (int*)BTreeIterator_First( iterator ); 
00674                 currNode!=NULL; 
00675                 currNode = (int*)BTreeIterator_Next( iterator )) {
00676         
00677             mesh->receiver[*currNode-1] = *currNode;
00678             dhmin = 0.0;
00679             dhmax = -1.0e32;
00680             
00681             for (i=0; i<mesh->numNeigh[*currNode-1]; i++) {
00682                 dh = 0.0;
00683                 neighbourNode = mesh->nodeNeighbours[*currNode -1][i];
00684                 if ( (mesh->h[neighbourNode-1] < mesh->h[*currNode-1]) ) {
00685                     dh = mesh->h[neighbourNode-1] - mesh->h[*currNode-1];
00686                     dd = sqrt( pow((mesh->x[neighbourNode-1] - mesh->x[*currNode-1]),2) + pow((mesh->y[neighbourNode-1] - mesh->y[*currNode-1]),2) );
00687                     if (dh/dd < dhmin) {
00688                         dhmin = dh/dd;  
00689                         mesh->receiver[*currNode-1] = neighbourNode;
00690                     }
00691                     
00692                     if (dh/dd > dhmax) {
00693                         dhmax = dh/dd;  
00694                         mesh->highestNeighbour[*currNode-1] = neighbourNode;
00695                     }
00696                 }
00697             }   
00698 
00699             neighbourNode = mesh->receiver[*currNode-1];
00700             if (*currNode != neighbourNode) {   
00701                 /* insert currNode to list of neighbourNode's providers */
00702                 LinkedList_InsertNode(mesh->nodeProviders[neighbourNode-1], (void*)(currNode), sizeof(int));
00703             }
00704             else {
00705             }
00706         }
00707         Stg_Class_Delete( iterator );
00708         Stg_Class_Delete( startPtr );
00709     }
00710     else{
00711         
00712     }
00713 }
00714 
00715 void _SurfaceMesh_FindDonors( SurfaceMesh* mesh ) {
00716     int i; 
00717     float hmax, dx, dy;
00718     
00719     if( mesh->rank == 0 ){
00720         
00721         /* Find slope and length , as done in find_donors.f in Cascade2002 */
00722         for (i=0; i<mesh->numNodes; i++) {
00723             if (mesh->receiver[i] == mesh->id[i]) {
00724                 mesh->slope[i] = 0;
00725                 mesh->length[i] = 0;
00726             }
00727             else {
00728                 dx = mesh->x[i] - mesh->x[mesh->receiver[i]-1];
00729                 dy = mesh->y[i] - mesh->y[mesh->receiver[i]-1];
00730                 mesh->length[i] = sqrt(dx*dx+dy*dy);
00731                 mesh->slope[i] = (mesh->h[mesh->receiver[i]-1] - mesh->h[i]) / mesh->length[i];
00732             }
00733         }
00734     
00735         hmax = mesh->h[0];
00736         for (i=0; i<mesh->numNodes; i++) {
00737             mesh->nLake[i] = 0;   /* initialize nlake to 0 for all nodes (no lake nodes) */
00738             if (mesh->h[i] > hmax){
00739                 hmax = mesh->h[i];
00740             }
00741         }
00742     }
00743     else{
00744         
00745     }
00746     return;
00747 }
00748 
00749 int getSize( SurfaceMesh *mesh, int nodeId )
00750 {
00751     int i = 0;
00752     int *result = NULL;
00753     int size = 1;
00754     int q;
00755     
00756     assert( mesh );
00757 
00758     /* recursively going uptree to calculate the number of nodes above node n */
00759     for (i=0; i<mesh->nodeProviders[nodeId-1]->nodeCount; i++) {
00760         result = LinkedList_ReturnNodeDataAt( mesh->nodeProviders[nodeId-1], i, int );
00761         assert( *result );
00762 
00763         q = *result;
00764         if( mesh->h[nodeId-1] < mesh->h[q-1] ){
00765             size += getSize( mesh, q );
00766         }
00767     }
00768 
00769     return size;
00770 }
00771 
00772 void _SurfaceMesh_DetermineCatchmentSizes( SurfaceMesh* mesh ) {
00773     int *p;
00774     int i = 0;
00775     BTreeIterator *iter_i = NULL;
00776     BTree *startPtr = NULL;
00777 
00778     
00779 
00780     if( mesh->rank == 0 ){
00781             startPtr = BTree_New(
00782                     heightCompareFunc,
00783                     NULL,
00784                     NULL,
00785                     nodePrintFunction,
00786                     BTREE_ALLOW_DUPLICATES);
00787     
00788         for (i=0; i<mesh->numNodes; i++) {
00789             mesh->size[i] = 0;
00790             BTree_InsertNode( startPtr, &mesh->id[i], sizeof( int ) );
00791         }   
00792 
00793         iter_i = BTreeIterator_New( startPtr );
00794         for( p = (int*)BTreeIterator_First( iter_i ); p != NULL; p = (int*)BTreeIterator_Next( iter_i ) ){
00795             mesh->size[*p-1] = getSize( mesh, *p );
00796         }
00797         
00798         Stg_Class_Delete( (void*)iter_i );
00799         Stg_Class_Delete( (void*) startPtr );
00800     }
00801     else{
00802         
00803     }
00804 }
00805 
00806 void SurfaceMesh_UpdateFlows( SurfaceMesh *mesh )
00807 {
00808     int i, j;
00809     int newReceiver, oldReceiver;
00810     int localIndex = 0;
00811     int neighbour = 0;
00812     int haloNodeProc, haloNodeIndex;
00813     providerSyncStruct data, *ps;
00814     float dx, dy, dh, dd, dhmin, gradient, dhmax;
00815     LinkedListIterator *iterator = NULL;
00816     void *result = NULL;
00817 
00818     assert( mesh );
00819     
00820     for( i=0; i<mesh->numProcs; i++ ){
00821         if( mesh->rank != i ){
00822             LinkedList_DeleteAllNodes( mesh->foreignProviderChanges[i] );
00823             LinkedList_DeleteAllNodes( mesh->localProviderChanges[i] );
00824         }
00825     }
00826     
00827     for( i=0; i<mesh->myLoad; i++ ){
00828         
00829         if( mesh->boundaryConditions[mesh->id[i]-1] != 0 ){
00830             
00831             dhmin = 0;
00832             dhmax = -1.0e32;
00833             newReceiver = mesh->id[i];
00834 
00835             for( j=0; j<mesh->numNeigh[mesh->id[i] - 1]; j++ ){
00836                 neighbour = mesh->nodeNeighbours[i][j];
00837             
00838                 if( neighbour != mesh->id[i] ){
00839 
00840                     localIndex = mesh->mapGlobalToLocal[neighbour-1];
00841 
00842                     if( localIndex == mesh->myLoad ){
00843                         SurfaceMesh_FindHaloNode( mesh, neighbour, &haloNodeProc, &haloNodeIndex );
00844 
00845                         dx = mesh->haloNodes[haloNodeProc].x[haloNodeIndex] - mesh->x[i];
00846                         dy = mesh->haloNodes[haloNodeProc].y[haloNodeIndex] - mesh->y[i];
00847                         dh = mesh->haloNodes[haloNodeProc].h[haloNodeIndex] - mesh->h[i];
00848                     }
00849                     else{
00850                         dx = mesh->x[localIndex] - mesh->x[i];
00851                         dy = mesh->y[localIndex] - mesh->y[i];
00852                         dh = mesh->h[localIndex] - mesh->h[i];
00853                     }
00854             
00855                     dd = sqrt( dx*dx + dy*dy );
00856                     gradient = dh/dd;
00857             
00858                     if( gradient < dhmin ){
00859                         dhmin = gradient;
00860                         newReceiver = neighbour;
00861                     }
00862                     
00863                     if( gradient > dhmax ){
00864                         dhmax = gradient;
00865                         mesh->highestNeighbour[i] = neighbour;
00866                     }
00867                 }
00868             }
00869 
00870             if( newReceiver != mesh->receiver[i] ){
00871                 
00872                 oldReceiver = mesh->receiver[i];
00873                 mesh->receiver[i] = newReceiver;
00874             
00875                 localIndex = mesh->mapGlobalToLocal[oldReceiver-1];
00876 
00877                 if( localIndex == mesh->myLoad ){
00878                     SurfaceMesh_FindHaloNode( mesh, oldReceiver, &haloNodeProc, &haloNodeIndex );
00879                     
00880                     data.nodeId = oldReceiver;
00881                     data.provider = -mesh->id[i];
00882                     LinkedList_InsertNode( mesh->foreignProviderChanges[haloNodeProc], (void*)&data, sizeof(providerSyncStruct) );
00883                 }
00884                 else{
00885                     LinkedList_DeleteNode( mesh->nodeProviders[localIndex], (void*)&(mesh->id[i]) );
00886                 }
00887             
00888                 localIndex = mesh->mapGlobalToLocal[newReceiver-1];
00889             
00890                 if( localIndex == mesh->myLoad ){
00891                     SurfaceMesh_FindHaloNode( mesh, newReceiver, &haloNodeProc, &haloNodeIndex );
00892                     
00893                     data.nodeId = newReceiver;
00894                     data.provider = mesh->id[i];
00895                     LinkedList_InsertNode( mesh->foreignProviderChanges[haloNodeProc], (void*)&data, sizeof(providerSyncStruct) );
00896                 }
00897                 else{
00898                     LinkedList_InsertNode( mesh->nodeProviders[localIndex], (void*)&(mesh->id[i]), sizeof(int) );
00899                 }
00900             }
00901         }
00902     }
00903 
00904 #ifdef USE_DEBUG
00905     printf( "Foreign Provider Changes\n" );
00906     for( i=0; i<mesh->numProcs; i++ ){
00907         if( mesh->rank != i ){
00908             Print( mesh->foreignProviderChanges[i], Journal_Register( Info_Type, "stream" ) );
00909         }
00910     }
00911 #endif
00912     
00913     for( i=0; i<mesh->numProcs; i++ ){
00914         if( mesh->rank != i ){
00915             NonBlockingSend_LinkedList( mesh->foreignProviderChanges[i], sizeof( providerSyncStruct ), i,
00916                   PROVIDERS_SYNC_TAG, MPI_COMM_WORLD, mesh->requestTable, mesh->statusTable );
00917         }
00918     }
00919     
00920     for( i=0; i<mesh->numProcs; i++ ){
00921         if( mesh->rank != i ){
00922             NonBlockingReceive_LinkedList( mesh->localProviderChanges[i], sizeof( providerSyncStruct ), i,
00923                   PROVIDERS_SYNC_TAG, MPI_COMM_WORLD, mesh->requestTable, mesh->statusTable );
00924         }
00925     }
00926 
00927 #ifdef USE_DEBUG
00928     printf( "Local Provider Changes\n" );
00929     for( i=0; i<mesh->numProcs; i++ ){
00930         if( mesh->rank != i ){
00931             Print( mesh->localProviderChanges[i], Journal_Register( Info_Type, "stream" ) );
00932         }
00933     }
00934 #endif
00935 
00936     for( i=0; i<mesh->numProcs; i++ ){
00937         
00938         if( i != mesh->rank ){
00939             iterator = LinkedListIterator_New( mesh->localProviderChanges[i] );
00940 
00941             for( result = LinkedListIterator_First( iterator ); result; result = LinkedListIterator_Next( iterator ) ){
00942                 ps = (providerSyncStruct*)result;
00943                 localIndex = mesh->mapGlobalToLocal[ps->nodeId-1];
00944                 
00945                 assert( localIndex != mesh->myLoad );
00946                 assert( mesh->id[localIndex] == ps->nodeId );
00947 
00948                 if( ps->provider < 0 ){
00949                     ps->provider = -ps->provider;
00950                     LinkedList_DeleteNode( mesh->nodeProviders[localIndex], (void*)&(ps->provider) );
00951                 }
00952                 else{
00953                     LinkedList_InsertNode( mesh->nodeProviders[localIndex], (void*)&(ps->provider), sizeof( int ) );
00954                     if( mesh->nodeProviders[localIndex]->nodeCount > mesh->numNeigh[ps->nodeId-1] ){
00955                         fprintf( stderr, "Number of providers cannot be greater than number of natural neighbours..!\n" );
00956                         assert( 0 );
00957                     }
00958                 }
00959             }
00960                 
00961             Stg_Class_Delete( iterator );
00962         }
00963     }
00964     
00965     SurfaceMesh_OrderNodes( mesh );
00966     
00967     for( i=0; i<mesh->myLoad; i++ ){
00968         if( mesh->receiver[i] == mesh->id[i] ){
00969             mesh->slope[i] = 0;
00970             mesh->length[i] = 0;
00971         }
00972         else{
00973             localIndex = mesh->mapGlobalToLocal[mesh->receiver[i] - 1];
00974             
00975             if( localIndex == mesh->myLoad ){
00976                 SurfaceMesh_FindHaloNode( mesh, mesh->receiver[i], &haloNodeProc, &haloNodeIndex );
00977                 
00978                 dx = mesh->x[i] - mesh->haloNodes[haloNodeProc].x[haloNodeIndex];
00979                 dy = mesh->y[i] - mesh->haloNodes[haloNodeProc].y[haloNodeIndex];
00980                 mesh->length[i] = sqrt( dx*dx + dy*dy );
00981                 mesh->slope[i] = (mesh->haloNodes[haloNodeProc].h[haloNodeIndex] - mesh->h[i]) / mesh->length[i];
00982             }
00983             else{
00984                 dx = mesh->x[i] - mesh->x[localIndex];
00985                 dy = mesh->y[i] - mesh->y[localIndex];
00986                 mesh->length[i] = sqrt( dx*dx + dy*dy );
00987                 mesh->slope[i] = (mesh->h[localIndex] - mesh->h[i]) / mesh->length[i];
00988             }
00989         }
00990     }
00991 
00992 #ifdef USE_DEBUG
00993     for( i=0; i<mesh->myLoad; i++ ){
00994         if( mesh->nodeProviders[i]->nodeCount > mesh->numNeigh[mesh->id[i]-1] ){
00995             printf( "node[%d]->providersCount = %d\n", mesh->id[i], mesh->nodeProviders[i]->nodeCount );
00996             Stg_Class_Print( mesh->nodeProviders[i], Journal_Register( Info_Type, "stream" ) );
00997             printf( "Number of providers cannot be greater than number of natural neighbours..!\n" );
00998             assert( 0 );
00999         }
01000     } 
01001 #endif 
01002 }
01003 
01004 void SurfaceMesh_OrderNodes( SurfaceMesh *mesh )
01005 {
01006     assert( mesh );
01007     
01008     localMesh = mesh;
01009     qsort( mesh->sortedId, mesh->myLoad, sizeof(int), idCompareFunc );
01010 }
01011 
01012 void SurfaceMesh_FindHaloNode( SurfaceMesh *mesh, int nodeId, int *processor, int *index )
01013 {
01014     int i, j;
01015     int found = 0;
01016     
01017     assert( mesh );
01018 
01019     for( i=0; i<mesh->numProcs; i++ ){
01020         for( j=0; j<mesh->haloNodes[i].numHaloNodes; j++ ){
01021             if( mesh->haloNodes[i].haloId[j] == nodeId ){
01022                 *processor = i;
01023                 *index = j;
01024                 found = 1;
01025                 break;
01026             }
01027         }
01028     }
01029     assert( found );
01030 }
01031 
01032 void SurfaceMesh_BoundaryConditions( SurfaceMesh *surfaceMesh )
01033 {
01034 
01035 }
01036