Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17974) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17975) @@ -277,7 +277,4 @@ virtual void Gradj(Vector* gradient,int control_type,int control_index)=0; virtual IssmDouble ThicknessAbsMisfit(void)=0; - virtual IssmDouble SurfaceAbsVelMisfit(void)=0; - virtual IssmDouble SurfaceRelVelMisfit(void)=0; - virtual IssmDouble SurfaceLogVelMisfit(void)=0; virtual IssmDouble SurfaceLogVxVyMisfit(void)=0; virtual IssmDouble SurfaceAverageVelMisfit(void)=0; @@ -287,5 +284,4 @@ virtual IssmDouble BalancethicknessMisfit(void)=0; virtual IssmDouble RheologyBbarAbsGradient(void)=0; - virtual IssmDouble DragCoefficientAbsGradient(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 17974) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17975) @@ -3365,76 +3365,4 @@ } /*}}}*/ -IssmDouble Penta::SurfaceAbsVelMisfit(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 SurfaceAbsVelMisfit*/ - tria=(Tria*)SpawnTria(0,1,2); - J=tria->SurfaceAbsVelMisfit(); - delete tria->material; delete tria; - return J; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); - J=tria->SurfaceAbsVelMisfit(); - delete tria->material; delete tria; - return J; - } -} -/*}}}*/ -IssmDouble Penta::SurfaceLogVelMisfit(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 SurfaceLogVelMisfit*/ - tria=(Tria*)SpawnTria(0,1,2); //lower face is 0, upper face is 1. - J=tria->SurfaceLogVelMisfit(); - delete tria->material; delete tria; - return J; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); //lower face is 0, upper face is 1. - J=tria->SurfaceLogVelMisfit(); - delete tria->material; delete tria; - return J; - } -} -/*}}}*/ IssmDouble Penta::SurfaceLogVxVyMisfit(void){/*{{{*/ @@ -3475,5 +3403,10 @@ } /*}}}*/ -IssmDouble Penta::SurfaceRelVelMisfit(void){/*{{{*/ +IssmDouble Penta::ThicknessAbsGradient(void){/*{{{*/ + + _error_("Not implemented yet"); +} +/*}}}*/ +IssmDouble Penta::ThicknessAbsMisfit(void){/*{{{*/ int approximation; @@ -3486,63 +3419,8 @@ /*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 SurfaceRelVelMisfit*/ - tria=(Tria*)SpawnTria(0,1,2); - J=tria->SurfaceRelVelMisfit(); - delete tria->material; delete tria; - return J; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); - J=tria->SurfaceRelVelMisfit(); - 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::DragCoefficientAbsGradient(void){/*{{{*/ - - IssmDouble J; - Tria* tria=NULL; - - /*If on water, on shelf or not on bed, skip: */ - if(!IsIceInElement()|| IsFloating() || !IsOnBase()) return 0; - - tria=(Tria*)SpawnTria(0,1,2); //lower face is 0, upper face is 1 - J=tria->DragCoefficientAbsGradient(); delete tria->material; delete tria; return J; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17974) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17975) @@ -132,5 +132,4 @@ #endif - IssmDouble DragCoefficientAbsGradient(void); void GradientIndexing(int* indexing,int control_index); void Gradj(Vector* gradient,int control_type,int control_index); @@ -149,7 +148,4 @@ IssmDouble RheologyBbarAbsGradient(void); IssmDouble ThicknessAbsMisfit(void); - IssmDouble SurfaceAbsVelMisfit(void); - IssmDouble SurfaceRelVelMisfit(void); - IssmDouble SurfaceLogVelMisfit(void); IssmDouble SurfaceLogVxVyMisfit(void); IssmDouble SurfaceAverageVelMisfit(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17974) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17975) @@ -168,5 +168,4 @@ #endif - IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");}; void GradientIndexing(int* indexing,int control_index); void Gradj(Vector* gradient,int control_type,int control_index){_error_("not implemented yet");}; @@ -190,11 +189,8 @@ IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");}; IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");}; - IssmDouble SurfaceAbsVelMisfit(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 SurfaceRelVelMisfit(void){_error_("not implemented yet");}; - IssmDouble SurfaceLogVelMisfit(void){_error_("not implemented yet");}; IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");}; IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17974) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17975) @@ -194,11 +194,8 @@ IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");}; IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");}; - IssmDouble SurfaceAbsVelMisfit(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 SurfaceRelVelMisfit(void){_error_("not implemented yet");}; - IssmDouble SurfaceLogVelMisfit(void){_error_("not implemented yet");}; IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");}; IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17974) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17975) @@ -3571,63 +3571,4 @@ } /*}}}*/ -IssmDouble Tria::SurfaceLogVelMisfit(void){/*{{{*/ - - IssmDouble Jelem=0.; - IssmDouble misfit,Jdet; - IssmDouble epsvel=2.220446049250313e-16; - IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ - IssmDouble velocity_mag,obs_velocity_mag; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble vx,vy,vxobs,vyobs,weight; - 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,SurfaceLogVelMisfitEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - - /*Compute SurfaceLogVelMisfit: - * 4 [ vel + eps ] 2 - * J = 4 \bar{v}^2 | log ( ----------- ) | - * [ vel + eps ] - * obs - */ - velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; - obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; - misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2); - - /*Add to cost function*/ - Jelem+=misfit*weight*Jdet*gauss->weight; - } - - /*clean-up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ IssmDouble Tria::SurfaceLogVxVyMisfit(void){/*{{{*/ @@ -3689,119 +3630,4 @@ } /*}}}*/ -IssmDouble Tria::SurfaceAbsVelMisfit(void){/*{{{*/ - - IssmDouble Jelem=0; - 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(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,SurfaceAbsVelMisfitEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - - /*Compute SurfaceAbsVelMisfitEnum: - * - * 1 [ 2 2 ] - * J = --- | (u - u ) + (v - v ) | - * 2 [ obs obs ] - * - */ - misfit=0.5*( 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::SurfaceRelVelMisfit(void){/*{{{*/ - - IssmDouble Jelem=0; - IssmDouble scalex=1,scaley=1; - IssmDouble misfit,Jdet; - IssmDouble epsvel=2.220446049250313e-16; - IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ - 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,SurfaceRelVelMisfitEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - - /*Compute SurfaceRelVelMisfit: - * - * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] - * J = --- | ------------- (u - u ) + ------------- (v - v ) | - * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] - * obs obs - */ - scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; - scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; - misfit=0.5*(scalex*pow((vx-vxobs),2)+scaley*pow((vy-vyobs),2)); - - /*Add to cost function*/ - Jelem+=misfit*weight*Jdet*gauss->weight; - } - - /*clean up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ IssmDouble Tria::ThicknessAbsGradient(void){/*{{{*/ @@ -3985,46 +3811,4 @@ /* clean up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ -IssmDouble Tria::DragCoefficientAbsGradient(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* drag_input =inputs->GetInput(FrictionCoefficientEnum); _assert_(drag_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,DragCoefficientAbsGradientEnum); - drag_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; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17974) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17975) @@ -137,5 +137,4 @@ #endif - IssmDouble DragCoefficientAbsGradient(void); void GradientIndexing(int* indexing,int control_index); void Gradj(Vector* gradient,int control_type,int control_index); @@ -159,11 +158,8 @@ IssmDouble RheologyBbarAbsGradient(void); IssmDouble ThicknessAbsMisfit(void); - IssmDouble SurfaceAbsVelMisfit(void); IssmDouble ThicknessAbsGradient(void); IssmDouble ThicknessAlongGradient(void); IssmDouble ThicknessAcrossGradient(void); IssmDouble BalancethicknessMisfit(void); - IssmDouble SurfaceRelVelMisfit(void); - IssmDouble SurfaceLogVelMisfit(void); IssmDouble SurfaceLogVxVyMisfit(void); IssmDouble SurfaceAverageVelMisfit(void); Index: /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp (revision 17974) +++ /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp (revision 17975) @@ -10,16 +10,12 @@ void SurfaceLogVelMisfitx( 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->SurfaceLogVelMisfit(); + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + J+=SurfaceLogVelMisfit(element); } @@ -32,2 +28,101 @@ *pJ=J; } + +IssmDouble SurfaceLogVelMisfit(Element* element){ + + int domaintype,numcomponents; + bool surfaceintegral; + IssmDouble Jelem=0.; + IssmDouble epsvel=2.220446049250313e-16; + IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ + IssmDouble velocity_mag,obs_velocity_mag; + 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: + surfaceintegral = true; + numcomponents = 1; + break; + case Domain3DEnum: + surfaceintegral = true; + numcomponents = 2; + break; + case Domain2DhorizontalEnum: + surfaceintegral = false; + 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,SurfaceLogVelMisfitEnum); + vx_input->GetInputValue(&vx,gauss); + vxobs_input->GetInputValue(&vxobs,gauss); + if(numcomponents==2){ + vy_input->GetInputValue(&vy,gauss); + vyobs_input->GetInputValue(&vyobs,gauss); + } + + /*Compute SurfaceLogVelMisfit: + * [ vel + eps ] 2 + * J = 4 \bar{v}^2 | log ( ----------- ) | + * [ vel + eps ] + * obs + */ + if(numcomponents==1){ + velocity_mag =fabs(vx)+epsvel; + obs_velocity_mag=fabs(vxobs)+epsvel; + } + else{ + velocity_mag =sqrt(vx*vx+vy*vy)+epsvel; + obs_velocity_mag=sqrt(vxobs*vxobs+vyobs*vyobs)+epsvel; + } + + misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),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/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h (revision 17974) +++ /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h (revision 17975) @@ -10,4 +10,5 @@ /* local prototypes: */ void SurfaceLogVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); +IssmDouble SurfaceLogVelMisfit(Element* element); #endif Index: /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp (revision 17974) +++ /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp (revision 17975) @@ -10,16 +10,12 @@ void SurfaceRelVelMisfitx( 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->SurfaceRelVelMisfit(); + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + J+=SurfaceRelVelMisfit(element); } @@ -32,2 +28,101 @@ *pJ=J; } + +IssmDouble SurfaceRelVelMisfit(Element* element){ + + int domaintype,numcomponents; + bool surfaceintegral; + IssmDouble Jelem=0.; + IssmDouble misfit,Jdet,scalex,scaley; + IssmDouble epsvel=2.220446049250313e-16; + IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ + 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: + surfaceintegral = true; + numcomponents = 1; + break; + case Domain3DEnum: + surfaceintegral = true; + numcomponents = 2; + break; + case Domain2DhorizontalEnum: + surfaceintegral = false; + 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,SurfaceRelVelMisfitEnum); + 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 [ \bar{v}^2 2 \bar{v}^2 2 ] + * J = --- | ------------- (u - u ) + ------------- (v - v ) | + * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] + * obs obs + */ + + if(numcomponents==2){ + scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; + scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; + misfit=0.5*(scalex*pow((vx-vxobs),2)+scaley*pow((vy-vyobs),2)); + } + else{ + scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; + misfit=0.5*(scalex*pow((vx-vxobs),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/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h (revision 17974) +++ /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h (revision 17975) @@ -10,4 +10,5 @@ /* local prototypes: */ void SurfaceRelVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); +IssmDouble SurfaceRelVelMisfit(Element* element); #endif