Changeset 12927


Ignore:
Timestamp:
08/07/12 10:45:59 (13 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added 2 new cost functions for balance thickness (ThicknessAcrossGradient and ThicknessAlongGradient)

Location:
issm/trunk-jpl/src/c
Files:
6 added
9 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h

    r12773 r12927  
    401401        StressTensoryzEnum,
    402402        StressTensorzzEnum,
    403         IceVolumeEnum,
    404         TotalSmbEnum,
     403        IceVolumeEnum, //FIXME reposition
     404        TotalSmbEnum,  //FIXME reposition
     405        ThicknessAlongGradientEnum,
     406        ThicknessAcrossGradientEnum,
    405407        /*}}}*/
    406408        /*Element Interpolations{{{1*/

  • issm/trunk-jpl/src/c/Makefile.am

    r12922 r12927  
    461461                                          ./modules/ThicknessAbsGradientx/ThicknessAbsGradientx.cpp\
    462462                                          ./modules/ThicknessAbsGradientx/ThicknessAbsGradientx.h\
     463                                          ./modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp\
     464                                          ./modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h\
     465                                          ./modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp\
     466                                          ./modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h\
    463467                                          ./modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.cpp\
    464468                                          ./modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h\

  • issm/trunk-jpl/src/c/classes/objects/Elements/Element.h

    r12832 r12927  
    105105                virtual IssmDouble SurfaceAverageVelMisfit(bool process_units,int weight_index)=0;
    106106                virtual IssmDouble ThicknessAbsGradient(bool process_units,int weight_index)=0;
     107                virtual IssmDouble ThicknessAlongGradient(bool process_units,int weight_index)=0;
     108                virtual IssmDouble ThicknessAcrossGradient(bool process_units,int weight_index)=0;
    107109                virtual IssmDouble RheologyBbarAbsGradient(bool process_units,int weight_index)=0;
    108110                virtual IssmDouble DragCoefficientAbsGradient(bool process_units,int weight_index)=0;

  • issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h

    r12832 r12927  
    162162                IssmDouble SurfaceAverageVelMisfit(bool process_units,int weight_index);
    163163                IssmDouble ThicknessAbsGradient(bool process_units,int weight_index);
     164                IssmDouble ThicknessAlongGradient( bool process_units,int weight_index){_error2_("not supported");};
     165                IssmDouble ThicknessAcrossGradient(bool process_units,int weight_index){_error2_("not supported");};
    164166                void   InputControlUpdate(IssmDouble scalar,bool save_parameter);
    165167                #endif

  • issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp

    r12893 r12927  
    35243524                case ThicknessAbsMisfitEnum:
    35253525                case ThicknessAbsGradientEnum:
     3526                case ThicknessAlongGradientEnum:
     3527                case ThicknessAcrossGradientEnum:
    35263528                case SurfaceAbsVelMisfitEnum:
    35273529                case SurfaceRelVelMisfitEnum:
    … …  
    42934295}
    42944296/*}}}*/
     4297/*FUNCTION Tria::ThicknessAlongGradient{{{*/
     4298IssmDouble Tria::ThicknessAlongGradient(bool process_units,int weight_index){
     4299
     4300        /* Intermediaries */
     4301        int         ig;
     4302        IssmDouble  Jelem = 0;
     4303        IssmDouble  weight;
     4304        IssmDouble  Jdet;
     4305        IssmDouble  xyz_list[NUMVERTICES][3];
     4306        IssmDouble  dp[NDOF2];
     4307        IssmDouble  vx,vy;
     4308        GaussTria  *gauss                    = NULL;
     4309
     4310        /*retrieve parameters and inputs*/
     4311
     4312        /*If on water, return 0: */
     4313        if(IsOnWater()) return 0;
     4314
     4315        /*Retrieve all inputs we will be needing: */
     4316        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4317        Input* weights_input  = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     4318        Input* thickness_input= inputs->GetInput(ThicknessEnum);                          _assert_(thickness_input);
     4319        Input* vx_input       = inputs->GetInput(VxEnum);                                 _assert_(vx_input);
     4320        Input* vy_input       = inputs->GetInput(VyEnum);                                 _assert_(vy_input);
     4321
     4322        /* Start looping on the number of gaussian points: */
     4323        gauss=new GaussTria(2);
     4324        for (ig=gauss->begin();ig<gauss->end();ig++){
     4325
     4326                gauss->GaussPoint(ig);
     4327
     4328                /* Get Jacobian determinant: */
     4329                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     4330
     4331                /*Get all parameters at gaussian point*/
     4332                weights_input->GetInputValue(&weight,gauss,weight_index);
     4333                thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
     4334                vx_input->GetInputValue(&vx,gauss);
     4335                vy_input->GetInputValue(&vy,gauss);
     4336
     4337                /*J = 1/2 ( vx*dH/dx + vy*dH/dy )^2 */
     4338                Jelem+=weight*1/2*(vx*dp[0] + vy*dp[1])*(vx*dp[0] + vy*dp[1])*Jdet*gauss->weight;
     4339        }
     4340
     4341        /*Clean up and return*/
     4342        delete gauss;
     4343        return Jelem;
     4344}
     4345/*}}}*/
     4346/*FUNCTION Tria::ThicknessAcrossGradient{{{*/
     4347IssmDouble Tria::ThicknessAcrossGradient(bool process_units,int weight_index){
     4348
     4349        /* Intermediaries */
     4350        int         ig;
     4351        IssmDouble  Jelem = 0;
     4352        IssmDouble  weight;
     4353        IssmDouble  Jdet;
     4354        IssmDouble  xyz_list[NUMVERTICES][3];
     4355        IssmDouble  dp[NDOF2];
     4356        IssmDouble  vx,vy;
     4357        GaussTria  *gauss                    = NULL;
     4358
     4359        /*retrieve parameters and inputs*/
     4360
     4361        /*If on water, return 0: */
     4362        if(IsOnWater()) return 0;
     4363
     4364        /*Retrieve all inputs we will be needing: */
     4365        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4366        Input* weights_input  = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     4367        Input* thickness_input= inputs->GetInput(ThicknessEnum);                          _assert_(thickness_input);
     4368        Input* vx_input       = inputs->GetInput(VxEnum);                                 _assert_(vx_input);
     4369        Input* vy_input       = inputs->GetInput(VyEnum);                                 _assert_(vy_input);
     4370
     4371        /* Start looping on the number of gaussian points: */
     4372        gauss=new GaussTria(2);
     4373        for (ig=gauss->begin();ig<gauss->end();ig++){
     4374
     4375                gauss->GaussPoint(ig);
     4376
     4377                /* Get Jacobian determinant: */
     4378                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     4379
     4380                /*Get all parameters at gaussian point*/
     4381                weights_input->GetInputValue(&weight,gauss,weight_index);
     4382                thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
     4383                vx_input->GetInputValue(&vx,gauss);
     4384                vy_input->GetInputValue(&vy,gauss);
     4385
     4386                /*J = 1/2 ( -vy*dH/dx + vx*dH/dy )^2 */
     4387                Jelem+=weight*1/2*(-vy*dp[0] + vx*dp[1])*(-vy*dp[0] + vx*dp[1])*Jdet*gauss->weight;
     4388        }
     4389
     4390        /*Clean up and return*/
     4391        delete gauss;
     4392        return Jelem;
     4393}
     4394/*}}}*/
    42954395/*FUNCTION Tria::ThicknessAbsMisfit {{{*/
    42964396IssmDouble Tria::ThicknessAbsMisfit(bool process_units,int weight_index){
    … …  
    43464446        /*Intermediaries */
    43474447        int         i,ig,resp;
    4348         IssmDouble      Jdet;
    4349         IssmDouble      thickness,thicknessobs,weight;
    4350         int        *responses = NULL;
     4448        IssmDouble  Jdet;
     4449        IssmDouble  thickness,thicknessobs,weight;
    43514450        int         num_responses;
    4352         IssmDouble      xyz_list[NUMVERTICES][3];
    4353         IssmDouble      basis[3];
    4354         IssmDouble      dbasis[NDOF2][NUMVERTICES];
    4355         IssmDouble      dH[2];
    4356         GaussTria*  gauss=NULL;
     4451        IssmDouble  xyz_list[NUMVERTICES][3];
     4452        IssmDouble  basis[3];
     4453        IssmDouble  dbasis[NDOF2][NUMVERTICES];
     4454        IssmDouble  dH[2];
     4455        IssmDouble  v[2];
     4456        GaussTria *gauss     = NULL;
     4457        int       *responses = NULL;
    43574458
    43584459        /*Initialize Element vector*/
    … …  
    43634464        this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    43644465        this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
    4365         Input* thickness_input    = inputs->GetInput(ThicknessEnum);   _assert_(thickness_input);
    4366         Input* thicknessobs_input = inputs->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
    4367         Input* weights_input      = inputs->GetInput(InversionCostFunctionsCoefficientsEnum);     _assert_(weights_input);
     4466        Input* thickness_input    = inputs->GetInput(ThicknessEnum);                          _assert_(thickness_input);
     4467        Input* thicknessobs_input = inputs->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
     4468        Input* weights_input      = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     4469        Input* vx_input           = inputs->GetInput(VxEnum);                                 _assert_(vx_input);
     4470        Input* vy_input           = inputs->GetInput(VyEnum);                                 _assert_(vy_input);
    43684471
    43694472        /* Start  looping on the number of gaussian points: */
    … …  
    43924495                                for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[0]*dbasis[0][i]*Jdet*gauss->weight;
    43934496                                for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[1]*dbasis[1][i]*Jdet*gauss->weight;
     4497                                break;
     4498                        case ThicknessAlongGradientEnum:
     4499                                weights_input->GetInputValue(&weight, gauss,resp);
     4500                                vx_input->GetInputValue(&v[0],gauss);
     4501                                vy_input->GetInputValue(&v[1],gauss);
     4502                                for(i=0;i<numdof;i++) pe->values[i]+= - weight*(dH[0]*v[0]+dH[1]*v[1])*(dbasis[0][i]*v[0]+dbasis[1][i]*v[1])*Jdet*gauss->weight;
     4503                                break;
     4504                        case ThicknessAcrossGradientEnum:
     4505                                weights_input->GetInputValue(&weight, gauss,resp);
     4506                                vx_input->GetInputValue(&v[0],gauss);
     4507                                vy_input->GetInputValue(&v[1],gauss);
     4508                                for(i=0;i<numdof;i++) pe->values[i]+= - weight*(dH[0]*(-v[1])+dH[1]*v[0])*(dbasis[0][i]*(-v[1])+dbasis[1][i]*v[0])*Jdet*gauss->weight;
    43944509                                break;
    43954510                        default:
    … …  
    45684683                                        /*Nothing in P vector*/
    45694684                                        break;
     4685                                case ThicknessAlongGradientEnum:
     4686                                        /*Nothing in P vector*/
     4687                                        break;
     4688                                case ThicknessAcrossGradientEnum:
     4689                                        /*Nothing in P vector*/
     4690                                        break;
    45704691                                case RheologyBbarAbsGradientEnum:
    45714692                                        /*Nothing in P vector*/
    … …  
    47454866                                        /*Nothing in P vector*/
    47464867                                        break;
     4868                                case ThicknessAcrossGradientEnum:
     4869                                        /*Nothing in P vector*/
     4870                                        break;
     4871                                case ThicknessAlongGradientEnum:
     4872                                        /*Nothing in P vector*/
     4873                                        break;
    47474874                                case RheologyBbarAbsGradientEnum:
    47484875                                        /*Nothing in P vector*/

  • issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h

    r12832 r12927  
    158158                IssmDouble ThicknessAbsMisfit(     bool process_units,int weight_index);
    159159                IssmDouble SurfaceAbsVelMisfit(    bool process_units,int weight_index);
    160                 IssmDouble ThicknessAbsGradient(bool process_units,int weight_index);
     160                IssmDouble ThicknessAbsGradient(   bool process_units,int weight_index);
     161                IssmDouble ThicknessAlongGradient( bool process_units,int weight_index);
     162                IssmDouble ThicknessAcrossGradient(bool process_units,int weight_index);
    161163                IssmDouble SurfaceRelVelMisfit(    bool process_units,int weight_index);
    162164                IssmDouble SurfaceLogVelMisfit(    bool process_units,int weight_index);

  • issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r12773 r12927  
    394394                case IceVolumeEnum : return "IceVolume";
    395395                case TotalSmbEnum : return "TotalSmb";
     396                case ThicknessAlongGradientEnum : return "ThicknessAlongGradient";
     397                case ThicknessAcrossGradientEnum : return "ThicknessAcrossGradient";
    396398                case P0Enum : return "P0";
    397399                case P1Enum : return "P1";

  • issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r12773 r12927  
    404404              else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
    405405              else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
     406              else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
     407              else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum;
    406408              else if (strcmp(name,"P0")==0) return P0Enum;
    407409              else if (strcmp(name,"P1")==0) return P1Enum;

  • issm/trunk-jpl/src/c/modules/modules.h

    r12773 r12927  
    119119#include "./ThicknessAbsMisfitx/ThicknessAbsMisfitx.h"
    120120#include "./ThicknessAbsGradientx/ThicknessAbsGradientx.h"
     121#include "./ThicknessAlongGradientx/ThicknessAlongGradientx.h"
     122#include "./ThicknessAcrossGradientx/ThicknessAcrossGradientx.h"
    121123#include "./UpdateVertexPositionsx/UpdateVertexPositionsx.h"
    122124#include "./UpdateConstraintsx/UpdateConstraintsx.h"
Note: See TracChangeset for help on using the changeset viewer.