Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 12926) +++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 12927) @@ -401,6 +401,8 @@ StressTensoryzEnum, StressTensorzzEnum, - IceVolumeEnum, - TotalSmbEnum, + IceVolumeEnum, //FIXME reposition + TotalSmbEnum, //FIXME reposition + ThicknessAlongGradientEnum, + ThicknessAcrossGradientEnum, /*}}}*/ /*Element Interpolations{{{1*/ Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 12926) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 12927) @@ -461,4 +461,8 @@ ./modules/ThicknessAbsGradientx/ThicknessAbsGradientx.cpp\ ./modules/ThicknessAbsGradientx/ThicknessAbsGradientx.h\ + ./modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp\ + ./modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h\ + ./modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp\ + ./modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h\ ./modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.cpp\ ./modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h\ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h (revision 12926) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h (revision 12927) @@ -105,4 +105,6 @@ virtual IssmDouble SurfaceAverageVelMisfit(bool process_units,int weight_index)=0; virtual IssmDouble ThicknessAbsGradient(bool process_units,int weight_index)=0; + virtual IssmDouble ThicknessAlongGradient(bool process_units,int weight_index)=0; + virtual IssmDouble ThicknessAcrossGradient(bool process_units,int weight_index)=0; virtual IssmDouble RheologyBbarAbsGradient(bool process_units,int weight_index)=0; virtual IssmDouble DragCoefficientAbsGradient(bool process_units,int weight_index)=0; Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h (revision 12926) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h (revision 12927) @@ -162,4 +162,6 @@ IssmDouble SurfaceAverageVelMisfit(bool process_units,int weight_index); IssmDouble ThicknessAbsGradient(bool process_units,int weight_index); + IssmDouble ThicknessAlongGradient( bool process_units,int weight_index){_error2_("not supported");}; + IssmDouble ThicknessAcrossGradient(bool process_units,int weight_index){_error2_("not supported");}; void InputControlUpdate(IssmDouble scalar,bool save_parameter); #endif Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 12926) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 12927) @@ -3524,4 +3524,6 @@ case ThicknessAbsMisfitEnum: case ThicknessAbsGradientEnum: + case ThicknessAlongGradientEnum: + case ThicknessAcrossGradientEnum: case SurfaceAbsVelMisfitEnum: case SurfaceRelVelMisfitEnum: @@ -4293,4 +4295,102 @@ } /*}}}*/ +/*FUNCTION Tria::ThicknessAlongGradient{{{*/ +IssmDouble Tria::ThicknessAlongGradient(bool process_units,int weight_index){ + + /* Intermediaries */ + int ig; + IssmDouble Jelem = 0; + IssmDouble weight; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dp[NDOF2]; + IssmDouble vx,vy; + GaussTria *gauss = NULL; + + /*retrieve parameters and inputs*/ + + /*If on water, return 0: */ + if(IsOnWater()) return 0; + + /*Retrieve all inputs we will be needing: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input= inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,weight_index); + thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + + /*J = 1/2 ( vx*dH/dx + vy*dH/dy )^2 */ + Jelem+=weight*1/2*(vx*dp[0] + vy*dp[1])*(vx*dp[0] + vy*dp[1])*Jdet*gauss->weight; + } + + /*Clean up and return*/ + delete gauss; + return Jelem; +} +/*}}}*/ +/*FUNCTION Tria::ThicknessAcrossGradient{{{*/ +IssmDouble Tria::ThicknessAcrossGradient(bool process_units,int weight_index){ + + /* Intermediaries */ + int ig; + IssmDouble Jelem = 0; + IssmDouble weight; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dp[NDOF2]; + IssmDouble vx,vy; + GaussTria *gauss = NULL; + + /*retrieve parameters and inputs*/ + + /*If on water, return 0: */ + if(IsOnWater()) return 0; + + /*Retrieve all inputs we will be needing: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input= inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,weight_index); + thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + + /*J = 1/2 ( -vy*dH/dx + vx*dH/dy )^2 */ + Jelem+=weight*1/2*(-vy*dp[0] + vx*dp[1])*(-vy*dp[0] + vx*dp[1])*Jdet*gauss->weight; + } + + /*Clean up and return*/ + delete gauss; + return Jelem; +} +/*}}}*/ /*FUNCTION Tria::ThicknessAbsMisfit {{{*/ IssmDouble Tria::ThicknessAbsMisfit(bool process_units,int weight_index){ @@ -4346,13 +4446,14 @@ /*Intermediaries */ int i,ig,resp; - IssmDouble Jdet; - IssmDouble thickness,thicknessobs,weight; - int *responses = NULL; + IssmDouble Jdet; + IssmDouble thickness,thicknessobs,weight; int num_responses; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[3]; - IssmDouble dbasis[NDOF2][NUMVERTICES]; - IssmDouble dH[2]; - GaussTria* gauss=NULL; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[3]; + IssmDouble dbasis[NDOF2][NUMVERTICES]; + IssmDouble dH[2]; + IssmDouble v[2]; + GaussTria *gauss = NULL; + int *responses = NULL; /*Initialize Element vector*/ @@ -4363,7 +4464,9 @@ this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); - Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* thicknessobs_input = inputs->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input); - Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* thicknessobs_input = inputs->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); + Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); /* Start looping on the number of gaussian points: */ @@ -4392,4 +4495,16 @@ for(i=0;ivalues[i]+= - weight*dH[0]*dbasis[0][i]*Jdet*gauss->weight; for(i=0;ivalues[i]+= - weight*dH[1]*dbasis[1][i]*Jdet*gauss->weight; + break; + case ThicknessAlongGradientEnum: + weights_input->GetInputValue(&weight, gauss,resp); + vx_input->GetInputValue(&v[0],gauss); + vy_input->GetInputValue(&v[1],gauss); + for(i=0;ivalues[i]+= - weight*(dH[0]*v[0]+dH[1]*v[1])*(dbasis[0][i]*v[0]+dbasis[1][i]*v[1])*Jdet*gauss->weight; + break; + case ThicknessAcrossGradientEnum: + weights_input->GetInputValue(&weight, gauss,resp); + vx_input->GetInputValue(&v[0],gauss); + vy_input->GetInputValue(&v[1],gauss); + for(i=0;ivalues[i]+= - weight*(dH[0]*(-v[1])+dH[1]*v[0])*(dbasis[0][i]*(-v[1])+dbasis[1][i]*v[0])*Jdet*gauss->weight; break; default: @@ -4568,4 +4683,10 @@ /*Nothing in P vector*/ break; + case ThicknessAlongGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAcrossGradientEnum: + /*Nothing in P vector*/ + break; case RheologyBbarAbsGradientEnum: /*Nothing in P vector*/ @@ -4745,4 +4866,10 @@ /*Nothing in P vector*/ break; + case ThicknessAcrossGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAlongGradientEnum: + /*Nothing in P vector*/ + break; case RheologyBbarAbsGradientEnum: /*Nothing in P vector*/ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 12926) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 12927) @@ -158,5 +158,7 @@ IssmDouble ThicknessAbsMisfit( bool process_units,int weight_index); IssmDouble SurfaceAbsVelMisfit( bool process_units,int weight_index); - IssmDouble ThicknessAbsGradient(bool process_units,int weight_index); + IssmDouble ThicknessAbsGradient( bool process_units,int weight_index); + IssmDouble ThicknessAlongGradient( bool process_units,int weight_index); + IssmDouble ThicknessAcrossGradient(bool process_units,int weight_index); IssmDouble SurfaceRelVelMisfit( bool process_units,int weight_index); IssmDouble SurfaceLogVelMisfit( bool process_units,int weight_index); Index: /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12926) +++ /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12927) @@ -394,4 +394,6 @@ case IceVolumeEnum : return "IceVolume"; case TotalSmbEnum : return "TotalSmb"; + case ThicknessAlongGradientEnum : return "ThicknessAlongGradient"; + case ThicknessAcrossGradientEnum : return "ThicknessAcrossGradient"; case P0Enum : return "P0"; case P1Enum : return "P1"; Index: /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12926) +++ /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12927) @@ -404,4 +404,6 @@ else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; + else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; + else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; else if (strcmp(name,"P0")==0) return P0Enum; else if (strcmp(name,"P1")==0) return P1Enum; Index: /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp (revision 12927) +++ /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp (revision 12927) @@ -0,0 +1,37 @@ +/*!\file ThicknessAcrossGradientx + * \brief: compute misfit between observations and model + */ + +#include "./ThicknessAcrossGradientx.h" + +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void ThicknessAcrossGradientx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units, int weight_index){ + + /*Intermediary*/ + int i; + Element* element=NULL; + + /*output: */ + double J=0; + double J_sum; + + /*Compute Misfit: */ + for (i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + J+=element->ThicknessAcrossGradient(process_units,weight_index); + } + + /*Sum all J from all cpus of the cluster:*/ + #ifdef _HAVE_MPI_ + MPI_Reduce (&J,&J_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&J_sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + J=J_sum; + #endif + + /*Assign output pointers: */ + *pJ=J; +} Index: /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h (revision 12927) +++ /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h (revision 12927) @@ -0,0 +1,14 @@ +/*!\file: ThicknessAcrossGradientx.h + * \brief header file for inverse methods misfit computation + */ + +#ifndef _THICKNESSACROSSGRADIENT_H +#define _THICKNESSACROSSGRADIENT_H + +#include "../../Container/Container.h" +#include "../../classes/objects/objects.h" + +/* local prototypes: */ +void ThicknessAcrossGradientx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index); + +#endif Index: /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp (revision 12927) +++ /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp (revision 12927) @@ -0,0 +1,37 @@ +/*!\file ThicknessAlongGradientx + * \brief: compute misfit between observations and model + */ + +#include "./ThicknessAlongGradientx.h" + +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void ThicknessAlongGradientx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units, int weight_index){ + + /*Intermediary*/ + int i; + Element* element=NULL; + + /*output: */ + double J=0; + double J_sum; + + /*Compute Misfit: */ + for (i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + J+=element->ThicknessAlongGradient(process_units,weight_index); + } + + /*Sum all J from all cpus of the cluster:*/ + #ifdef _HAVE_MPI_ + MPI_Reduce (&J,&J_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&J_sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + J=J_sum; + #endif + + /*Assign output pointers: */ + *pJ=J; +} Index: /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h (revision 12927) +++ /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h (revision 12927) @@ -0,0 +1,14 @@ +/*!\file: ThicknessAlongGradientx.h + * \brief header file for inverse methods misfit computation + */ + +#ifndef _THICKNESSALONGGRADIENT_H +#define _THICKNESSALONGGRADIENT_H + +#include "../../Container/Container.h" +#include "../../classes/objects/objects.h" + +/* local prototypes: */ +void ThicknessAlongGradientx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index); + +#endif Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 12926) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 12927) @@ -119,4 +119,6 @@ #include "./ThicknessAbsMisfitx/ThicknessAbsMisfitx.h" #include "./ThicknessAbsGradientx/ThicknessAbsGradientx.h" +#include "./ThicknessAlongGradientx/ThicknessAlongGradientx.h" +#include "./ThicknessAcrossGradientx/ThicknessAcrossGradientx.h" #include "./UpdateVertexPositionsx/UpdateVertexPositionsx.h" #include "./UpdateConstraintsx/UpdateConstraintsx.h"