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

SplineInterpolator.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 */
00030 #include <mpi.h>
00031 #include <StGermain/StGermain.h>
00032 #include <Cascade/cascade.h>
00033 
00034 #include "types.h"
00035 
00036 #include <stdio.h>
00037 #include <stdlib.h>
00038 #include <string.h>
00039 #include <assert.h>
00040 #include <math.h>
00041 #include "Catchment.h"
00042 #include "SurfaceMesh.h"
00043 #include "SurfaceRegularMesh.h"
00044 #include "Context.h"
00045 #include "_Interpolator.h"
00046 #include "SplineInterpolator.h"
00047 
00048 /* Textual name of this class */
00049 const Type SplineInterpolator_Type = "SplineInterpolator";
00050 
00051 SplineInterpolator *SplineInterpolator_DefaultNew( Name name )
00052 {
00053     return _SplineInterpolator_New( sizeof( SplineInterpolator ),
00054                                 SplineInterpolator_Type,
00055                                 _SplineInterpolator_Delete,
00056                                 _SplineInterpolator_Print,
00057                                 NULL,
00058                                 (void*)SplineInterpolator_DefaultNew,
00059                                 _SplineInterpolator_Construct,
00060                                 _SplineInterpolator_Build,
00061                                 _SplineInterpolator_Initialise,
00062                                 _SplineInterpolator_Execute,
00063                                 _SplineInterpolator_Destroy,
00064                                 name,
00065                                 False,
00066                                 NULL,
00067                                 NULL,
00068                                 NULL );
00069     
00070 }
00071 
00072 SplineInterpolator *SplineInterpolator_New( Name name, Dictionary *dictionary, SurfaceMesh *mesh, SurfaceRegularMesh *regularMesh )
00073 {
00074     return _SplineInterpolator_New( sizeof( SplineInterpolator ),
00075                                 SplineInterpolator_Type,
00076                                 _SplineInterpolator_Delete,
00077                                 _SplineInterpolator_Print,
00078                                 NULL,
00079                                 (void*)SplineInterpolator_DefaultNew,
00080                                 _SplineInterpolator_Construct,
00081                                 _SplineInterpolator_Build,
00082                                 _SplineInterpolator_Initialise,
00083                                 _SplineInterpolator_Execute,
00084                                 _SplineInterpolator_Destroy,
00085                                 name,
00086                                 True,
00087                                 dictionary,
00088                                 mesh,
00089                                 regularMesh );
00090     
00091 }
00092 
00093 SplineInterpolator          *_SplineInterpolator_New( SizeT                         _sizeOfSelf,
00094                                     Type                            type,
00095                                     Stg_Class_DeleteFunction*               _delete,
00096                                     Stg_Class_PrintFunction*                _print,
00097                                     Stg_Class_CopyFunction*             _copy, 
00098                                     Stg_Component_DefaultConstructorFunction*   _defaultConstructor,
00099                                     Stg_Component_ConstructFunction*            _construct,
00100                                     Stg_Component_BuildFunction*        _build,
00101                                     Stg_Component_InitialiseFunction*       _initialise,
00102                                     Stg_Component_ExecuteFunction*      _execute,
00103                                     Stg_Component_DestroyFunction*      _destroy,
00104                                     Name                            name,
00105                                     Bool                            initFlag,
00106                                     Dictionary                      *dictionary,
00107                                     SurfaceMesh                     *mesh,
00108                                     SurfaceRegularMesh              *regularMesh )
00109 {
00110     SplineInterpolator* self;
00111     
00112     /* Allocate memory */
00113     assert( _sizeOfSelf >= sizeof(SplineInterpolator) );
00114     self = (SplineInterpolator*)_Interpolator_New( _sizeOfSelf, type, _delete, _print, _copy, _defaultConstructor, _construct, _build, 
00115             _initialise, _execute, _destroy, name, NON_GLOBAL, dictionary, mesh, regularMesh, SplineInterpolator_InterpolateFromGridToMesh );
00116 
00117     if( initFlag ){
00118         _SplineInterpolator_Init( self );
00119     }
00120 
00121     return self;
00122 }
00123 
00124 void _SplineInterpolator_Init( SplineInterpolator *self )
00125 {
00126     _Interpolator_Init( (_Interpolator*)self );
00127 }
00128 
00129 void _SplineInterpolator_Print( void *splineInterpolator, Stream *stream )
00130 {
00131     int numNodes;
00132     
00133     SplineInterpolator *self = (SplineInterpolator*) splineInterpolator;
00134 
00135     assert( self );
00136     assert( stream );
00137     
00138     if( self->mesh->rank == MASTER_PROC ){
00139         numNodes = self->mesh->numNodes;
00140     }
00141     else{
00142         numNodes = self->mesh->myLoad;
00143     }
00144 
00145     /* print parent */
00146     _Stg_Component_Print( (void*) self, stream );
00147     
00148     _Interpolator_Print( (void*) self, stream );
00149 
00150     /* general info */
00151     Journal_Printf( stream, "SplineInterpolator (ptr): (%p)\n", self );
00152 
00153     /* Virtual Info */
00154 
00155     /* SplineInterpolator Info */
00156     Journal_Printf( stream, "Rank - %d\n", self->mesh->rank );
00157 }
00158 
00159 void _SplineInterpolator_Delete( void *splineInterpolator )
00160 {
00161     SplineInterpolator *self = (SplineInterpolator*)splineInterpolator;
00162 
00163     assert( self );
00164 
00165     free( self->um );
00166     free( self->un );
00167     free( self->x1a );
00168     free( self->x2a );
00169     free( self->y2a[0] );
00170     free( self->y2a );
00171     free( self->ytmp );
00172     free( self->yytmp );
00173     free( self->y2a_t );
00174     free( self->ya_t );
00175     
00176     _Interpolator_Delete( self );
00177 }
00178 
00179 void _SplineInterpolator_Construct( void *splineInterpolator, Stg_ComponentFactory *cf )
00180 {
00181 }
00182 
00183 void _SplineInterpolator_Build( void *splineInterpolator, void *data )
00184 {
00185     int i;
00186     SplineInterpolator *self = (SplineInterpolator*)splineInterpolator;
00187 
00188     assert( self );
00189     _Interpolator_Build( splineInterpolator, data );
00190     
00191     self->um = malloc( sizeof(double) * (self->regularMesh->nx-1) );
00192     memset( self->um, 0, sizeof( double ) * (self->regularMesh->nx-1) );
00193     
00194     self->un = malloc( sizeof(double) * (self->regularMesh->ny-1) );
00195     memset( self->un, 0, sizeof( double ) * (self->regularMesh->ny-1) );
00196     
00197     self->x1a = malloc( sizeof( float ) * self->regularMesh->nx );
00198     memset( self->x1a, 0, sizeof( float ) * self->regularMesh->nx );
00199     
00200     self->x2a = malloc( sizeof( float ) * self->regularMesh->ny );
00201     memset( self->x2a, 0, sizeof( float ) * self->regularMesh->ny );
00202     
00203     self->y2a = (float**)malloc( sizeof( float* ) * self->regularMesh->nx );
00204     self->y2a[0] = ( float* )malloc( sizeof( float ) * self->regularMesh->nx * self->regularMesh->ny );
00205 
00206     self->ya_t = malloc( sizeof( float ) * self->regularMesh->ny );
00207     memset( self->ya_t, 0, sizeof( float ) * self->regularMesh->ny );
00208     
00209     self->y2a_t = malloc( sizeof( float ) * self->regularMesh->ny );
00210     memset( self->y2a_t, 0, sizeof( float ) *self->regularMesh->ny );
00211     
00212     self->yytmp = malloc( sizeof( float ) * self->regularMesh->nx );
00213     memset( self->yytmp, 0, sizeof( float ) * self->regularMesh->nx );
00214     
00215     self->ytmp = malloc( sizeof( float ) * self->regularMesh->nx );
00216     memset( self->ytmp, 0, sizeof( float ) * self->regularMesh->nx );
00217 
00218     for( i=0; i<self->regularMesh->nx; i++ ){
00219         self->y2a[i] = self->y2a[0] + i * self->regularMesh->ny;
00220     }
00221 
00222     for( i=0; i<self->regularMesh->nx; i++ ){
00223         self->x1a[i] = self->regularMesh->regularNodes[i][0].x;
00224     }
00225     
00226     for( i=0; i<self->regularMesh->ny; i++ ){
00227         self->x2a[i] = self->regularMesh->regularNodes[0][i].y;
00228     }
00229 }
00230 
00231 void _SplineInterpolator_Initialise( void *splineInterpolator, void *data )
00232 {
00233 
00234 }
00235 
00236 void _SplineInterpolator_Execute( void *splineInterpolator, void *data )
00237 {
00238     
00239 }
00240 
00241 void _SplineInterpolator_Destroy( void *splineInterpolator, void *data )
00242 {
00243     
00244 }
00245 
00246 /* public functions */
00247 
00248 void SplineInterpolator_InterpolateFromGridToMesh( _Interpolator *interpolator,
00249         float **gridDeflectionArray,
00250         float *meshDeflectionArray )
00251 {
00252     int i = 0;
00253     float result;
00254     SurfaceMesh *mesh = NULL;
00255     SplineInterpolator *splineInterpolator;
00256 
00257     assert( interpolator );
00258     
00259     assert( gridDeflectionArray );
00260     assert( meshDeflectionArray );
00261 
00262     splineInterpolator = (SplineInterpolator*)interpolator;
00263     
00264     mesh = interpolator->mesh;
00265     
00266     memset( meshDeflectionArray, 0, sizeof( float ) * interpolator->mesh->numNodes );
00267     
00268     splie2( splineInterpolator, gridDeflectionArray );
00269     
00270     for( i=0; i<mesh->myLoad; i++ ){
00271         splin2( splineInterpolator,
00272                 gridDeflectionArray,
00273                 mesh->x[i], mesh->y[i], &result );
00274 
00275         meshDeflectionArray[mesh->id[i]-1] = result;
00276     }
00277 }
00278 
00279 /* private functions */
00280 
00281 /* The following spline interpolation routines have 
00282  * been taken from the Numerical Recipes for C++ book */
00283 
00284 void spline( int numElements, float *x, float *y, float yp1, float ypn, float *y2, double *u )
00285 {
00286     int i, k, n;
00287     double p, qn, sig, un;
00288 
00289     n = numElements;
00290 
00291     if( yp1 > 0.99e30 ){
00292         y2[0] = u[0] = 0.0f;
00293     }
00294     else{
00295         y2[0] = -0.5;
00296         u[0] = (3.0 / ( x[1]-x[0] )) * (( y[1]-y[0] )/( x[1]-x[0] )-yp1);
00297     }
00298 
00299     for( i=1; i<n-1; i++ ){
00300         sig = ( x[i]-x[i-1] ) / ( x[i+1]-x[i-1] );
00301         p = sig*y2[i-1] + 2.0f;
00302         y2[i] = ( sig - 1.0f ) / p;
00303 
00304         u[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
00305         u[i] = (6.0*u[i]/(x[i+1]-x[i-1]) - sig*u[i-1])/p;
00306     }
00307 
00308     if( ypn > 0.99e30 ){
00309         qn = un = 0.0f;
00310     }
00311     else{
00312         qn = 0.5;
00313         un = (3.0/(x[n-1]-x[n-2])) * (ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
00314     }
00315     
00316     y2[n-1] = (un-qn*u[n-2])/(qn*y2[n-2] + 1.0);
00317 
00318     for( k=n-2; k>=0; k-- ){
00319         y2[k] = y2[k]*y2[k+1] + u[k];
00320     }
00321 }
00322 
00323 void splint( int numElements, float *xa, float *ya, float *y2a, float x, float *y )
00324 {
00325     int k, n, klo, khi;
00326     double h, b, a;
00327     
00328     n = numElements;
00329     klo = 0;
00330     khi = n-1;
00331 
00332     while( khi-klo > 1 ){
00333         k = ( khi+klo ) >> 1;
00334 
00335         if( xa[k] > x ){
00336             khi = k;
00337         }
00338         else{
00339             klo = k;
00340         }
00341     }
00342 
00343     h = xa[khi] - xa[klo];
00344     if( h == 0.0 ){
00345         fprintf( stderr, "Bad input in spline interpolation routine..! Aborting..!\n" );
00346         assert( 0 );
00347     }
00348 
00349     a = ( xa[khi]-x ) / h;
00350     b = ( x - xa[klo] ) / h;
00351 
00352     *y = a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi])*(h*h)/6.0;
00353 }
00354 
00355 void splie2( SplineInterpolator *splInt, float **ya )
00356 {
00357     int m, n, j, k;
00358 
00359     assert( splInt );
00360     
00361     m = splInt->regularMesh->nx;
00362     n = splInt->regularMesh->ny;
00363 
00364     for( j=0; j<m; j++ ){
00365         for( k=0; k<n; k++ ){
00366             splInt->ya_t[k] = ya[j][k];
00367         }
00368 
00369         spline( n, splInt->x2a, splInt->ya_t, 1.0e30, 1.0e30, splInt->y2a_t, splInt->un );
00370 
00371         for( k=0; k<n; k++ ){
00372             splInt->y2a[j][k] = splInt->y2a_t[k];
00373         }
00374     }
00375 }
00376 
00377 void splin2( SplineInterpolator *splInt, float **ya, float x1, float x2, float *y )
00378 {
00379     int j, k, m, n;
00380     
00381     assert( splInt );
00382     
00383     m = splInt->regularMesh->nx;
00384     n = splInt->regularMesh->ny;
00385 
00386     for( j=0; j<m; j++ ){
00387         for( k=0; k<n; k++ ){
00388             splInt->ya_t[k] = ya[j][k];
00389             splInt->y2a_t[k] = splInt->y2a[j][k];
00390         }
00391         splint( n, splInt->x2a, splInt->ya_t, splInt->y2a_t, x2, &(splInt->yytmp[j]) );
00392     }
00393 
00394     spline( m, splInt->x1a, splInt->yytmp, 1.0e30, 1.0e30, splInt->ytmp, splInt->um );
00395     splint( m, splInt->x1a, splInt->yytmp, splInt->ytmp, x1, y );
00396 }
00397