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

_Interpolator.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 */
00032 #include <mpi.h>
00033 #include <StGermain/StGermain.h>
00034 #include <Cascade/cascade.h>
00035 
00036 #include "types.h"
00037 
00038 #include <stdio.h>
00039 #include <stdlib.h>
00040 #include <string.h>
00041 #include <assert.h>
00042 #include <math.h>
00043 #include "Catchment.h"
00044 #include "SurfaceMesh.h"
00045 #include "SurfaceRegularMesh.h"
00046 #include "Context.h"
00047 #include "_Interpolator.h"
00048 
00049 /* Textual name of this class */
00050 const Type _Interpolator_Type = "_Interpolator";
00051 
00052 _Interpolator                   *_Interpolator_New( SizeT                           _sizeOfSelf,
00053                                     Type                            type,
00054                                     Stg_Class_DeleteFunction*               _delete,
00055                                     Stg_Class_PrintFunction*                _print,
00056                                     Stg_Class_CopyFunction*             _copy, 
00057                                     Stg_Component_DefaultConstructorFunction*   _defaultConstructor,
00058                                     Stg_Component_ConstructFunction*            _construct,
00059                                     Stg_Component_BuildFunction*        _build,
00060                                     Stg_Component_InitialiseFunction*       _initialise,
00061                                     Stg_Component_ExecuteFunction*      _execute,
00062                                     Stg_Component_DestroyFunction*      _destroy,
00063                                     Name                            name,
00064                                     Bool                            initFlag,
00065                                     Dictionary                      *dictionary,
00066                                     SurfaceMesh                     *mesh,
00067                                     SurfaceRegularMesh              *regularMesh,
00068                                     Interpolator_InterpolateFromGridToMeshFunction *interpolateG2MFunc)
00069 {
00070     _Interpolator* self;
00071     
00072     /* Allocate memory */
00073     assert( _sizeOfSelf >= sizeof(_Interpolator) );
00074     self = (_Interpolator*)_Stg_Component_New( _sizeOfSelf, type, _delete, _print, _copy, _defaultConstructor, _construct, _build, 
00075             _initialise, _execute, _destroy, name, NON_GLOBAL );
00076 
00077     assert( dictionary );
00078     assert( mesh );
00079     assert( regularMesh );
00080     assert( interpolateG2MFunc );
00081         
00082     self->dictionary = dictionary;
00083     self->mesh = mesh; 
00084     self->regularMesh = regularMesh;
00085     self->interpolateG2MFunc = interpolateG2MFunc;
00086 
00087     Journal_Firewall( (self->regularMesh->nx > 0), Journal_Register( Error_Type, "errorStream" ), "0 nodes in x direction..!\n" );
00088     Journal_Firewall( (self->regularMesh->ny > 0), Journal_Register( Error_Type, "errorStream" ), "0 nodes in y direction..!\n" );
00089     
00090     if( initFlag ){
00091         _Interpolator_Init( self );
00092     }
00093 
00094     return self;    
00095 }
00096 
00097 void _Interpolator_Init( _Interpolator *self )
00098 {
00099     
00100 }
00101 
00102 void _Interpolator_Print( void *_interpolator, Stream *stream )
00103 {
00104     int numNodes = 0;
00105     int i, j, count;
00106     _Interpolator *self = (_Interpolator*) _interpolator;
00107 
00108     assert( self );
00109     assert( stream );
00110     
00111     if( self->mesh->rank == MASTER_PROC ){
00112         numNodes = self->mesh->numNodes;
00113     }
00114     else{
00115         numNodes = self->mesh->myLoad;
00116     }
00117 
00118     /* print parent */
00119     _Stg_Component_Print( (void*) self, stream );
00120 
00121     /* general info */
00122     Journal_Printf( stream, "_Interpolator (ptr): (%p)\n", self );
00123 
00124     /* Virtual Info */
00125 
00126     /* _Interpolator Info */
00127     Journal_Printf( stream, "Rank - %d\n", self->mesh->rank );
00128     
00129     Journal_Printf( stream, "Node to grid point mapping lists:\n-----------------------------\n" );
00130     count = 0;
00131     for( i=0; i<self->regularMesh->nx; i++ ){
00132         for( j=0; j<self->regularMesh->ny; j++ ){
00133             Journal_Printf( stream, "Irregular nodes associated:\n" );
00134             Print( self->nodesForGridPoint[i*self->regularMesh->nx+j], stream );
00135             if( self->nodesForGridPoint[i*self->regularMesh->nx+j]->nodeCount > count ){
00136                 count = self->nodesForGridPoint[i*self->regularMesh->nx+j]->nodeCount;
00137             }
00138         }
00139     }
00140     
00141     Journal_Printf( stream, "Grid point to node mapping tables:\n-----------------------------\n" );
00142     for( i=0; i<self->mesh->numNodes; i++ ){
00143         Journal_Printf( stream, "Node id %d\n", i + 1 );
00144         Journal_Printf( stream, "Grid points:\n" );
00145         for( j=0; j<4; j++ ){
00146             Journal_Printf( stream, "\t%d, %d\n", self->gridPointsForNode[i][j].i, self->gridPointsForNode[i][j].j );
00147         }
00148     }
00149 
00150     Journal_Printf( stream, "Max number of nodes for a grid cell %d\n", count );
00151 }
00152 
00153 void _Interpolator_Delete( void *interpolator )
00154 {
00155     int i = 0, j = 0;
00156     
00157     _Interpolator *self = (_Interpolator*)interpolator;
00158 
00159     assert( self );
00160 
00161     for( i=0; i<self->regularMesh->nx; i++ ){
00162         for( j=0; j<self->regularMesh->ny; j++ ){
00163             Stg_Class_Delete( self->nodesForGridPoint[i*self->regularMesh->nx+j] );
00164         }
00165     }
00166     
00167     free( self->gridPointsForNode[0] );
00168     free( self->gridPointsForNode );
00169     
00170     if( self->mappingTable ){
00171         free( self->mappingTable[0][0] );
00172         free( self->mappingTable[0] );
00173         free( self->mappingTable );
00174     }
00175     
00176     _Stg_Component_Delete( self );
00177 }
00178 
00179 void _Interpolator_Construct( void *interpolator, Stg_ComponentFactory *cf )
00180 {
00181     
00182 }
00183 
00184 void _Interpolator_Build( void *interpolator, void *data )
00185 {
00186     float minX, minY, maxX, maxY;
00187     float minI, minJ, maxI, maxJ;
00188     float dx, dy;
00189     int i, j, k, id;
00190     GridPoint gridPoints[4];
00191     SurfaceMesh *mesh = NULL;
00192     SurfaceMesh *globalMesh = NULL;
00193     SurfaceRegularMesh *regularMesh = NULL;
00194 
00195     _Interpolator *self = (_Interpolator*)interpolator;
00196     
00197     assert( self );
00198 
00199     mesh = self->mesh;
00200     regularMesh = self->regularMesh;
00201     globalMesh = ( SurfaceMesh* )data;
00202     assert( globalMesh );
00203 
00204     Journal_Firewall( (self->regularMesh->nx > 0), Journal_Register( Error_Type, "errorStream" ), "0 nodes in x direction..!\n" );
00205     Journal_Firewall( (self->regularMesh->ny > 0), Journal_Register( Error_Type, "errorStream" ), "0 nodes in y direction..!\n" );
00206     
00207     if( mesh->rank == MASTER_PROC ){
00208         self->mappingTable = (int***)malloc( sizeof( int** ) * regularMesh->nx );
00209         self->mappingTable[0] = (int**)malloc( sizeof( int* ) * regularMesh->nx * regularMesh->ny );
00210         
00211     
00212         for( i=0; i<regularMesh->nx; i++ ){
00213             self->mappingTable[i] = self->mappingTable[0] + i * regularMesh->ny;
00214         }
00215 
00216         self->mappingTable[0][0] = (int*)malloc( sizeof( int ) * 
00217                 regularMesh->nx *
00218                 regularMesh->ny * 
00219                 mesh->numProcs );
00220         
00221         memset( self->mappingTable[0][0], 0, sizeof( int ) *
00222                 regularMesh->nx *
00223                 regularMesh->ny *
00224                 mesh->numProcs );
00225 
00226         for( i=0; i<regularMesh->nx; i++ ){
00227             for( j=0; j<regularMesh->ny; j++ ){
00228                 self->mappingTable[i][j] = self->mappingTable[0][0] +
00229                     (i*regularMesh->nx + j)*mesh->numProcs;
00230             }
00231         }
00232     }
00233     
00234     self->nodesForGridPoint = (LinkedList**)malloc( sizeof(LinkedList) * self->regularMesh->nx * self->regularMesh->ny );
00235     for( i=0; i<self->regularMesh->nx*self->regularMesh->ny; i++ ){
00236         self->nodesForGridPoint[i] = LinkedList_New( interpolatorCompareFunction,
00237                                                         interpolatorDataCopyFunction,
00238                                                         interpolatorPrintFunction,
00239                                                         interpolatorDataDeleteFunction,
00240                                                         LINKEDLIST_UNSORTED);
00241     }
00242 
00243     self->gridPointsForNode = (GridPoint**)malloc( sizeof(GridPoint*) * self->mesh->numNodes );
00244     self->gridPointsForNode[0] = (GridPoint*)malloc( sizeof(GridPoint) * 4 * self->mesh->numNodes );
00245     for( i=0; i<self->mesh->numNodes; i++ ){
00246         self->gridPointsForNode[i] = (self->gridPointsForNode[0] + i * 4 );
00247         for( j=0; j<4; j++ ){
00248             self->gridPointsForNode[i][j].i = -1;
00249             self->gridPointsForNode[i][j].j = -1;
00250         }
00251     }
00252     
00253     dx = self->mesh->sideX / (self->regularMesh->nx-1);
00254     dy = self->mesh->sideY / (self->regularMesh->ny-1);
00255     
00256     for( i=0; i<self->regularMesh->nx; i++ ){
00257         for( j=0; j<self->regularMesh->ny; j++ ){
00258         
00259             minX = self->regularMesh->regularNodes[i][j].x;
00260             minI = i;
00261             
00262             minY = self->regularMesh->regularNodes[i][j].y;
00263             minJ = j;
00264             
00265             if( self->regularMesh->regularNodes[i][j].x + dx  <= self->mesh->sideX ){
00266                 maxX = self->regularMesh->regularNodes[i][j].x + dx;
00267                 maxI = i + 1;
00268             }
00269             else{
00270                 maxX = self->regularMesh->regularNodes[i][j].x;
00271                 minX = self->regularMesh->regularNodes[i][j].x - dx;
00272                 maxI = i;
00273                 minI = i-1;
00274             }
00275             
00276             if( self->regularMesh->regularNodes[i][j].y + dy  <= self->mesh->sideY ){
00277                 maxY = self->regularMesh->regularNodes[i][j].y + dy;
00278                 maxJ = j + 1;
00279             }
00280             else{
00281                 maxY = self->regularMesh->regularNodes[i][j].y;
00282                 minY = self->regularMesh->regularNodes[i][j].y - dy;
00283                 maxJ = j;
00284                 minJ = j-1;
00285             }
00286 
00287             if( mesh->rank == MASTER_PROC ){
00288                 for( k=0; k<globalMesh->numNodes; k++ ){
00289                     if( globalMesh->x[k] >=minX && globalMesh->x[k] <= maxX ){
00290                         if( globalMesh->y[k] >=minY && globalMesh->y[k] <= maxY ){
00291                             self->mappingTable[i][j][globalMesh->processor[globalMesh->id[k]-1]] = 1;
00292                         }
00293                     }
00294                 }
00295             }
00296             
00297             for( k=0; k<self->mesh->myLoad; k++ ){
00298                 if( self->mesh->x[k] >=minX && self->mesh->x[k] <= maxX ){
00299                     if( self->mesh->y[k] >=minY && self->mesh->y[k] <= maxY ){
00300                         id = self->mesh->id[k];
00301                         LinkedList_InsertNode( self->nodesForGridPoint[i*self->regularMesh->nx + j], (void*)(&id), sizeof( int ) );
00302                         
00303                         gridPoints[0].i = minI+1; gridPoints[0].j = minJ+1;
00304                         gridPoints[1].i = maxI+1; gridPoints[1].j = minJ+1;
00305                         gridPoints[2].i = maxI+1; gridPoints[2].j = maxJ+1;
00306                         gridPoints[3].i = minI+1; gridPoints[3].j = maxJ+1;
00307 
00308                         /*printf( "node id %d, x %lf, y %lf\n grid pints:\n", id, self->mesh->x[k], self->mesh->y[k] );
00309                         for( p=0; p<4; p++ ){
00310                             printf( "\t grid point: i %d j %d\n", gridPoints[p].i, gridPoints[p].j );
00311                         }*/
00312                         
00313                         if( self->gridPointsForNode[id - 1][0].i == -1 ){
00314                             memcpy( self->gridPointsForNode[id - 1], gridPoints, sizeof(GridPoint)*4 );
00315                         }
00316                     }
00317                 }
00318             }
00319         }
00320     }
00321 }
00322 
00323 void _Interpolator_Initialise( void *interpolator, void *data )
00324 {
00325 
00326 }
00327 
00328 void _Interpolator_Execute( void *interpolator, void *data )
00329 {
00330     
00331 }
00332 
00333 void _Interpolator_Destroy( void *interpolator, void *data )
00334 {
00335     
00336 }
00337 
00338 int interpolatorCompareFunction(void *data1, void *data2)
00339 {
00340     int *d1 = NULL, *d2 = NULL;
00341 
00342     d1 = (int*)data1;
00343     d2 = (int*)data2;
00344 
00345     if (d1 == NULL || d2 == NULL){
00346         return 0;   
00347     }
00348     
00349     if (*d1 > *d2){
00350         return  1;
00351     }
00352     else if (*d1 == *d2){
00353         return 0;
00354     }
00355     else{
00356         return -1;
00357     }
00358 }
00359 
00360 void interpolatorPrintFunction( void *data, void *stream )
00361 {
00362     Stream *myStream = (Stream*)stream;
00363 
00364     Journal_Printf( myStream, "%d\n", *((int*)data) );
00365 }
00366 
00367 void interpolatorDataCopyFunction( void **nodeData, void *newData, SizeT dataSize)
00368 {
00369     *nodeData = malloc(dataSize);
00370     memset(*nodeData, 0, dataSize);
00371 
00372     memcpy(*nodeData, newData, dataSize);
00373 }
00374 
00375 void interpolatorDataDeleteFunction( void *data )
00376 {
00377     
00378 }
00379 
00380 /* public functions */
00381 
00382 void _Interpolator_InterpolateFromMeshToGrid( _Interpolator *interpolator,
00383                                             Interpolator_InterpolateFunction *interpolateFunction,
00384                                             void *arguments )
00385 {
00386     int i = 0, j = 0;
00387     double sum = 0.0;
00388     void *result = NULL;
00389     SurfaceMesh *mesh = NULL;
00390     SurfaceRegularMesh *regularMesh = NULL;
00391     int localIndex = 0;
00392     LinkedListIterator *p;
00393     float **dataArray = NULL;
00394     
00395     assert( interpolator );
00396     assert( interpolateFunction );
00397     assert( arguments );
00398     
00399     regularMesh = interpolator->regularMesh;
00400     assert( regularMesh );
00401     dataArray = regularMesh->dataArray;
00402     mesh = interpolator->mesh;
00403     assert( mesh );
00404     
00405     p = LinkedListIterator_New( interpolator->nodesForGridPoint[0] );
00406     
00407     for( i=0; i<interpolator->regularMesh->nx; i++ ){
00408         for( j=0; j<interpolator->regularMesh->ny; j++ ){
00409 
00410             p->list = interpolator->nodesForGridPoint[i*interpolator->regularMesh->nx + j];
00411             
00412             sum = 0.0;
00413             for( result = LinkedListIterator_First( p ); result; result = LinkedListIterator_Next( p ) ){
00414                 localIndex = mesh->mapGlobalToLocal[*((int*)result) - 1];
00415                 
00416                 if( localIndex != mesh->myLoad ){
00417                     sum += interpolator->mesh->surface[localIndex];
00418                 }
00419             }
00420             
00421             dataArray[i][j] = 0.0f;
00422             result = NULL;
00423                 
00424             for( result = LinkedListIterator_First( p ); result; result = LinkedListIterator_Next( p ) ){
00425                 localIndex = mesh->mapGlobalToLocal[*((int*)result) - 1];
00426 
00427                 if( localIndex != mesh->myLoad ){
00428                 
00429                     dataArray[i][j] += (interpolator->mesh->surface[localIndex]
00430                                         / ( sum )
00431                                         * interpolateFunction( localIndex, arguments ));
00432 
00433                     /*printf ( "localIndex %d, sum %lf, result %lf\n", localIndex, sum, dataArray[i][j] );*/
00434                 }
00435             }
00436         }
00437     }
00438     Stg_Class_Delete( p );
00439 }
00440 
00441 void _Interpolator_InterpolateFromGridToMesh( _Interpolator *interpolator,
00442         float **gridDeflectionArray,
00443         float *meshDeflectionArray )
00444 {
00445     assert( interpolator );
00446     
00447     assert( gridDeflectionArray );
00448     assert( meshDeflectionArray );
00449     
00450     assert( interpolator->interpolateG2MFunc );
00451     interpolator->interpolateG2MFunc( interpolator, gridDeflectionArray, meshDeflectionArray );
00452 }
00453