Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17975) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17976) @@ -276,12 +276,4 @@ virtual void Gradj(Vector* gradient,int control_type,int control_index)=0; - virtual IssmDouble ThicknessAbsMisfit(void)=0; - virtual IssmDouble SurfaceLogVxVyMisfit(void)=0; - virtual IssmDouble SurfaceAverageVelMisfit(void)=0; - virtual IssmDouble ThicknessAbsGradient(void)=0; - virtual IssmDouble ThicknessAlongGradient(void)=0; - virtual IssmDouble ThicknessAcrossGradient(void)=0; - virtual IssmDouble BalancethicknessMisfit(void)=0; - virtual IssmDouble RheologyBbarAbsGradient(void)=0; virtual void ControlInputGetGradient(Vector* gradient,int enum_type,int control_index)=0; virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17976) @@ -3329,116 +3329,4 @@ } /*}}}*/ -IssmDouble Penta::SurfaceAverageVelMisfit(void){/*{{{*/ - - int approximation; - IssmDouble J; - Tria* tria=NULL; - - /*retrieve inputs :*/ - inputs->GetInputValue(&approximation,ApproximationEnum); - - /*If on water, return 0: */ - if(!IsIceInElement())return 0; - - /*Bail out if this element if: - * -> Non SSA and not on the surface - * -> SSA (2d model) and not on bed) */ - if ((approximation!=SSAApproximationEnum && !IsOnSurface()) || (approximation==SSAApproximationEnum && !IsOnBase())){ - return 0; - } - else if (approximation==SSAApproximationEnum){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute SurfaceAverageVelMisfit*/ - tria=(Tria*)SpawnTria(0,1,2); - J=tria->SurfaceAverageVelMisfit(); - delete tria->material; delete tria; - return J; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); - J=tria->SurfaceAverageVelMisfit(); - delete tria->material; delete tria; - return J; - } -} -/*}}}*/ -IssmDouble Penta::SurfaceLogVxVyMisfit(void){/*{{{*/ - - IssmDouble J; - Tria* tria=NULL; - - /*inputs: */ - int approximation; - - /*retrieve inputs :*/ - inputs->GetInputValue(&approximation,ApproximationEnum); - - /*If on water, return 0: */ - if(!IsIceInElement())return 0; - - /*Bail out if this element if: - * -> Non SSA and not on the surface - * -> SSA (2d model) and not on bed) */ - if ((approximation!=SSAApproximationEnum && !IsOnSurface()) || (approximation==SSAApproximationEnum && !IsOnBase())){ - return 0; - } - else if (approximation==SSAApproximationEnum){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute SurfaceLogVxVyMisfit*/ - tria=(Tria*)SpawnTria(0,1,2); - J=tria->SurfaceLogVxVyMisfit(); - delete tria->material; delete tria; - return J; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); - J=tria->SurfaceLogVxVyMisfit(); - delete tria->material; delete tria; - return J; - } -} -/*}}}*/ -IssmDouble Penta::ThicknessAbsGradient(void){/*{{{*/ - - _error_("Not implemented yet"); -} -/*}}}*/ -IssmDouble Penta::ThicknessAbsMisfit(void){/*{{{*/ - - int approximation; - IssmDouble J; - Tria* tria=NULL; - - /*retrieve inputs :*/ - inputs->GetInputValue(&approximation,ApproximationEnum); - - /*If on water, return 0: */ - if(!IsIceInElement())return 0; - _error_("Not implemented yet"); - - tria=(Tria*)SpawnTria(0,1,2); - J=tria->ThicknessAbsMisfit(); - delete tria->material; delete tria; - return J; -} -/*}}}*/ -IssmDouble Penta::RheologyBbarAbsGradient(void){/*{{{*/ - - IssmDouble J; - Tria* tria=NULL; - - /*If on water, on shelf or not on bed, skip: */ - if(!IsIceInElement() || !IsOnBase()) return 0; - - tria=(Tria*)SpawnTria(0,1,2); - J=tria->RheologyBbarAbsGradient(); - delete tria->material; delete tria; - return J; -} -/*}}}*/ void Penta::GetVectorFromControlInputs(Vector* vector,int control_enum,int control_index,const char* data){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17975) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17976) @@ -146,12 +146,4 @@ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); - IssmDouble RheologyBbarAbsGradient(void); - IssmDouble ThicknessAbsMisfit(void); - IssmDouble SurfaceLogVxVyMisfit(void); - IssmDouble SurfaceAverageVelMisfit(void); - IssmDouble ThicknessAbsGradient(void); - IssmDouble ThicknessAlongGradient(void){_error_("not supported");}; - IssmDouble ThicknessAcrossGradient(void){_error_("not supported");}; - IssmDouble BalancethicknessMisfit(void){_error_("not supported");}; void InputControlUpdate(IssmDouble scalar,bool save_parameter); IssmDouble Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17975) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17976) @@ -187,12 +187,4 @@ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");}; void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum){_error_("not implemented yet");}; - IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");}; - IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");}; - IssmDouble ThicknessAbsGradient(void){_error_("not implemented yet");}; - IssmDouble ThicknessAlongGradient(void){_error_("not implemented yet");}; - IssmDouble ThicknessAcrossGradient(void){_error_("not implemented yet");}; - IssmDouble BalancethicknessMisfit(void){_error_("not implemented yet");}; - IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");}; - IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");}; void InputControlUpdate(IssmDouble scalar,bool save_parameter){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17975) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17976) @@ -192,12 +192,4 @@ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");}; void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum){_error_("not implemented yet");}; - IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");}; - IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");}; - IssmDouble ThicknessAbsGradient(void){_error_("not implemented yet");}; - IssmDouble ThicknessAlongGradient(void){_error_("not implemented yet");}; - IssmDouble ThicknessAcrossGradient(void){_error_("not implemented yet");}; - IssmDouble BalancethicknessMisfit(void){_error_("not implemented yet");}; - IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");}; - IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");}; void InputControlUpdate(IssmDouble scalar,bool save_parameter){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17976) @@ -2708,61 +2708,4 @@ #endif -IssmDouble Tria::BalancethicknessMisfit(void){/*{{{*/ - - /* Intermediaries */ - IssmDouble Jelem = 0; - IssmDouble weight; - IssmDouble Jdet,temp; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dH[2]; - IssmDouble vx,vy,H; - IssmDouble dvx[2],dvy[2]; - IssmDouble dhdt,basal_melting,surface_mass_balance; - GaussTria *gauss = NULL; - - /*If on water, return 0: */ - if(!IsIceInElement()) return 0; - - /*Retrieve all inputs we will be needing: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,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); - Input* surface_mass_balance_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input); - Input* basal_melting_input = inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); - Input* dhdt_input = inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get all parameters at gaussian point*/ - weights_input->GetInputValue(&weight,gauss,BalancethicknessMisfitEnum); - thickness_input->GetInputValue(&H, gauss); - thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss); - surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss); - basal_melting_input->GetInputValue(&basal_melting,gauss); - dhdt_input->GetInputValue(&dhdt,gauss); - vx_input->GetInputValue(&vx,gauss); - vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); - vy_input->GetInputValue(&vy,gauss); - vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); - - /*Balance thickness soft constraint J = 1/2 (div(Hv)-a)^2*/ - temp = vx*dH[0]+vy*dH[1]+H*(dvx[0]+dvy[1]) - (surface_mass_balance-basal_melting-dhdt); - Jelem+=weight*1/2*temp*temp*Jdet*gauss->weight; - } - - /*Clean up and return*/ - delete gauss; - return Jelem; -} -/*}}}*/ void Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/ @@ -3474,345 +3417,4 @@ } /*}}}*/ -IssmDouble Tria::RheologyBbarAbsGradient(void){/*{{{*/ - - /* Intermediaries */ - IssmDouble Jelem = 0; - IssmDouble weight; - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dp[NDOF2]; - GaussTria *gauss = NULL; - - /*retrieve parameters and inputs*/ - - /*If on water, return 0: */ - if(!IsIceInElement()) return 0; - - /*Retrieve all inputs we will be needing: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* weights_input =inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* rheologyb_input=inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get all parameters at gaussian point*/ - weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum); - rheologyb_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss); - - /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ - Jelem+=weight*1/2*(dp[0]*dp[0] + dp[1]*dp[1])*Jdet*gauss->weight; - } - - /*Clean up and return*/ - delete gauss; - return Jelem; -} -/*}}}*/ -IssmDouble Tria::SurfaceAverageVelMisfit(void){/*{{{*/ - - IssmDouble Jelem=0,S,Jdet; - IssmDouble misfit; - IssmDouble vx,vy,vxobs,vyobs,weight; - IssmDouble xyz_list[NUMVERTICES][3]; - GaussTria *gauss=NULL; - - /*If on water, return 0: */ - if(!IsIceInElement())return 0; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*Retrieve all inputs we will be needing: */ - inputs->GetInputValue(&S,SurfaceAreaEnum); - 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); - Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(3); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get all parameters at gaussian point*/ - weights_input->GetInputValue(&weight,gauss,SurfaceAverageVelMisfitEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - - /*Compute SurfaceAverageVelMisfitEnum: - * - * 1 2 2 - * J = --- sqrt( (u - u ) + (v - v ) ) - * S obs obs - */ - misfit=1/S*sqrt( pow(vx-vxobs,2) + pow(vy-vyobs,2)); - - /*Add to cost function*/ - Jelem+=misfit*weight*Jdet*gauss->weight; - } - - /*clean-up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ -IssmDouble Tria::SurfaceLogVxVyMisfit(void){/*{{{*/ - - IssmDouble Jelem=0, S=0; - IssmDouble epsvel=2.220446049250313e-16; - IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ - IssmDouble misfit, Jdet; - IssmDouble vx,vy,vxobs,vyobs,weight; - IssmDouble xyz_list[NUMVERTICES][3]; - GaussTria *gauss=NULL; - - /*If on water, return 0: */ - if(!IsIceInElement())return 0; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*Retrieve all inputs we will be needing: */ - 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); - Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(4); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get all parameters at gaussian point*/ - weights_input->GetInputValue(&weight,gauss,SurfaceLogVxVyMisfitEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - - /*Compute SurfaceRelVelMisfit: - * - * 1 [ |u| + eps 2 |v| + eps 2 ] - * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | - * 2 [ |u |+ eps |v |+ eps ] - * obs obs - */ - misfit=0.5*pow(meanvel,2)*( - pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2) + - pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),2) ); - - /*Add to cost function*/ - Jelem+=misfit*weight*Jdet*gauss->weight; - } - - /*clean-up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ -IssmDouble Tria::ThicknessAbsGradient(void){/*{{{*/ - - /* Intermediaries */ - IssmDouble Jelem = 0; - IssmDouble weight; - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dp[NDOF2]; - GaussTria *gauss = NULL; - - /*retrieve parameters and inputs*/ - - /*If on water, return 0: */ - if(!IsIceInElement()) return 0; - - /*Retrieve all inputs we will be needing: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* weights_input =inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get all parameters at gaussian point*/ - weights_input->GetInputValue(&weight,gauss,ThicknessAbsGradientEnum); - thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss); - - /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ - Jelem+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight; - } - - /*Clean up and return*/ - delete gauss; - return Jelem; -} -/*}}}*/ -IssmDouble Tria::ThicknessAlongGradient(void){/*{{{*/ - - /* Intermediaries */ - IssmDouble Jelem = 0; - IssmDouble weight; - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dp[NDOF2]; - IssmDouble vx,vy,vel; - GaussTria *gauss = NULL; - - /*retrieve parameters and inputs*/ - - /*If on water, return 0: */ - if(!IsIceInElement()) return 0; - - /*Retrieve all inputs we will be needing: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,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(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get all parameters at gaussian point*/ - weights_input->GetInputValue(&weight,gauss,ThicknessAlongGradientEnum); - thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vel = sqrt(vx*vx+vy*vy); - vx = vx/(vel+1.e-9); - vy = vy/(vel+1.e-9); - - /*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; -} -/*}}}*/ -IssmDouble Tria::ThicknessAcrossGradient(void){/*{{{*/ - - /* Intermediaries */ - IssmDouble Jelem = 0; - IssmDouble weight; - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dp[NDOF2]; - IssmDouble vx,vy,vel; - GaussTria *gauss = NULL; - - /*retrieve parameters and inputs*/ - - /*If on water, return 0: */ - if(!IsIceInElement()) return 0; - - /*Retrieve all inputs we will be needing: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,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(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get all parameters at gaussian point*/ - weights_input->GetInputValue(&weight,gauss,ThicknessAcrossGradientEnum); - thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vel = sqrt(vx*vx+vy*vy); - vx = vx/(vel+1.e-9); - vy = vy/(vel+1.e-9); - - /*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; -} -/*}}}*/ -IssmDouble Tria::ThicknessAbsMisfit(void){/*{{{*/ - - /*Intermediaries*/ - IssmDouble thickness,thicknessobs,weight; - IssmDouble Jdet; - IssmDouble Jelem = 0; - IssmDouble xyz_list[NUMVERTICES][3]; - GaussTria *gauss = NULL; - IssmDouble dH[2]; - - /*If on water, return 0: */ - if(!IsIceInElement())return 0; - - /*Retrieve all inputs we will be needing: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - 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); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get parameters at gauss point*/ - thickness_input->GetInputValue(&thickness,gauss); - thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss); - thicknessobs_input->GetInputValue(&thicknessobs,gauss); - weights_input->GetInputValue(&weight,gauss,ThicknessAbsMisfitEnum); - - /*compute ThicknessAbsMisfit*/ - Jelem+=0.5*(thickness-thicknessobs)*(thickness-thicknessobs)*weight*Jdet*gauss->weight; - } - - /* clean up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ void Tria::GetVectorFromControlInputs(Vector* vector,int control_enum,int control_index,const char* data){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17975) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17976) @@ -156,12 +156,4 @@ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); - IssmDouble RheologyBbarAbsGradient(void); - IssmDouble ThicknessAbsMisfit(void); - IssmDouble ThicknessAbsGradient(void); - IssmDouble ThicknessAlongGradient(void); - IssmDouble ThicknessAcrossGradient(void); - IssmDouble BalancethicknessMisfit(void); - IssmDouble SurfaceLogVxVyMisfit(void); - IssmDouble SurfaceAverageVelMisfit(void); void InputControlUpdate(IssmDouble scalar,bool save_parameter); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 17976) @@ -1134,11 +1134,62 @@ void FemModel::BalancethicknessMisfitx(IssmDouble* presponse){/*{{{*/ - IssmDouble J = 0; + /*output: */ + IssmDouble J=0.; IssmDouble J_sum; - for(int i=0;ielements->Size();i++){ - Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); - J+=element->BalancethicknessMisfit(); - } + IssmDouble weight,vx,vy,H,dvx[2],dvy[2],dH[2]; + IssmDouble temp,Jdet,dhdt,basal_melting,surface_mass_balance; + IssmDouble* xyz_list = NULL; + IssmDouble dp[3]; + + /*Compute Misfit: */ + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + + /*If on water, return 0: */ + if(!element->IsIceInElement()) continue; + + /* Get node coordinates*/ + element->GetVerticesCoordinates(&xyz_list); + Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* surface_mass_balance_input = element->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input); + Input* basal_melting_input = element->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); + Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,BalancethicknessMisfitEnum); + thickness_input->GetInputValue(&H, gauss); + thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss); + surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss); + basal_melting_input->GetInputValue(&basal_melting,gauss); + dhdt_input->GetInputValue(&dhdt,gauss); + vx_input->GetInputValue(&vx,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputValue(&vy,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + + /*Balance thickness soft constraint J = 1/2 (div(Hv)-a)^2*/ + temp = vx*dH[0]+vy*dH[1]+H*(dvx[0]+dvy[1]) - (surface_mass_balance-basal_melting-dhdt); + J +=weight*1/2*temp*temp*Jdet*gauss->weight; + } + + /*clean up and Return: */ + xDelete(xyz_list); + delete gauss; + } + + /*Sum all J from all cpus of the cluster:*/ ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); @@ -1151,16 +1202,47 @@ void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/ - /*Intermediary*/ - int i; - Element* element=NULL; - /*output: */ - IssmDouble J=0; + IssmDouble J=0.; IssmDouble J_sum; + IssmDouble thickness,thicknessobs,weight; + IssmDouble Jdet; + IssmDouble* xyz_list = NULL; + IssmDouble dp[3]; + /*Compute Misfit: */ - for (i=0;iSize();i++){ - element=dynamic_cast(elements->GetObjectByOffset(i)); - J+=element->ThicknessAbsGradient(); + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + + /*If on water, return 0: */ + if(!element->IsIceInElement()) continue; + + /* Get node coordinates*/ + element->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + Input* weights_input =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input =element->GetInput(ThicknessEnum); _assert_(thickness_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,ThicknessAcrossGradientEnum); + thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); + + /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ + J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight; + } + + /*clean up and Return: */ + xDelete(xyz_list); + delete gauss; } Index: /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp (revision 17976) @@ -32,5 +32,4 @@ int domaintype,numcomponents; - bool surfaceintegral; IssmDouble Jelem=0.; IssmDouble misfit,Jdet; @@ -47,16 +46,7 @@ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ - case Domain2DverticalEnum: - surfaceintegral = true; - numcomponents = 1; - break; - case Domain3DEnum: - surfaceintegral = true; - numcomponents = 2; - break; - case Domain2DhorizontalEnum: - surfaceintegral = false; - numcomponents = 2; - break; + case Domain2DverticalEnum: numcomponents = 1; break; + case Domain3DEnum: numcomponents = 2; break; + case Domain2DhorizontalEnum: numcomponents = 2; break; default: _error_("not supported yet"); } Index: /issm/trunk-jpl/src/c/modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.cpp (revision 17976) @@ -10,16 +10,12 @@ void RheologyBbarAbsGradientx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ - /*Intermediary*/ - int i; - Element* element=NULL; - /*output: */ - IssmDouble J=0; + IssmDouble J=0.; IssmDouble J_sum; /*Compute Misfit: */ - for (i=0;iSize();i++){ - element=dynamic_cast(elements->GetObjectByOffset(i)); - J+=element->RheologyBbarAbsGradient(); + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + J+=RheologyBbarAbsGradient(element); } @@ -32,2 +28,59 @@ *pJ=J; } + +IssmDouble RheologyBbarAbsGradient(Element* element){ + + int domaintype,numcomponents; + IssmDouble Jelem=0.; + IssmDouble misfit,Jdet; + IssmDouble dp[2],weight; + IssmDouble* xyz_list = NULL; + + /*Get basal element*/ + if(!element->IsOnBase()) return 0.; + + /*If on water, return 0: */ + if(!element->IsIceInElement()) return 0.; + + /*Get problem dimension*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DverticalEnum: numcomponents = 1; break; + case Domain3DEnum: numcomponents = 2; break; + case Domain2DhorizontalEnum: numcomponents = 2; break; + default: _error_("not supported yet"); + } + + /*Spawn basal element*/ + Element* basalelement = element->SpawnBasalElement(); + + /* Get node coordinates*/ + basalelement->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + Input* weights_input=basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* rheologyb_input=basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum); + rheologyb_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); + + /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ + Jelem+=weight*1/2*(dp[0]*dp[0] + dp[1]*dp[1])*Jdet*gauss->weight; + } + + /*clean up and Return: */ + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + xDelete(xyz_list); + delete gauss; + return Jelem; +} Index: /issm/trunk-jpl/src/c/modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h (revision 17975) +++ /issm/trunk-jpl/src/c/modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h (revision 17976) @@ -10,4 +10,5 @@ /* local prototypes: */ void RheologyBbarAbsGradientx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); +IssmDouble RheologyBbarAbsGradient(Element* element); #endif Index: /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp (revision 17976) @@ -32,5 +32,4 @@ int domaintype,numcomponents; - bool surfaceintegral; IssmDouble Jelem=0.; IssmDouble misfit,Jdet; @@ -47,16 +46,7 @@ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ - case Domain2DverticalEnum: - surfaceintegral = true; - numcomponents = 1; - break; - case Domain3DEnum: - surfaceintegral = true; - numcomponents = 2; - break; - case Domain2DhorizontalEnum: - surfaceintegral = false; - numcomponents = 2; - break; + case Domain2DverticalEnum: numcomponents = 1; break; + case Domain3DEnum: numcomponents = 2; break; + case Domain2DhorizontalEnum: numcomponents = 2; break; default: _error_("not supported yet"); } Index: /issm/trunk-jpl/src/c/modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp (revision 17976) @@ -11,18 +11,12 @@ void SurfaceAverageVelMisfitx(IssmDouble* pJ,FemModel* femmodel){ - /*Intermediary*/ - Element* element=NULL; - /*output: */ - IssmDouble J = 0.; + IssmDouble J=0.; IssmDouble J_sum; - - /*Compute surface area and add to elements inputs */ - SurfaceAreax(NULL,femmodel); /*Compute Misfit: */ for(int i=0;ielements->Size();i++){ - element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); - J+=element->SurfaceAverageVelMisfit(); + Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + J+=SurfaceAverageVelMisfit(element); } @@ -35,2 +29,91 @@ *pJ=J; } + +IssmDouble SurfaceAverageVelMisfit(Element* element){ + + int domaintype,numcomponents; + IssmDouble Jelem=0.; + IssmDouble misfit,S,Jdet; + IssmDouble vx,vy,vxobs,vyobs,weight; + IssmDouble* xyz_list = NULL; + + /*Get basal element*/ + if(!element->IsOnSurface()) return 0.; + + /*If on water, return 0: */ + if(!element->IsIceInElement()) return 0.; + + /*Get problem dimension*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DverticalEnum: + numcomponents = 1; + break; + case Domain3DEnum: + numcomponents = 2; + break; + case Domain2DhorizontalEnum: + numcomponents = 2; + break; + default: _error_("not supported yet"); + } + + /*Spawn surface element*/ + Element* topelement = element->SpawnTopElement(); + + /* Get node coordinates*/ + topelement->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + topelement->GetInputValue(&S,SurfaceAreaEnum); + Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* vx_input =topelement->GetInput(VxEnum); _assert_(vx_input); + Input* vxobs_input =topelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input); + Input* vy_input = NULL; + Input* vyobs_input = NULL; + if(numcomponents==2){ + vy_input =topelement->GetInput(VyEnum); _assert_(vy_input); + vyobs_input =topelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + } + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=topelement->NewGauss(3); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + topelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,SurfaceAverageVelMisfitEnum); + vx_input->GetInputValue(&vx,gauss); + vxobs_input->GetInputValue(&vxobs,gauss); + if(numcomponents==2){ + vy_input->GetInputValue(&vy,gauss); + vyobs_input->GetInputValue(&vyobs,gauss); + } + + /*Compute SurfaceAverageVelMisfitEnum: + * + * 1 2 2 + * J = --- sqrt( (u - u ) + (v - v ) ) + * S obs obs + */ + if(numcomponents==1){ + misfit=1/S*(vx-vxobs)*(vx-vxobs); + } + else{ + misfit=1/S*sqrt( pow(vx-vxobs,2) + pow(vy-vyobs,2)); + } + + /*Add to cost function*/ + Jelem+=misfit*weight*Jdet*gauss->weight; + } + + /*clean up and Return: */ + if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;}; + xDelete(xyz_list); + delete gauss; + return Jelem; +} Index: /issm/trunk-jpl/src/c/modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h (revision 17975) +++ /issm/trunk-jpl/src/c/modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h (revision 17976) @@ -10,4 +10,5 @@ /* local prototypes: */ void SurfaceAverageVelMisfitx(IssmDouble* pJ,FemModel* femmodel); +IssmDouble SurfaceAverageVelMisfit(Element* element); #endif Index: /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp (revision 17976) @@ -32,5 +32,4 @@ int domaintype,numcomponents; - bool surfaceintegral; IssmDouble Jelem=0.; IssmDouble epsvel=2.220446049250313e-16; @@ -50,16 +49,7 @@ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ - case Domain2DverticalEnum: - surfaceintegral = true; - numcomponents = 1; - break; - case Domain3DEnum: - surfaceintegral = true; - numcomponents = 2; - break; - case Domain2DhorizontalEnum: - surfaceintegral = false; - numcomponents = 2; - break; + case Domain2DverticalEnum: numcomponents = 1; break; + case Domain3DEnum: numcomponents = 2; break; + case Domain2DhorizontalEnum: numcomponents = 2; break; default: _error_("not supported yet"); } Index: /issm/trunk-jpl/src/c/modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.cpp (revision 17976) @@ -10,16 +10,12 @@ void SurfaceLogVxVyMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ - /*Intermediary*/ - int i; - Element* element=NULL; - /*output: */ - IssmDouble J=0; + IssmDouble J=0.; IssmDouble J_sum; /*Compute Misfit: */ - for (i=0;iSize();i++){ - element=dynamic_cast(elements->GetObjectByOffset(i)); - J+=element->SurfaceLogVxVyMisfit(); + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + J+=SurfaceLogVxVyMisfit(element); } @@ -32,2 +28,90 @@ *pJ=J; } + +IssmDouble SurfaceLogVxVyMisfit(Element* element){ + + int domaintype,numcomponents; + IssmDouble Jelem=0.; + IssmDouble epsvel=2.220446049250313e-16; + IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ + IssmDouble misfit,Jdet; + IssmDouble vx,vy,vxobs,vyobs,weight; + IssmDouble* xyz_list = NULL; + + /*Get basal element*/ + if(!element->IsOnSurface()) return 0.; + + /*If on water, return 0: */ + if(!element->IsIceInElement()) return 0.; + + /*Get problem dimension*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DverticalEnum: numcomponents = 1; break; + case Domain3DEnum: numcomponents = 2; break; + case Domain2DhorizontalEnum: numcomponents = 2; break; + default: _error_("not supported yet"); + } + + /*Spawn surface element*/ + Element* topelement = element->SpawnTopElement(); + + /* Get node coordinates*/ + topelement->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* vx_input =topelement->GetInput(VxEnum); _assert_(vx_input); + Input* vxobs_input =topelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input); + Input* vy_input = NULL; + Input* vyobs_input = NULL; + if(numcomponents==2){ + vy_input =topelement->GetInput(VyEnum); _assert_(vy_input); + vyobs_input =topelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + } + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=topelement->NewGauss(4); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + topelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,SurfaceLogVxVyMisfitEnum); + vx_input->GetInputValue(&vx,gauss); + vxobs_input->GetInputValue(&vxobs,gauss); + if(numcomponents==2){ + vy_input->GetInputValue(&vy,gauss); + vyobs_input->GetInputValue(&vyobs,gauss); + } + + /*Compute SurfaceRelVelMisfit: + * + * 1 [ |u| + eps 2 |v| + eps 2 ] + * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | + * 2 [ |u |+ eps |v |+ eps ] + * obs obs + */ + + if(numcomponents==1){ + misfit=0.5*meanvel*meanvel*pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2); + } + else{ + misfit=0.5*meanvel*meanvel*( + pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2) + + pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),2) ); + } + + /*Add to cost function*/ + Jelem+=misfit*weight*Jdet*gauss->weight; + } + + /*clean up and Return: */ + if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;}; + xDelete(xyz_list); + delete gauss; + return Jelem; +} Index: /issm/trunk-jpl/src/c/modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h (revision 17975) +++ /issm/trunk-jpl/src/c/modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h (revision 17976) @@ -10,4 +10,5 @@ /* local prototypes: */ void SurfaceLogVxVyMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); +IssmDouble SurfaceLogVxVyMisfit(Element* element); #endif Index: /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp (revision 17976) @@ -32,5 +32,4 @@ int domaintype,numcomponents; - bool surfaceintegral; IssmDouble Jelem=0.; IssmDouble misfit,Jdet,scalex,scaley; @@ -49,16 +48,7 @@ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ - case Domain2DverticalEnum: - surfaceintegral = true; - numcomponents = 1; - break; - case Domain3DEnum: - surfaceintegral = true; - numcomponents = 2; - break; - case Domain2DhorizontalEnum: - surfaceintegral = false; - numcomponents = 2; - break; + case Domain2DverticalEnum: numcomponents = 1; break; + case Domain3DEnum: numcomponents = 2; break; + case Domain2DhorizontalEnum: numcomponents = 2; break; default: _error_("not supported yet"); } Index: /issm/trunk-jpl/src/c/modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.cpp (revision 17976) @@ -10,16 +10,12 @@ void ThicknessAbsMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ - /*Intermediary*/ - int i; - Element* element=NULL; - /*output: */ - IssmDouble J=0; + IssmDouble J=0.; IssmDouble J_sum; /*Compute Misfit: */ - for (i=0;iSize();i++){ - element=dynamic_cast(elements->GetObjectByOffset(i)); - J+=element->ThicknessAbsMisfit(); + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + J+=ThicknessAbsMisfit(element); } @@ -32,2 +28,44 @@ *pJ=J; } + +IssmDouble ThicknessAbsMisfit(Element* element){ + + IssmDouble thickness,thicknessobs,weight; + IssmDouble Jelem=0.; + IssmDouble misfit,Jdet; + IssmDouble* xyz_list = NULL; + + /*If on water, return 0: */ + if(!element->IsIceInElement()) return 0.; + + /* Get node coordinates*/ + element->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + Input* weights_input =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input =element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,ThicknessAbsMisfitEnum); + thickness_input->GetInputValue(&thickness,gauss); + thicknessobs_input->GetInputValue(&thicknessobs,gauss); + + /*Compute ThicknessAbsMisfitEnum*/ + Jelem+=0.5*(thickness-thicknessobs)*(thickness-thicknessobs)*weight*Jdet*gauss->weight; + } + + /*clean up and Return: */ + xDelete(xyz_list); + delete gauss; + return Jelem; +} Index: /issm/trunk-jpl/src/c/modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h (revision 17975) +++ /issm/trunk-jpl/src/c/modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h (revision 17976) @@ -10,4 +10,5 @@ /* local prototypes: */ void ThicknessAbsMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); +IssmDouble ThicknessAbsMisfit(Element* element); #endif Index: /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.cpp (revision 17976) @@ -10,16 +10,12 @@ void ThicknessAcrossGradientx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ - /*Intermediary*/ - int i; - Element* element=NULL; - /*output: */ - IssmDouble J=0; + IssmDouble J=0.; IssmDouble J_sum; /*Compute Misfit: */ - for (i=0;iSize();i++){ - element=dynamic_cast(elements->GetObjectByOffset(i)); - J+=element->ThicknessAcrossGradient(); + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + J+=ThicknessAcrossGradient(element); } @@ -32,2 +28,52 @@ *pJ=J; } + +IssmDouble ThicknessAcrossGradient(Element* element){ + + IssmDouble thickness,thicknessobs,weight; + IssmDouble Jelem=0.; + IssmDouble misfit,Jdet; + IssmDouble* xyz_list = NULL; + IssmDouble dp[3]; + IssmDouble vx,vy,vel; + + /*If on water, return 0: */ + if(!element->IsIceInElement()) return 0.; + + /* Get node coordinates*/ + element->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + Input* weights_input =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input =element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,ThicknessAcrossGradientEnum); + thickness_input->GetInputValue(&thickness,gauss); + thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel = sqrt(vx*vx+vy*vy); + vx = vx/(vel+1.e-9); + vy = vy/(vel+1.e-9); + + /*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: */ + xDelete(xyz_list); + delete gauss; + return Jelem; +} Index: /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h (revision 17975) +++ /issm/trunk-jpl/src/c/modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h (revision 17976) @@ -10,4 +10,5 @@ /* local prototypes: */ void ThicknessAcrossGradientx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); +IssmDouble ThicknessAcrossGradient(Element* element); #endif Index: /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp (revision 17975) +++ /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.cpp (revision 17976) @@ -10,16 +10,12 @@ void ThicknessAlongGradientx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ - /*Intermediary*/ - int i; - Element* element=NULL; - /*output: */ - IssmDouble J=0; + IssmDouble J=0.; IssmDouble J_sum; /*Compute Misfit: */ - for (i=0;iSize();i++){ - element=dynamic_cast(elements->GetObjectByOffset(i)); - J+=element->ThicknessAlongGradient(); + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + J+=ThicknessAlongGradient(element); } @@ -32,2 +28,52 @@ *pJ=J; } + +IssmDouble ThicknessAlongGradient(Element* element){ + + IssmDouble thickness,thicknessobs,weight; + IssmDouble Jelem=0.; + IssmDouble misfit,Jdet; + IssmDouble* xyz_list = NULL; + IssmDouble dp[3]; + IssmDouble vx,vy,vel; + + /*If on water, return 0: */ + if(!element->IsIceInElement()) return 0.; + + /* Get node coordinates*/ + element->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + Input* weights_input =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input =element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,ThicknessAlongGradientEnum); + thickness_input->GetInputValue(&thickness,gauss); + thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel = sqrt(vx*vx+vy*vy); + vx = vx/(vel+1.e-9); + vy = vy/(vel+1.e-9); + + /*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: */ + xDelete(xyz_list); + delete gauss; + return Jelem; +} Index: /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h (revision 17975) +++ /issm/trunk-jpl/src/c/modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h (revision 17976) @@ -10,4 +10,5 @@ /* local prototypes: */ void ThicknessAlongGradientx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); +IssmDouble ThicknessAlongGradient(Element* element); #endif