VPAC - Computational Software Development

Main | SPModel | StGermain FrameWork |

IsostaticFlexure.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.
**
** $Id: IsostaticFlexure.c 225 2005-12-22 00:01:19Z AlanLo $
**
**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/


#include <mpi.h>
#include <StGermain/StGermain.h>
#include "SPModel/SPModel.h"
#include "flex2d.h"
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include "IsostaticFlexure.h"

/* Textual name of this module */
const Type SPModelIsostaticFlexure_Type = "SPModelIsostaticFlexure";

ExtensionInfo_Index SPModel_Simulation_ContextExtHandle;

void* _SPModelIsostaticFlexure_DefaultNew( Name name ) {
    return Codelet_New(
            SPModelIsostaticFlexure_Type,
            _SPModelIsostaticFlexure_DefaultNew,
            _SPModelIsostaticFlexure_Construct,
            _Codelet_Build,
            _Codelet_Initialise,
            _Codelet_Execute,
            _Codelet_Destroy,
            name );
}

void _SPModelIsostaticFlexure_Construct( void* component, Stg_ComponentFactory* data ) {    
    SPModel_Context*                context;

    context = (SPModel_Context*)Stg_ComponentFactory_ConstructByName( data, "context", SPModel_Context, True );
    
    assert( context );

    SPModel_Simulation_ContextExtHandle = ExtensionManager_GetHandle( context->extensionMgr, "Simulation" );

    _SPModelIsostaticFlexure_Init( context );
    
    Journal_DPrintf( context->debug, "In %s():\n", __func__ );
    ContextEP_Append( context, AbstractContext_EP_Solve, SPModelIsostaticFlexure_Solve );
}   

void _SPModelIsostaticFlexure_Init( void* _context ) {

}

float interpolateFunction ( int meshLocalIndex, void *args )
{
    float result = 0.0f;
    SPModel_Context *context = (SPModel_Context *) args;
    SPModelSimulationContextExtension *simulationExt = NULL;

    assert( context );

    simulationExt = (SPModelSimulationContextExtension*)ExtensionManager_Get
                ( context->extensionMgr, context, SPModel_Simulation_ContextExtHandle );

    assert( simulationExt );

    result = 0.0f;
    
        result = -simulationExt->erosion[meshLocalIndex] * context->bedRockDensity*context->gravity
                 +simulationExt->sedimentHistory[meshLocalIndex] * context->sedimentDensity*context->gravity
                 +simulationExt->uplift[meshLocalIndex] * context->bedRockDensity*context->gravity;
    
    return result;
}

void SPModelIsostaticFlexure_Solve( void *_context ) {
    double **loadg;
    double **deflec;
    double ***loadt;
    double **speq;
    SPModelSimulationContextExtension *simulationExt = NULL;
    int i, j;

    SPModel_Context *context = NULL;
    SurfaceMesh *globalMesh = NULL, *localMesh = NULL;
    SurfaceRegularMesh *regularMesh = NULL;
    _Interpolator *interpolator = NULL;
    
    context = (SPModel_Context*)_context;
    assert( context );
    
    if( ( context->timeStep % context->flexureInterval ) ){
        return;
    }

    simulationExt = (SPModelSimulationContextExtension*)ExtensionManager_Get
                    ( context->extensionMgr, context, SPModel_Simulation_ContextExtHandle );
    
    globalMesh = context->globalMesh;
    localMesh = context->localMesh;
    interpolator = context->interpolator;
    regularMesh = context->regularMesh;
    
    _Interpolator_InterpolateFromMeshToGrid( interpolator, interpolateFunction, context );
    
    /*Stg_Class_Print( interpolator, Journal_Register( Info_Type, "lota" ) );*/

    SurfaceRegularMesh_GatherData( regularMesh, interpolator );

    if( localMesh->rank == MASTER_PROC ){

        loadg = dmatrix(1,regularMesh->ny,1,regularMesh->nx);
        deflec = dmatrix(1,regularMesh->ny,1,regularMesh->nx);
        loadt = d3tensor(1,1,1,2*regularMesh->nx,1,2*regularMesh->ny);
        speq = dmatrix(1,1,1,2*2*regularMesh->nx);

        for( i=1; i<=regularMesh->nx; i++ ){
            for( j=1; j<=regularMesh->ny; j++ ){
                loadg[i][j] = regularMesh->dataArray[i-1][j-1];
            }
        }

        cosfilt( context->elasticPlateLength, regularMesh->nx, context->elasticPlateLength, 
                regularMesh->ny, loadg, context->flexuralRigidity,
                context->asthenosphereDensity, context->gravity, deflec, loadt, speq );

        for( i=1; i<=regularMesh->nx; i++ ){
            for( j=1; j<=regularMesh->nx; j++ ){
                regularMesh->dataArray[i-1][j-1] = deflec[i][j];
            }
        }
    
        free_dmatrix(loadg,1,regularMesh->ny,1,regularMesh->nx);
        free_dmatrix(deflec,1,regularMesh->ny,1,regularMesh->nx);
        free_d3tensor(loadt,1,1,1,2*regularMesh->nx,1,2*regularMesh->ny);
        free_dmatrix(speq,1,1,1,2*2*regularMesh->nx);
    }
    
    SurfaceRegularMesh_BroadcastData( regularMesh );
    
    _Interpolator_InterpolateFromGridToMesh( interpolator, regularMesh->dataArray, simulationExt->hiso );

    for( i=0; i<localMesh->numNodes; i++ ){
        simulationExt->hisotot[i] = -simulationExt->hisoPrev[i];
        simulationExt->hisotot[i] += simulationExt->hiso[i];
        simulationExt->hisoPrev[i] = simulationExt->hiso[i];
    }

    for( i=0; i<localMesh->myLoad; i++ ){
        localMesh->h[i] += simulationExt->hisotot[localMesh->id[i]-1];
    }
}

Index SPModelIsostaticFlexure_Register( PluginsManager* pluginsManager ) {
    return PluginsManager_Submit( pluginsManager, SPModelIsostaticFlexure_Type, "0", _SPModelIsostaticFlexure_DefaultNew );
}

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