VPAC - Computational Software Development

Main | SPModel | StGermain FrameWork |

SplineInterpolator.c

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
**
** Copyright (C), 2004, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street, Melbourne, 3053, Australia.
**
** Authors:
**  Ogar R. Widjaja, Computational Scientist, VPAC.
**  Raquibul Hassan, Software Engineer, VPAC. (raq@vpac.org)
**  Keith Hsuan, Computational Scientist, VPAC (keith@vpac.org)
**  William F. Appelbe, Director, VPAC. (bill@vpac.org)
**  Stevan M. Quenette, Senior Software Engineer, VPAC. (steve@vpac.org)
**  Patrick D. Sunter, Software Engineer, VPAC. (patrick@vpac.org)
**
** This file may be distributed under the terms of the VPAC Public License
** as defined by VPAC of Australia and appearing in the file
** LICENSE.VPL included in the packaging of this file.
**
** This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
** WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
**
*/
#include <mpi.h>
#include <StGermain/StGermain.h>
#include <Cascade/cascade.h>

#include "types.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include "Catchment.h"
#include "SurfaceMesh.h"
#include "SurfaceRegularMesh.h"
#include "Context.h"
#include "_Interpolator.h"
#include "SplineInterpolator.h"

/* Textual name of this class */
const Type SplineInterpolator_Type = "SplineInterpolator";

SplineInterpolator *SplineInterpolator_DefaultNew( Name name )
{
    return _SplineInterpolator_New( sizeof( SplineInterpolator ),
                                SplineInterpolator_Type,
                                _SplineInterpolator_Delete,
                                _SplineInterpolator_Print,
                                NULL,
                                (void*)SplineInterpolator_DefaultNew,
                                _SplineInterpolator_Construct,
                                _SplineInterpolator_Build,
                                _SplineInterpolator_Initialise,
                                _SplineInterpolator_Execute,
                                _SplineInterpolator_Destroy,
                                name,
                                False,
                                NULL,
                                NULL,
                                NULL );
    
}

SplineInterpolator *SplineInterpolator_New( Name name, Dictionary *dictionary, SurfaceMesh *mesh, SurfaceRegularMesh *regularMesh )
{
    return _SplineInterpolator_New( sizeof( SplineInterpolator ),
                                SplineInterpolator_Type,
                                _SplineInterpolator_Delete,
                                _SplineInterpolator_Print,
                                NULL,
                                (void*)SplineInterpolator_DefaultNew,
                                _SplineInterpolator_Construct,
                                _SplineInterpolator_Build,
                                _SplineInterpolator_Initialise,
                                _SplineInterpolator_Execute,
                                _SplineInterpolator_Destroy,
                                name,
                                True,
                                dictionary,
                                mesh,
                                regularMesh );
    
}

SplineInterpolator          *_SplineInterpolator_New( SizeT                         _sizeOfSelf,
                                    Type                            type,
                                    Stg_Class_DeleteFunction*               _delete,
                                    Stg_Class_PrintFunction*                _print,
                                    Stg_Class_CopyFunction*             _copy, 
                                    Stg_Component_DefaultConstructorFunction*   _defaultConstructor,
                                    Stg_Component_ConstructFunction*            _construct,
                                    Stg_Component_BuildFunction*        _build,
                                    Stg_Component_InitialiseFunction*       _initialise,
                                    Stg_Component_ExecuteFunction*      _execute,
                                    Stg_Component_DestroyFunction*      _destroy,
                                    Name                            name,
                                    Bool                            initFlag,
                                    Dictionary                      *dictionary,
                                    SurfaceMesh                     *mesh,
                                    SurfaceRegularMesh              *regularMesh )
{
    SplineInterpolator* self;
    
    /* Allocate memory */
    assert( _sizeOfSelf >= sizeof(SplineInterpolator) );
    self = (SplineInterpolator*)_Interpolator_New( _sizeOfSelf, type, _delete, _print, _copy, _defaultConstructor, _construct, _build, 
            _initialise, _execute, _destroy, name, NON_GLOBAL, dictionary, mesh, regularMesh, SplineInterpolator_InterpolateFromGridToMesh );

    if( initFlag ){
        _SplineInterpolator_Init( self );
    }

    return self;
}

void _SplineInterpolator_Init( SplineInterpolator *self )
{
    _Interpolator_Init( (_Interpolator*)self );
}

void _SplineInterpolator_Print( void *splineInterpolator, Stream *stream )
{
    int numNodes;
    
    SplineInterpolator *self = (SplineInterpolator*) splineInterpolator;

    assert( self );
    assert( stream );
    
    if( self->mesh->rank == MASTER_PROC ){
        numNodes = self->mesh->numNodes;
    }
    else{
        numNodes = self->mesh->myLoad;
    }

    /* print parent */
    _Stg_Component_Print( (void*) self, stream );
    
    _Interpolator_Print( (void*) self, stream );

    /* general info */
    Journal_Printf( stream, "SplineInterpolator (ptr): (%p)\n", self );

    /* Virtual Info */

    /* SplineInterpolator Info */
    Journal_Printf( stream, "Rank - %d\n", self->mesh->rank );
}

void _SplineInterpolator_Delete( void *splineInterpolator )
{
    SplineInterpolator *self = (SplineInterpolator*)splineInterpolator;

    assert( self );

    free( self->um );
    free( self->un );
    free( self->x1a );
    free( self->x2a );
    free( self->y2a[0] );
    free( self->y2a );
    free( self->ytmp );
    free( self->yytmp );
    free( self->y2a_t );
    free( self->ya_t );
    
    _Interpolator_Delete( self );
}

void _SplineInterpolator_Construct( void *splineInterpolator, Stg_ComponentFactory *cf )
{
}

void _SplineInterpolator_Build( void *splineInterpolator, void *data )
{
    int i;
    SplineInterpolator *self = (SplineInterpolator*)splineInterpolator;

    assert( self );
    _Interpolator_Build( splineInterpolator, data );
    
    self->um = malloc( sizeof(double) * (self->regularMesh->nx-1) );
    memset( self->um, 0, sizeof( double ) * (self->regularMesh->nx-1) );
    
    self->un = malloc( sizeof(double) * (self->regularMesh->ny-1) );
    memset( self->un, 0, sizeof( double ) * (self->regularMesh->ny-1) );
    
    self->x1a = malloc( sizeof( float ) * self->regularMesh->nx );
    memset( self->x1a, 0, sizeof( float ) * self->regularMesh->nx );
    
    self->x2a = malloc( sizeof( float ) * self->regularMesh->ny );
    memset( self->x2a, 0, sizeof( float ) * self->regularMesh->ny );
    
    self->y2a = (float**)malloc( sizeof( float* ) * self->regularMesh->nx );
    self->y2a[0] = ( float* )malloc( sizeof( float ) * self->regularMesh->nx * self->regularMesh->ny );

    self->ya_t = malloc( sizeof( float ) * self->regularMesh->ny );
    memset( self->ya_t, 0, sizeof( float ) * self->regularMesh->ny );
    
    self->y2a_t = malloc( sizeof( float ) * self->regularMesh->ny );
    memset( self->y2a_t, 0, sizeof( float ) *self->regularMesh->ny );
    
    self->yytmp = malloc( sizeof( float ) * self->regularMesh->nx );
    memset( self->yytmp, 0, sizeof( float ) * self->regularMesh->nx );
    
    self->ytmp = malloc( sizeof( float ) * self->regularMesh->nx );
    memset( self->ytmp, 0, sizeof( float ) * self->regularMesh->nx );

    for( i=0; i<self->regularMesh->nx; i++ ){
        self->y2a[i] = self->y2a[0] + i * self->regularMesh->ny;
    }

    for( i=0; i<self->regularMesh->nx; i++ ){
        self->x1a[i] = self->regularMesh->regularNodes[i][0].x;
    }
    
    for( i=0; i<self->regularMesh->ny; i++ ){
        self->x2a[i] = self->regularMesh->regularNodes[0][i].y;
    }
}

void _SplineInterpolator_Initialise( void *splineInterpolator, void *data )
{

}

void _SplineInterpolator_Execute( void *splineInterpolator, void *data )
{
    
}

void _SplineInterpolator_Destroy( void *splineInterpolator, void *data )
{
    
}

/* public functions */

void SplineInterpolator_InterpolateFromGridToMesh( _Interpolator *interpolator,
        float **gridDeflectionArray,
        float *meshDeflectionArray )
{
    int i = 0;
    float result;
    SurfaceMesh *mesh = NULL;
    SplineInterpolator *splineInterpolator;

    assert( interpolator );
    
    assert( gridDeflectionArray );
    assert( meshDeflectionArray );

    splineInterpolator = (SplineInterpolator*)interpolator;
    
    mesh = interpolator->mesh;
    
    memset( meshDeflectionArray, 0, sizeof( float ) * interpolator->mesh->numNodes );
    
    splie2( splineInterpolator, gridDeflectionArray );
    
    for( i=0; i<mesh->myLoad; i++ ){
        splin2( splineInterpolator,
                gridDeflectionArray,
                mesh->x[i], mesh->y[i], &result );

        meshDeflectionArray[mesh->id[i]-1] = result;
    }
}

/* private functions */

/* The following spline interpolation routines have 
 * been taken from the Numerical Recipes for C++ book */

void spline( int numElements, float *x, float *y, float yp1, float ypn, float *y2, double *u )
{
    int i, k, n;
    double p, qn, sig, un;

    n = numElements;

    if( yp1 > 0.99e30 ){
        y2[0] = u[0] = 0.0f;
    }
    else{
        y2[0] = -0.5;
        u[0] = (3.0 / ( x[1]-x[0] )) * (( y[1]-y[0] )/( x[1]-x[0] )-yp1);
    }

    for( i=1; i<n-1; i++ ){
        sig = ( x[i]-x[i-1] ) / ( x[i+1]-x[i-1] );
        p = sig*y2[i-1] + 2.0f;
        y2[i] = ( sig - 1.0f ) / p;

        u[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
        u[i] = (6.0*u[i]/(x[i+1]-x[i-1]) - sig*u[i-1])/p;
    }

    if( ypn > 0.99e30 ){
        qn = un = 0.0f;
    }
    else{
        qn = 0.5;
        un = (3.0/(x[n-1]-x[n-2])) * (ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
    }
    
    y2[n-1] = (un-qn*u[n-2])/(qn*y2[n-2] + 1.0);

    for( k=n-2; k>=0; k-- ){
        y2[k] = y2[k]*y2[k+1] + u[k];
    }
}

void splint( int numElements, float *xa, float *ya, float *y2a, float x, float *y )
{
    int k, n, klo, khi;
    double h, b, a;
    
    n = numElements;
    klo = 0;
    khi = n-1;

    while( khi-klo > 1 ){
        k = ( khi+klo ) >> 1;

        if( xa[k] > x ){
            khi = k;
        }
        else{
            klo = k;
        }
    }

    h = xa[khi] - xa[klo];
    if( h == 0.0 ){
        fprintf( stderr, "Bad input in spline interpolation routine..! Aborting..!\n" );
        assert( 0 );
    }

    a = ( xa[khi]-x ) / h;
    b = ( x - xa[klo] ) / h;

    *y = a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi])*(h*h)/6.0;
}

void splie2( SplineInterpolator *splInt, float **ya )
{
    int m, n, j, k;

    assert( splInt );
    
    m = splInt->regularMesh->nx;
    n = splInt->regularMesh->ny;

    for( j=0; j<m; j++ ){
        for( k=0; k<n; k++ ){
            splInt->ya_t[k] = ya[j][k];
        }

        spline( n, splInt->x2a, splInt->ya_t, 1.0e30, 1.0e30, splInt->y2a_t, splInt->un );

        for( k=0; k<n; k++ ){
            splInt->y2a[j][k] = splInt->y2a_t[k];
        }
    }
}

void splin2( SplineInterpolator *splInt, float **ya, float x1, float x2, float *y )
{
    int j, k, m, n;
    
    assert( splInt );
    
    m = splInt->regularMesh->nx;
    n = splInt->regularMesh->ny;

    for( j=0; j<m; j++ ){
        for( k=0; k<n; k++ ){
            splInt->ya_t[k] = ya[j][k];
            splInt->y2a_t[k] = splInt->y2a[j][k];
        }
        splint( n, splInt->x2a, splInt->ya_t, splInt->y2a_t, x2, &(splInt->yytmp[j]) );
    }

    spline( m, splInt->x1a, splInt->yytmp, 1.0e30, 1.0e30, splInt->ytmp, splInt->um );
    splint( m, splInt->x1a, splInt->yytmp, splInt->ytmp, x1, y );
}

SPModel documentation generated by doxygen at Fri Mar 24 01:41:45 2006