SPModel: plugins/faultModel/FaultModel.c Source 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: FaultModel.c 225 2005-12-22 00:01:19Z AlanLo $
00021 **
00022 **~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
00023 
00024 
00025 #include <mpi.h>
00026 #include <StGermain/StGermain.h>
00027 #include "SPModel/SPModel.h"
00028 #include <stdlib.h>
00029 #include <string.h>
00030 #include <assert.h>
00031 
00032 
00033 #include "FaultModel.h"
00034 
00035 /* Textual name of this module */
00036 const Type SPModelFaultModel_Type = "SPModelFaultModel";
00037 ExtensionInfo_Index SPModel_FaultModel_ContextExtHandle;
00038 ExtensionInfo_Index SPModel_Simulation_ContextExtHandle;
00039 
00040 void* _SPModelFaultModel_DefaultNew( Name name ) {
00041     return Codelet_New(
00042             SPModelFaultModel_Type,
00043             _SPModelFaultModel_DefaultNew,
00044             _SPModelFaultModel_Construct,
00045             _Codelet_Build,
00046             _Codelet_Initialise,
00047             _Codelet_Execute,
00048             _Codelet_Destroy,
00049             name );
00050 }
00051 
00052 void _SPModelFaultModel_Construct( void* component, Stg_ComponentFactory* data ) {  
00053     SPModel_Context*                context;
00054 
00055     context = (SPModel_Context*)Stg_ComponentFactory_ConstructByName( data, "context", SPModel_Context, True );
00056     
00057     SPModel_FaultModel_ContextExtHandle = ExtensionManager_Add( context->extensionMgr, "FaultModel",
00058                                             sizeof(SPModelFaultModelContextExtension) );
00059     
00060     SPModel_Simulation_ContextExtHandle = ExtensionManager_GetHandle( context->extensionMgr, "Simulation" );
00061 
00062     _SPModelFaultModel_Init( context );
00063     
00064     Journal_DPrintf( context->debug, "In %s():\n", __func__ );
00065     ContextEP_Append( context, AbstractContext_EP_Solve, SPModelFaultModel_Solve );
00066 }   
00067 
00068 #define EPSILON 1e-06
00069 int colinearityCheck( Coord a, Coord b, Coord c )
00070 {
00071     Coord ab;
00072     Coord ac;
00073 
00074     StGermain_VectorSubtraction( ab, b, a, 3 );
00075     StGermain_VectorSubtraction( ac, c, a, 3 );
00076 
00077     StGermain_VectorNormalise( ab, 3 );
00078     StGermain_VectorNormalise( ac, 3 );
00079     
00080     /*printf( "angle  %lf\n", StGermain_VectorDotProduct( ab, ac, 3 ) );*/
00081     if( (1.0f - fabs(StGermain_VectorDotProduct( ab, ac, 3 ))) < EPSILON ){
00082         return 1;
00083     }
00084     else{
00085         return 0;
00086     }
00087 }
00088 
00089 void _SPModelFaultModel_Init( void* _context ) {
00090     SPModel_Context *context = NULL;
00091     SPModelFaultModelContextExtension *faultModel = NULL;
00092     Dictionary *dict = NULL, *subDict = NULL;
00093     Dictionary_Entry    *subDictEntry = NULL;
00094     FILE *f;
00095     char *fileName = NULL;
00096     int i = 0, j = 0;
00097     char buffer[512] = { 0 };
00098     Coord midChord1, midChord2;
00099     Coord direcChord1, direcChord2;
00100     Coord cp1, cp2;
00101     Coord rhs, lhs;
00102     Coord nodeCoord;
00103     double temp = 0.0f, lambda = 0.0f;
00104     SurfaceMesh *mesh = NULL;
00105 
00106     context = (SPModel_Context*)_context;
00107     assert( context );
00108     
00109     faultModel = (SPModelFaultModelContextExtension*)ExtensionManager_Get
00110                     ( context->extensionMgr, context, SPModel_FaultModel_ContextExtHandle );
00111     assert( faultModel );
00112 
00113     mesh = context->localMesh;
00114     assert( mesh );
00115     
00116     dict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( context->dictionary, "faultModel" ) );
00117 
00118     if( dict ){
00119 
00120         faultModel->numFaultModel = Dictionary_GetCount( dict );
00121 
00122         faultModel->coordType = Memory_Alloc_Array_Unnamed( CoordinateType, faultModel->numFaultModel );
00123         faultModel->arcPoints = Memory_Alloc_2DArray_Unnamed( Coord, faultModel->numFaultModel, 4);
00124         faultModel->arcs = Memory_Alloc_2DArray_Unnamed( Arc, faultModel->numFaultModel, 2 );
00125         faultModel->circleCentrePoints = Memory_Alloc_2DArray_Unnamed( Coord, faultModel->numFaultModel, 2 );
00126         faultModel->radiusSq = Memory_Alloc_2DArray_Unnamed( double, faultModel->numFaultModel, 2 );
00127         faultModel->faultBounds = Memory_Alloc_Array_Unnamed( int, faultModel->numFaultModel );
00128 
00129         for( i=0; i<faultModel->numFaultModel; i++ ){
00130             
00131             subDictEntry = dict->entryPtr[ i ];
00132             subDict = Dictionary_Entry_Value_AsDictionary( subDictEntry->value );
00133 
00134             assert( subDict );
00135 
00136             fileName = Dictionary_GetString( subDict, "faultModelBoundaryFile" );
00137             f = fopen( fileName, "r+" );
00138 
00139             if( f ){
00140                 
00141                 memset( buffer, 0, sizeof( buffer ) );
00142                 
00143                 fgets( buffer, sizeof( buffer ), f );
00144                 if( strcasecmp(buffer, "NORMALISED\n") == 0 ){
00145                     faultModel->coordType[i] = NORMALISED;
00146                 }
00147                 else if( strcasecmp(buffer, "NATURAL\n") == 0 ){
00148                     faultModel->coordType[i] = NATURAL;
00149                 }
00150                 else{
00151                     Journal_Printf( context->info, "Error in file `%s`\n", fileName );
00152                     exit( EXIT_FAILURE );
00153                 }
00154 
00155                 fscanf( f, "%d", &(faultModel->faultBounds[i]) );
00156                 
00157                 for( j=0; j<4; j++ ){
00158                     fscanf( f, "%lf %lf", &(faultModel->arcPoints[i][j][0]), &(faultModel->arcPoints[i][j][1]) );
00159                 }
00160 
00161                 for( j=0; j<2; j++ ){
00162                     fscanf( f, "%d %d %d", &(faultModel->arcs[i][j].index[0]), &(faultModel->arcs[i][j].index[1]),
00163                             &(faultModel->arcs[i][j].index[2]) );
00164 
00165                     faultModel->arcs[i][j].index[0] -= 1;
00166                     faultModel->arcs[i][j].index[1] -= 1;
00167                     faultModel->arcs[i][j].index[2] -= 1;
00168                 }
00169                 fclose( f );
00170             }
00171             else{
00172                 Journal_Printf( context->info, "Error opening file `%s`..!\n", fileName );
00173                 exit( EXIT_FAILURE );
00174             }
00175         }
00176     
00177 #if 0
00178         printf( (faultModel->coordType==NORMALISED)?"normalised\n":"natural\n" );
00179         printf( "%d\n", faultModel->numFaultModel );
00180         for( i=0; i<faultModel->numFaultModel; i++ ){
00181             for( j=0; j<4; j++ ){
00182                 printf( "%lf %lf\n", faultModel->arcPoints[i][j][0], faultModel->arcPoints[i][j][1] );
00183             }
00184             printf( "%d %d %d\n", faultModel->arcs[i][0].index[0], faultModel->arcs[i][0].index[1],
00185                     faultModel->arcs[i][0].index[2] );
00186             printf( "%d %d %d\n", faultModel->arcs[i][1].index[0], faultModel->arcs[i][1].index[1],
00187                     faultModel->arcs[i][1].index[2] );
00188             printf( "\n\n" );
00189         }
00190 #endif
00191 
00192         /*colinearity check*/
00193         for( i=0; i<faultModel->numFaultModel; i++ ){
00194             for( j=0; j<2; j++ ){
00195                 if( colinearityCheck( faultModel->arcPoints[i][faultModel->arcs[i][j].index[0]],
00196                             faultModel->arcPoints[i][faultModel->arcs[i][j].index[1]],
00197                             faultModel->arcPoints[i][faultModel->arcs[i][j].index[2]]) ){
00198                     Journal_Printf( context->info, "Arc points for faultModel %d are colinear..!\n", i+1 );
00199                     exit( EXIT_FAILURE );
00200                 }
00201             }
00202         }
00203         
00204         /*TODO verify natural/normalised coordinates*/
00205         /*for( i=0; i<faultModel->numFaultModel; i++ ){
00206             for( j=0; j<4; j++ ){
00207                 if( faultModel->coordType == NORMALISED ){
00208                     for( k=0; k<3; k++ ){
00209                         if( faultModel->arcPoints[i][j][k] > 1.0f ){
00210                             Journal_Printf( context->info, 
00211                             "Arc Coordinates for faultModel %d not in normalised coordinates..!\n", i+1 );
00212                             exit( EXIT_FAILURE );
00213                         }
00214                     }
00215                 }
00216             }
00217         }*/
00218             
00219         faultModel->nodeIndices = Memory_Alloc_2DArray_Unnamed( int, faultModel->numFaultModel, mesh->myLoad );
00220         faultModel->distFromBoundary = Memory_Alloc_2DArray_Unnamed( double, faultModel->numFaultModel, mesh->myLoad );
00221 
00222         for( i=0; i<faultModel->numFaultModel; i++ ){
00223 
00224             for( j=0; j<2; j++ ){
00225             
00226                 memset( midChord1, 0, sizeof(Coord) );
00227                 midChord1[0] = (faultModel->arcPoints[i][faultModel->arcs[i][j].index[0]][0] +
00228                                 faultModel->arcPoints[i][faultModel->arcs[i][j].index[1]][0]) / 2.0f;
00229                 midChord1[1] = (faultModel->arcPoints[i][faultModel->arcs[i][j].index[0]][1] +
00230                                 faultModel->arcPoints[i][faultModel->arcs[i][j].index[1]][1]) / 2.0f;
00231             
00232                 memset( midChord2, 0, sizeof(Coord) );
00233                 midChord2[0] = (faultModel->arcPoints[i][faultModel->arcs[i][j].index[0]][0] +
00234                                 faultModel->arcPoints[i][faultModel->arcs[i][j].index[2]][0]) / 2.0f;
00235                 midChord2[1] = (faultModel->arcPoints[i][faultModel->arcs[i][j].index[0]][1] +
00236                                 faultModel->arcPoints[i][faultModel->arcs[i][j].index[2]][1]) / 2.0f;
00237 
00238                 memset( direcChord1, 0, sizeof( Coord ) );
00239                 direcChord1[0] = (faultModel->arcPoints[i][faultModel->arcs[i][j].index[1]][0] - 
00240                                     faultModel->arcPoints[i][faultModel->arcs[i][j].index[0]][0] );
00241                 direcChord1[1] = (faultModel->arcPoints[i][faultModel->arcs[i][j].index[1]][1] - 
00242                                     faultModel->arcPoints[i][faultModel->arcs[i][j].index[0]][1] );
00243                 
00244                 memset( direcChord2, 0, sizeof( Coord ) );
00245                 direcChord2[0] = (faultModel->arcPoints[i][faultModel->arcs[i][j].index[2]][0] -
00246                                     faultModel->arcPoints[i][faultModel->arcs[i][j].index[0]][0] );
00247                 direcChord2[1] = (faultModel->arcPoints[i][faultModel->arcs[i][j].index[2]][1] -
00248                                     faultModel->arcPoints[i][faultModel->arcs[i][j].index[0]][1] );
00249 
00250                 StGermain_VectorNormalise( direcChord1, 3 );
00251                 StGermain_VectorNormalise( direcChord2, 3 );
00252 
00253                 /* taking the perpendiculars */
00254                 temp = direcChord1[0];
00255                 direcChord1[0] = -direcChord1[1];
00256                 direcChord1[1] = temp;
00257                 
00258                 temp = direcChord2[0];
00259                 direcChord2[0] = -direcChord2[1];
00260                 direcChord2[1] = temp;
00261 
00262                 StGermain_VectorCrossProduct( cp1, midChord1, direcChord2 );
00263                 StGermain_VectorCrossProduct( cp2, midChord2, direcChord2 );
00264                 StGermain_VectorSubtraction( rhs, cp2, cp1, 3 );
00265                 
00266                 StGermain_VectorCrossProduct( cp2, direcChord1, direcChord2 );
00267 
00268                 memcpy( lhs, cp2, sizeof( Coord ) );
00269                 lambda = rhs[2] / lhs[2];
00270 
00271                 faultModel->circleCentrePoints[i][j][0] = midChord1[0] + direcChord1[0]*lambda;
00272                 faultModel->circleCentrePoints[i][j][1] = midChord1[1] + direcChord1[1]*lambda;
00273                 faultModel->circleCentrePoints[i][j][2] = 0;
00274 
00275                 
00276                 faultModel->radiusSq[i][j] = pow( StGermain_DistanceBetweenPoints( 
00277                             faultModel->arcPoints[i][faultModel->arcs[i][j].index[0]], 
00278                             faultModel->circleCentrePoints[i][j], 3), 2.0f);
00279                 /*printf( "centre x %lf, y %lf\n", faultModel->circleCentrePoints[i][j][0], faultModel->circleCentrePoints[i][j][1] );
00280                 printf( "%lf\n", StGermain_DistanceBetweenPoints( faultModel->arcPoints[i][faultModel->arcs[i][j].index[0]], faultModel->circleCentrePoints[i][j], 3) );
00281                 printf( "%lf\n", StGermain_DistanceBetweenPoints( faultModel->arcPoints[i][faultModel->arcs[i][j].index[1]], faultModel->circleCentrePoints[i][j], 3) );
00282                 printf( "%lf\n", StGermain_DistanceBetweenPoints( faultModel->arcPoints[i][faultModel->arcs[i][j].index[2]], faultModel->circleCentrePoints[i][j], 3) );
00283                 printf( "\n\n" );*/
00284             }
00285         }
00286 
00287         for( i=0; i<faultModel->numFaultModel; i++ ){
00288             for( j=0; j<2; j++ ){
00289                 if( faultModel->coordType[i] == NORMALISED ){
00290                     faultModel->circleCentrePoints[i][j][0] *= mesh->sideX;
00291                     faultModel->circleCentrePoints[i][j][1] *= mesh->sideY;
00292 
00293                     faultModel->radiusSq[i][j] *= (mesh->sideX * mesh->sideY);
00294                 }
00295             }
00296 
00297             switch( faultModel->faultBounds[i] ){
00298                 
00299                 case 1:
00300                 {
00301                     for( j=0; j<mesh->myLoad; j++ ){
00302                         temp = pow( mesh->x[j] - faultModel->circleCentrePoints[i][0][0], 2.0f ) +
00303                                 pow( mesh->y[j] - faultModel->circleCentrePoints[i][0][1], 2.0f );
00304                         
00305                         faultModel->nodeIndices[i][j] = UNDEFINED;
00306                         if( temp > faultModel->radiusSq[i][0] ){ /* inside circle 1*/
00307                             temp = pow( mesh->x[j] - faultModel->circleCentrePoints[i][1][0], 2.0f ) +
00308                                     pow( mesh->y[j] - faultModel->circleCentrePoints[i][1][1], 2.0f );
00309 
00310                             if( temp < faultModel->radiusSq[i][1] ){
00311                                 faultModel->nodeIndices[i][j] = mesh->id[j];
00312                             }
00313                         }
00314                     }
00315                 }
00316                 break;
00317                 
00318                 case 2:
00319                 {
00320                     for( j=0; j<mesh->myLoad; j++ ){
00321                         temp = pow( mesh->x[j] - faultModel->circleCentrePoints[i][0][0], 2.0f ) +
00322                                 pow( mesh->y[j] - faultModel->circleCentrePoints[i][0][1], 2.0f );
00323                         
00324                         faultModel->nodeIndices[i][j] = UNDEFINED;
00325                         if( temp < faultModel->radiusSq[i][0] ){ /* inside circle 1*/
00326                             temp = pow( mesh->x[j] - faultModel->circleCentrePoints[i][1][0], 2.0f ) +
00327                                     pow( mesh->y[j] - faultModel->circleCentrePoints[i][1][1], 2.0f );
00328 
00329                             if( temp < faultModel->radiusSq[i][1] ){
00330                                 faultModel->nodeIndices[i][j] = mesh->id[j];
00331                             }
00332                         }
00333                     }
00334                 }
00335                 break;
00336                 
00337                 default:
00338                 {
00339                 }
00340                 break;
00341             }
00342 
00343             /*check if circle 2 is entirely within circle 1*/
00344             if( StGermain_DistanceBetweenPoints( faultModel->circleCentrePoints[i][0],
00345                 faultModel->circleCentrePoints[i][1], 3 ) + faultModel->radiusSq[i][1] <= 
00346                 sqrt( faultModel->radiusSq[i][0] ) ){
00347 
00348                 Journal_Printf( context->info, "Circle 2 inside circle 1 in faultModel %d..!\n", i+1 );
00349                 exit( EXIT_FAILURE );
00350             }
00351             
00352             for( j=0; j<mesh->myLoad; j++ ){
00353                 if( faultModel->nodeIndices[i][j] != UNDEFINED ){
00354                     nodeCoord[0] = mesh->x[mesh->mapGlobalToLocal[faultModel->nodeIndices[i][j]-1]];
00355                     nodeCoord[1] = mesh->y[mesh->mapGlobalToLocal[faultModel->nodeIndices[i][j]-1]];
00356                     nodeCoord[2] = 0;
00357                     faultModel->distFromBoundary[i][j] = (sqrt(faultModel->radiusSq[i][1]) - 
00358                         StGermain_DistanceBetweenPoints( nodeCoord, faultModel->circleCentrePoints[i][1], 3 )) /
00359                         ( sqrt(faultModel->radiusSq[i][0]) - (sqrt(faultModel->radiusSq[i][0]) - sqrt(faultModel->radiusSq[i][1])) );
00360                 }
00361             }
00362         }
00363 
00364         faultModel->fmUpliftRate = (ParameterTimeSeries**) malloc( sizeof( void* ) * faultModel->numFaultModel );
00365 
00366         for( i=0; i<faultModel->numFaultModel; i++ ){
00367             subDictEntry = dict->entryPtr[ i ];
00368             subDict = Dictionary_Entry_Value_AsDictionary( subDictEntry->value );
00369 
00370             assert( subDict );
00371 
00372             faultModel->fmUpliftRate[i] = ParameterTimeSeries_New( "faultModelUpliftRate", 0.0f, subDict );
00373             Stg_Component_Build( faultModel->fmUpliftRate[i], NULL, False );
00374         }
00375     }
00376     else{
00377         Journal_Printf( context->info, "faultModel struct not found in xml file..!\n" );
00378         exit( EXIT_FAILURE );
00379     }
00380 }
00381 
00382 void SPModelFaultModel_RefreshNodeIndices( void *_context )
00383 {
00384     SPModel_Context *context = (SPModel_Context*)_context;
00385     SurfaceMesh *mesh = NULL;
00386     SPModelFaultModelContextExtension *faultModel = NULL;
00387     int i = 0, j = 0;
00388     float temp = 0;
00389     Coord nodeCoord;
00390 
00391     assert( context );
00392     
00393     faultModel = (SPModelFaultModelContextExtension*)ExtensionManager_Get
00394                     ( context->extensionMgr, context, SPModel_FaultModel_ContextExtHandle );
00395     assert( faultModel );
00396     
00397     mesh = context->localMesh;
00398     assert( mesh );
00399     
00400     Memory_Free( faultModel->nodeIndices );
00401     Memory_Free( faultModel->distFromBoundary );
00402 
00403     faultModel->nodeIndices = Memory_Alloc_2DArray_Unnamed( int, faultModel->numFaultModel, mesh->myLoad );
00404     faultModel->distFromBoundary = Memory_Alloc_2DArray_Unnamed( double, faultModel->numFaultModel, mesh->myLoad );
00405     
00406     for( i=0; i<faultModel->numFaultModel; i++ ){
00407         for( j=0; j<mesh->myLoad; j++ ){
00408             faultModel->nodeIndices[i][j] = UNDEFINED;
00409         }
00410     }
00411     
00412     for( i=0; i<faultModel->numFaultModel; i++ ){
00413         switch( faultModel->faultBounds[i] ){
00414                 
00415             case 1:
00416             {
00417                 for( j=0; j<mesh->myLoad; j++ ){
00418                     temp = pow( mesh->x[j] - faultModel->circleCentrePoints[i][0][0], 2.0f ) +
00419                             pow( mesh->y[j] - faultModel->circleCentrePoints[i][0][1], 2.0f );
00420                     
00421                     faultModel->nodeIndices[i][j] = UNDEFINED;
00422                     if( temp > faultModel->radiusSq[i][0] ){ /* inside circle 1*/
00423                         temp = pow( mesh->x[j] - faultModel->circleCentrePoints[i][1][0], 2.0f ) +
00424                                 pow( mesh->y[j] - faultModel->circleCentrePoints[i][1][1], 2.0f );
00425 
00426                         if( temp < faultModel->radiusSq[i][1] ){
00427                             faultModel->nodeIndices[i][j] = mesh->id[j];
00428                         }
00429                     }
00430                 }
00431             }
00432             break;
00433             
00434             case 2:
00435             {
00436                 for( j=0; j<mesh->myLoad; j++ ){
00437                     temp = pow( mesh->x[j] - faultModel->circleCentrePoints[i][0][0], 2.0f ) +
00438                             pow( mesh->y[j] - faultModel->circleCentrePoints[i][0][1], 2.0f );
00439                     
00440                     faultModel->nodeIndices[i][j] = UNDEFINED;
00441                     if( temp < faultModel->radiusSq[i][0] ){ /* inside circle 1*/
00442                         temp = pow( mesh->x[j] - faultModel->circleCentrePoints[i][1][0], 2.0f ) +
00443                                 pow( mesh->y[j] - faultModel->circleCentrePoints[i][1][1], 2.0f );
00444 
00445                         if( temp < faultModel->radiusSq[i][1] ){
00446                             faultModel->nodeIndices[i][j] = mesh->id[j];
00447                         }
00448                     }
00449                 }
00450             }
00451             break;
00452             
00453             default:
00454             {
00455             }
00456             break;
00457         }
00458         
00459         for( j=0; j<mesh->myLoad; j++ ){
00460             if( faultModel->nodeIndices[i][j] != UNDEFINED ){
00461                 nodeCoord[0] = mesh->x[mesh->mapGlobalToLocal[faultModel->nodeIndices[i][j]-1]];
00462                 nodeCoord[1] = mesh->y[mesh->mapGlobalToLocal[faultModel->nodeIndices[i][j]-1]];
00463                 nodeCoord[2] = 0;
00464                 faultModel->distFromBoundary[i][j] = (sqrt(faultModel->radiusSq[i][1]) - 
00465                     StGermain_DistanceBetweenPoints( nodeCoord, faultModel->circleCentrePoints[i][1], 3 )) /
00466                     ( sqrt(faultModel->radiusSq[i][0]) - (sqrt(faultModel->radiusSq[i][0]) - sqrt(faultModel->radiusSq[i][1])) );
00467             }
00468         }
00469     }
00470 }
00471 
00472 void SPModelFaultModel_Solve( void *_context ) {
00473     SPModel_Context *context = (SPModel_Context*)_context;
00474     SurfaceMesh *mesh = NULL;
00475     SPModelFaultModelContextExtension *faultModel = NULL;
00476     SPModelSimulationContextExtension *simulationExt = NULL;
00477     int i = 0, j = 0, localIndex = 0;
00478 
00479     assert( context );
00480     
00481     faultModel = (SPModelFaultModelContextExtension*)ExtensionManager_Get
00482                     ( context->extensionMgr, context, SPModel_FaultModel_ContextExtHandle );
00483     simulationExt = (SPModelSimulationContextExtension*)ExtensionManager_Get
00484                     ( context->extensionMgr, context, SPModel_Simulation_ContextExtHandle );
00485     
00486     assert( faultModel );
00487     assert( simulationExt );
00488 
00489     mesh = context->localMesh;
00490     assert( mesh );
00491 
00492 
00493     if( context->redistributeNodesFlag ){
00494         SPModelFaultModel_RefreshNodeIndices( _context );
00495         
00496     }
00497     
00498     for( i=0; i<faultModel->numFaultModel; i++ ){
00499         ParameterTimeSeries_Interpolate( faultModel->fmUpliftRate[i], context );
00500         for( j=0; j<mesh->myLoad; j++ ){
00501             if( faultModel->nodeIndices[i][j] != UNDEFINED ){
00502                 localIndex = mesh->mapGlobalToLocal[faultModel->nodeIndices[i][j]-1];
00503                 assert( localIndex != mesh->myLoad );
00504                 mesh->h[localIndex] += faultModel->fmUpliftRate[i]->value * faultModel->distFromBoundary[i][j];
00505                 /* Adding the bed-rock uplift*/
00506                 simulationExt->uplift[localIndex] += faultModel->fmUpliftRate[i]->value * faultModel->distFromBoundary[i][j];
00507             }
00508         }
00509     }
00510 }
00511 
00512 Index SPModelFaultModel_Register( PluginsManager* pluginsManager ) {
00513     return PluginsManager_Submit( pluginsManager, SPModelFaultModel_Type, "0", _SPModelFaultModel_DefaultNew );
00514 }
00515