Index: /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 18061) +++ /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 18062) @@ -480,5 +480,118 @@ }/*}}}*/ void BalancethicknessAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ - _error_("Not implemented yet"); + /* WARNING: this gradient is valid for Soft balance thickness only */ + + /*If on water, grad = 0: */ + if(!element->IsIceInElement()) return; + + /*Intermediaries*/ + IssmDouble Jdet,weight; + IssmDouble thickness,thicknessobs,dH[3],dp[3]; + IssmDouble vx,vy,vel,dvx[2],dvy[2],dhdt,basal_melting,surface_mass_balance; + IssmDouble *xyz_list= NULL; + + /*Get list of cost functions*/ + int *responses = NULL; + int num_responses,resp,solution; + element->FindParam(&num_responses,InversionNumCostFunctionsEnum); + element->FindParam(&responses,NULL,InversionCostFunctionsEnum); + element->FindParam(&solution,SolutionTypeEnum); + if(solution!=BalancethicknessSoftSolutionEnum) _error_("not implemented yet"); + if(control_type!=ThicknessEnum) _error_("Control "<GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); + IssmDouble* dbasis = xNew(2*numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Retrieve all inputs we will be needing: */ + element->GetVerticesCoordinates(&xyz_list); + element->GradientIndexing(&vertexpidlist[0],control_index); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* thicknessobs_input = element->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); + Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_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(4); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + thickness_input->GetInputValue(&thickness, gauss); + thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss); + thicknessobs_input->GetInputValue(&thicknessobs, gauss); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsP1(basis,gauss); + element->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss); + + /*Deal with first part (partial derivative a J with respect to k)*/ + for(resp=0;respGetInputValue(&weight,gauss,responses[resp]); + + switch(responses[resp]){ + case ThicknessAbsMisfitEnum: + for(int i=0;iweight*basis[i]; + break; + case ThicknessAbsGradientEnum: + for(int i=0;iweight; + for(int i=0;iweight; + break; + case ThicknessAlongGradientEnum: + 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); + for(int i=0;iweight; + break; + case ThicknessAcrossGradientEnum: + 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); + for(int i=0;iweight; + break; + case BalancethicknessMisfitEnum: + 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); + for(int i=0;iweight*( + (vx*dH[0]+vy*dH[1] + thickness*(dvx[0]+dvy[1]))*(vx*dbasis[0*numvertices+i]+ vy*dbasis[1*numvertices+i] + basis[i]*(dvx[0]+dvy[1])) + -(surface_mass_balance-basal_melting-dhdt)*(vx*dbasis[0*numvertices+i]+ vy*dbasis[1*numvertices+i] + basis[i]*(dvx[0]+dvy[1])) + ); + } + break; + default: + _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); + } + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + xDelete(responses); + delete gauss; + + }/*}}}*/ void BalancethicknessAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18061) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18062) @@ -277,5 +277,4 @@ #endif - virtual void Gradj(Vector* gradient,int control_type,int control_index)=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 18061) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18062) @@ -2895,324 +2895,4 @@ }/*}}}*/ -void Penta::Gradj(Vector* gradient,int control_type,int control_index){/*{{{*/ - /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/ - - int approximation; - Tria* tria=NULL; - - /*If on water, skip grad (=0): */ - if(!IsIceInElement())return; - - /*First deal with ∂/∂alpha(KU-F)*/ - switch(control_type){ - - case FrictionCoefficientEnum: - inputs->GetInputValue(&approximation,ApproximationEnum); - switch(approximation){ - case SSAApproximationEnum: - GradjDragSSA(gradient,control_index); - break; - case HOApproximationEnum: - GradjDragHO(gradient,control_index); - break; - case FSApproximationEnum: - GradjDragFS(gradient,control_index); - break; - case NoneApproximationEnum: - /*Gradient is 0*/ - break; - default: - _error_("approximation " << EnumToStringx(approximation) << " not supported yet"); - } - break; - - case MaterialsRheologyBbarEnum: - inputs->GetInputValue(&approximation,ApproximationEnum); - switch(approximation){ - case SSAApproximationEnum: - GradjBbarSSA(gradient,control_index); - break; - case HOApproximationEnum: - GradjBbarHO(gradient,control_index); - break; - case FSApproximationEnum: - GradjBbarFS(gradient,control_index); - break; - case NoneApproximationEnum: - /*Gradient is 0*/ - break; - default: - _error_("approximation " << EnumToStringx(approximation) << " not supported yet"); - } - break; - - default: - _error_("control type " << EnumToStringx(control_type) << " not supported yet: "); - } - - /*Now deal with ∂J/∂alpha*/ - int *responses = NULL; - int num_responses,resp; - this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); - this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum); - - for(resp=0;respGradjDragGradient(gradient,control_index); - delete tria->material; delete tria; - } - break; - case RheologyBbarAbsGradientEnum: - if(IsOnBase()){ - tria=(Tria*)SpawnTria(0,1,2); - tria->GradjBGradient(gradient,control_index); - delete tria->material; delete tria; - } - break; - default: - _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); - } - xDelete(responses); -} -/*}}}*/ -void Penta::GradjDragSSA(Vector* gradient,int control_index){/*{{{*/ - - /*Gradient is 0 if on shelf or not on bed*/ - if(IsFloating() || !IsOnBase()) return; - - /*Spawn tria*/ - Tria* tria=(Tria*)SpawnTria(0,1,2); - tria->GradjDragSSA(gradient,control_index); - delete tria->material; delete tria; - -} /*}}}*/ -void Penta::GradjDragHO(Vector* gradient,int control_index){/*{{{*/ - - printf("-------------- file: Penta.cpp line: %i\n",__LINE__); - int i,j; - int analysis_type; - int vertexpidlist[NUMVERTICES]; - IssmDouble vx,vy,lambda,mu,alpha_complement,Jdet; - IssmDouble bed,thickness,Neff,drag; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; - IssmDouble dk[NDOF3]; - IssmDouble grade_g[NUMVERTICES]={0.0}; - IssmDouble grade_g_gaussian[NUMVERTICES]; - IssmDouble basis[6]; - Friction* friction=NULL; - GaussPenta *gauss=NULL; - - /*Gradient is 0 if on shelf or not on bed*/ - if(IsFloating() || !IsOnBase()) return; - - /*Retrieve all inputs and parameters*/ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - GradientIndexing(&vertexpidlist[0],control_index); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;iGetInput(AdjointxEnum); _assert_(adjointx_input); - Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input); - - /*Build frictoin element, needed later: */ - friction=new Friction(this,2); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(0,1,2,4); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0],gauss); - GetNodalFunctionsP1(basis,gauss); - - /*Build alpha_complement_list: */ - friction->GetAlphaComplement(&alpha_complement,gauss); - - dragcoefficient_input->GetInputValue(&drag, gauss); - adjointx_input->GetInputValue(&lambda, gauss); - adjointy_input->GetInputValue(&mu, gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;iweight*basis[i]; - } - - /*Add gradje_g_gaussian vector to gradje_g: */ - for(i=0;i(grade_g[i])); - grade_g[i]+=grade_g_gaussian[i]; - } - } - gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); - - /*Clean up and return*/ - delete gauss; - delete friction; -} -/*}}}*/ -void Penta::GradjDragFS(Vector* gradient,int control_index){/*{{{*/ - - int i,j; - int analysis_type; - int vertexpidlist[NUMVERTICES]; - IssmDouble bed,thickness,Neff; - IssmDouble lambda,mu,xi,Jdet,vx,vy,vz; - IssmDouble alpha_complement,drag; - IssmDouble surface_normal[3],bed_normal[3]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; - IssmDouble dk[NDOF3]; - IssmDouble basis[6]; - IssmDouble grade_g[NUMVERTICES]={0.0}; - IssmDouble grade_g_gaussian[NUMVERTICES]; - Friction* friction=NULL; - GaussPenta* gauss=NULL; - - /*Gradient is 0 if on shelf or not on bed*/ - if(IsFloating() || !IsOnBase()) return; - - /*Retrieve all inputs and parameters*/ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;iGetInput(FrictionCoefficientEnum); _assert_(drag_input); - Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input =inputs->GetInput(VzEnum); _assert_(vz_input); - Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); - Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); - Input* adjointz_input=inputs->GetInput(AdjointzEnum); _assert_(adjointz_input); - - /*Build frictoin element, needed later: */ - friction=new Friction(this,3); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(0,1,2,4); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /*Recover alpha_complement and drag: */ - friction->GetAlphaComplement(&alpha_complement,gauss); - drag_input->GetInputValue(&drag,gauss); - - /*recover lambda mu and xi: */ - adjointx_input->GetInputValue(&lambda,gauss); - adjointy_input->GetInputValue(&mu ,gauss); - adjointz_input->GetInputValue(&xi ,gauss); - - /*recover vx vy and vz: */ - vx_input->GetInputValue(&vx, gauss); - vy_input->GetInputValue(&vy, gauss); - vz_input->GetInputValue(&vz, gauss); - - /*Get normal vector to the bed */ - NormalTop(&surface_normal[0],&xyz_list_tria[0][0]); - - bed_normal[0]=-surface_normal[0]; //Function is for upper surface, so the normal to the bed is the opposite of the result - bed_normal[1]=-surface_normal[1]; - bed_normal[2]=-surface_normal[2]; - - /* Get Jacobian determinant: */ - GetTriaJacobianDeterminant(&Jdet,&xyz_list_tria[0][0],gauss); - GetNodalFunctionsP1(basis, gauss); - - /*Get k derivative: dk/dx */ - drag_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;iweight*basis[i]; - } - - /*Add gradje_g_gaussian vector to gradje_g: */ - for( i=0; iSetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); - - delete friction; - delete gauss; -} -/*}}}*/ -void Penta::GradjBbarSSA(Vector* gradient,int control_index){/*{{{*/ - - /*This element should be collapsed into a tria element at its base*/ - if (!IsOnBase()) return; - - /*Depth Average B*/ - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); - if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum); - - /*Collapse element to the base*/ - Tria* tria=(Tria*)SpawnTria(0,1,2); - tria->GradjBSSA(gradient,control_index); - delete tria->material; delete tria; - - /*delete Average B*/ - this->inputs->DeleteInput(MaterialsRheologyBbarEnum); - this->inputs->DeleteInput(DamageDbarEnum); - -} /*}}}*/ -void Penta::GradjBbarHO(Vector* gradient,int control_index){/*{{{*/ - - /*Gradient is computed on bed only (Bbar)*/ - if (!IsOnBase()) return; - - /*Depth Average B and D*/ - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); - if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum); - - /*Collapse element to the base*/ - Tria* tria=(Tria*)SpawnTria(0,1,2); - tria->GradjBSSA(gradient,control_index); //We use SSA as an estimate for now - delete tria->material; delete tria; - - /*delete Average B*/ - this->inputs->DeleteInput(MaterialsRheologyBbarEnum); - this->inputs->DeleteInput(DamageDbarEnum); -} /*}}}*/ -void Penta::GradjBbarFS(Vector* gradient,int control_index){/*{{{*/ - - /*Gradient is computed on bed only (Bbar)*/ - if (!IsOnBase()) return; - - /*Depth Average B and D*/ - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); - if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum); - - /*Collapse element to the base*/ - Tria* tria=(Tria*)SpawnTria(0,1,2); - tria->GradjBSSA(gradient,control_index); //We use SSA as an estimate for now - delete tria->material; delete tria; - - /*delete Average B*/ - this->inputs->DeleteInput(MaterialsRheologyBbarEnum); - this->inputs->DeleteInput(DamageDbarEnum); -} /*}}}*/ void Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18061) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18062) @@ -131,11 +131,4 @@ #endif - void Gradj(Vector* gradient,int control_type,int control_index); - void GradjDragSSA(Vector* gradient,int control_index); - void GradjDragHO(Vector* gradient,int control_index); - void GradjDragFS(Vector* gradient,int control_index); - void GradjBbarSSA(Vector* gradient,int control_index); - void GradjBbarHO(Vector* gradient,int control_index); - void GradjBbarFS(Vector* gradient,int control_index); void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data); void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 18061) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 18062) @@ -116,113 +116,4 @@ }/*}}}*/ -void Seg::GradjDragFS(Vector* gradient,int control_index){/*{{{*/ - - int i; - int analysis_type; - int vertexpidlist[NUMVERTICES]; - int connectivity[NUMVERTICES]; - IssmDouble vx,lambda,alpha_complement,drag,Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dk[NDOF2]; - IssmDouble grade_g[NUMVERTICES]={0.0}; - IssmDouble grade_g_gaussian[NUMVERTICES]; - IssmDouble basis[3]; - Friction* friction=NULL; - GaussSeg *gauss=NULL; - - if(IsFloating())return; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GradientIndexing(&vertexpidlist[0],control_index); - this->GetVerticesConnectivityList(&connectivity[0]); - - /*Build frictoin element, needed later: */ - friction=new Friction(this,1); - - /*Retrieve all inputs we will be needing: */ - Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussSeg(4); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis, gauss); - - /*Build alpha_complement_list: */ - friction->GetAlphaComplement(&alpha_complement,gauss); - - dragcoefficient_input->GetInputValue(&drag, gauss); - adjointx_input->GetInputValue(&lambda, gauss); - vx_input->GetInputValue(&vx,gauss); - dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;iweight*basis[i]; - } - - /*Add gradje_g_gaussian vector to gradje_g: */ - for(i=0;i(grade_g[i])); - grade_g[i]+=grade_g_gaussian[i]; - } - } - gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); - - /*Clean up and return*/ - delete gauss; - delete friction; -} -/*}}}*/ -void Seg::GradjDragGradient(Vector* gradient,int control_index){/*{{{*/ - - int i; - int vertexpidlist[NUMVERTICES]; - IssmDouble Jdet,weight; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dbasis[NDOF2][NUMVERTICES]; - IssmDouble dk[NDOF2]; - IssmDouble grade_g[NUMVERTICES]={0.0}; - GaussSeg *gauss=NULL; - - /*Retrieve all inputs we will be needing: */ - if(IsFloating())return; - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GradientIndexing(&vertexpidlist[0],control_index); - Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input); - Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussSeg(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); - weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum); - - /*Build alpha_complement_list: */ - dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;iweight*dbasis[0][i]*dk[0]; - _assert_(!xIsNan(grade_g[i])); - } - } - gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ bool Seg::IsIcefront(void){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 18061) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 18062) @@ -167,16 +167,4 @@ #endif - void Gradj(Vector* gradient,int control_type,int control_index){_error_("not implemented yet");}; - void GradjBGradient(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjDGradient(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjBSSA(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjDSSA(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjDragSSA(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjDragFS(Vector* gradient,int control_index); - void GradjDragGradient(Vector* gradient,int control_index); - void GradjDhDtBalancedthickness(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjVxBalancedthickness(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjVyBalancedthickness(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjThicknessBalancethicknessSoft(Vector* gradient,int control_index){_error_("not implemented yet");}; void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data){_error_("not implemented yet");}; void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18061) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18062) @@ -173,16 +173,4 @@ IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");}; void GradientIndexing(int* indexing,int control_index){_error_("not implemented yet");}; - void Gradj(Vector* gradient,int control_type,int control_index){_error_("not implemented yet");}; - void GradjBGradient(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjDGradient(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjBSSA(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjDSSA(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjDragSSA(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjDragFS(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjDragGradient(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjDhDtBalancedthickness(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjVxBalancedthickness(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjVyBalancedthickness(Vector* gradient,int control_index){_error_("not implemented yet");}; - void GradjThicknessBalancethicknessSoft(Vector* gradient,int control_index){_error_("not implemented yet");}; void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data){_error_("not implemented yet");}; void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18061) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18062) @@ -2773,648 +2773,4 @@ }/*}}}*/ -void Tria::Gradj(Vector* gradient,int control_type,int control_index){/*{{{*/ - /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/ - - int approximation; - Seg* seg=NULL; - - /*If on water, grad = 0: */ - if(!IsIceInElement()) return; - - /*First deal with ∂/∂alpha(KU-F)*/ - switch(control_type){ - case FrictionCoefficientEnum: - inputs->GetInputValue(&approximation,ApproximationEnum); - switch(approximation){ - case SSAApproximationEnum: - GradjDragSSA(gradient,control_index); - break; - case FSApproximationEnum: - GradjDragFS(gradient,control_index); - break; - case NoneApproximationEnum: - /*Gradient is 0*/ - break; - default: - _error_("approximation " << EnumToStringx(approximation) << " not supported yet"); - } - break; - case MaterialsRheologyBbarEnum: - GradjBSSA(gradient,control_index); - break; - case DamageDbarEnum: - GradjDSSA(gradient,control_index); - break; - case BalancethicknessThickeningRateEnum: - GradjDhDtBalancedthickness(gradient,control_index); - break; - case VxEnum: - GradjVxBalancedthickness(gradient,control_index); - break; - case VyEnum: - GradjVyBalancedthickness(gradient,control_index); - break; - case ThicknessEnum: - GradjThicknessBalancethicknessSoft(gradient,control_index); - break; - case BalancethicknessApparentMassbalanceEnum: - GradjAdotBulancethickness2(gradient,control_index); - break; - default: - _error_("control type not supported yet: " << EnumToStringx(control_type)); - } - - /*Now deal with ∂J/∂alpha*/ - int *responses = NULL; - int num_responses,resp; - this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); - this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum); - - for(resp=0;respGetInputValue(&approximation,ApproximationEnum); - switch(approximation){ - case SSAApproximationEnum: - GradjDragGradient(gradient,control_index); - break; - case FSApproximationEnum:{ - if(IsFloating() || !IsOnBase()) return; - int index1,index2; - this->EdgeOnBaseIndices(&index1,&index2); - Seg* seg = SpawnSeg(index1,index2); - seg->GradjDragGradient(gradient,control_index); - seg->DeleteMaterials(); delete seg; - break; - } - case NoneApproximationEnum: - /*Gradient is 0*/ - break; - default: - _error_("approximation " << EnumToStringx(approximation) << " not supported yet"); - } - break; - case RheologyBbarAbsGradientEnum: - GradjBGradient(gradient,control_index); - break; - default: - _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); - } - - xDelete(responses); -} -/*}}}*/ -void Tria::GradjBGradient(Vector* gradient,int control_index){/*{{{*/ - - int i; - int vertexpidlist[NUMVERTICES]; - IssmDouble Jdet,weight; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dbasis[NDOF2][NUMVERTICES]; - IssmDouble dk[NDOF2]; - IssmDouble grade_g[NUMVERTICES]={0.0}; - GaussTria *gauss=NULL; - - /*Retrieve all inputs we will be needing: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GradientIndexing(&vertexpidlist[0],control_index); - Input* rheologyb_input=inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_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); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); - weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum); - - /*Build alpha_complement_list: */ - rheologyb_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;iweight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]); - } - gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ -void Tria::GradjBSSA(Vector* gradient,int control_index){/*{{{*/ - - /*Intermediaries*/ - int i; - int doflist[NUMVERTICES]; - IssmDouble vx,vy,lambda,mu,thickness,Jdet; - IssmDouble viscosity_complement; - IssmDouble dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[3],epsilon[3]; - IssmDouble grad[NUMVERTICES]={0.0}; - GaussTria *gauss = NULL; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GradientIndexing(&doflist[0],control_index); - - /*Retrieve all inputs*/ - 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* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); - Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); - Input* rheologyb_input=inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(4); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - thickness_input->GetInputValue(&thickness,gauss); - rheologyb_input->GetInputDerivativeValue(&dB[0],&xyz_list[0][0],gauss); - vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); - vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); - adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss); - adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss); - - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - material->GetViscosityComplement(&viscosity_complement,&epsilon[0]); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - /*standard gradient dJ/dki*/ - for (i=0;iweight*basis[i]; - } - - gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL); - - /*clean-up*/ - delete gauss; -} -/*}}}*/ -void Tria::GradjDSSA(Vector* gradient,int control_index){/*{{{*/ - - /*Intermediaries*/ - int i; - int doflist[NUMVERTICES]; - IssmDouble vx,vy,lambda,mu,thickness,Jdet; - IssmDouble viscosity_complement; - IssmDouble dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[3],epsilon[3]; - IssmDouble grad[NUMVERTICES]={0.0}; - GaussTria *gauss = NULL; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GradientIndexing(&doflist[0],control_index); - - /*Retrieve all inputs*/ - 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* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); - Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); - if(this->material->IsDamage()){ - Input* rheologyd_input=inputs->GetInput(DamageDbarEnum); _assert_(rheologyd_input); - } - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(4); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - thickness_input->GetInputValue(&thickness,gauss); - vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); - vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); - adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss); - adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss); - - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - material->GetViscosityDComplement(&viscosity_complement,&epsilon[0]); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - /*standard gradient dJ/dki*/ - for (i=0;iweight*basis[i]; - } - - gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL); - - /*clean-up*/ - delete gauss; -} -/*}}}*/ -void Tria::GradjDragSSA(Vector* gradient,int control_index){/*{{{*/ - - int i; - int analysis_type; - int vertexpidlist[NUMVERTICES]; - int connectivity[NUMVERTICES]; - IssmDouble vx,vy,lambda,mu,alpha_complement,Jdet; - IssmDouble bed,thickness,Neff,drag; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dk[NDOF2]; - IssmDouble grade_g[NUMVERTICES]={0.0}; - IssmDouble grade_g_gaussian[NUMVERTICES]; - IssmDouble basis[3]; - Friction* friction=NULL; - GaussTria *gauss=NULL; - - if(IsFloating())return; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GradientIndexing(&vertexpidlist[0],control_index); - this->GetVerticesConnectivityList(&connectivity[0]); - - /*Build frictoin element, needed later: */ - friction=new Friction(this,2); - - /*Retrieve all inputs we will be needing: */ - Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); - Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(4); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis, gauss); - - /*Build alpha_complement_list: */ - friction->GetAlphaComplement(&alpha_complement,gauss); - - dragcoefficient_input->GetInputValue(&drag, gauss); - adjointx_input->GetInputValue(&lambda, gauss); - adjointy_input->GetInputValue(&mu, gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;iweight*basis[i]; - } - - /*Add gradje_g_gaussian vector to gradje_g: */ - for(i=0;i(grade_g[i])); - grade_g[i]+=grade_g_gaussian[i]; - } - } - /*Analytical gradient*/ - //delete gauss; - //gauss=new GaussTria(); - //for (int iv=0;ivGaussVertex(iv); - // friction->GetAlphaComplement(&alpha_complement,gauss); - // dragcoefficient_input->GetInputValue(&drag, gauss); - // adjointx_input->GetInputValue(&lambda, gauss); - // adjointy_input->GetInputValue(&mu, gauss); - // vx_input->GetInputValue(&vx,gauss); - // vy_input->GetInputValue(&vy,gauss); - // grade_g[iv] = -2*1.e+7*drag*alpha_complement*(lambda*vx+mu*vy)/((IssmDouble)connectivity[iv]); - //} - /*End Analytical gradient*/ - - gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); - - /*Clean up and return*/ - delete gauss; - delete friction; -} -/*}}}*/ -void Tria::GradjDragFS(Vector* gradient,int control_index){/*{{{*/ - - /*Gradient is 0 if on shelf or not on bed*/ - if(IsFloating() || !IsOnBase()) return; - - int index1,index2; - this->EdgeOnBaseIndices(&index1,&index2); - Seg* seg = SpawnSeg(index1,index2); - seg->GradjDragFS(gradient,control_index); - seg->DeleteMaterials(); delete seg; -} -/*}}}*/ -void Tria::GradjDragGradient(Vector* gradient,int control_index){/*{{{*/ - - int i; - int vertexpidlist[NUMVERTICES]; - IssmDouble Jdet,weight; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dbasis[NDOF2][NUMVERTICES]; - IssmDouble dk[NDOF2]; - IssmDouble grade_g[NUMVERTICES]={0.0}; - GaussTria *gauss=NULL; - - /*Retrieve all inputs we will be needing: */ - if(IsFloating())return; - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GradientIndexing(&vertexpidlist[0],control_index); - Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_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); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); - weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum); - - /*Build alpha_complement_list: */ - dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;iweight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]); - _assert_(!xIsNan(grade_g[i])); - } - } - gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ -void Tria::GradjDhDtBalancedthickness(Vector* gradient,int control_index){/*{{{*/ - - /*Intermediaries*/ - int vertexpidlist[NUMVERTICES]; - IssmDouble lambda[NUMVERTICES]; - IssmDouble gradient_g[NUMVERTICES]; - - /*Compute Gradient*/ - GradientIndexing(&vertexpidlist[0],control_index); - GetInputListOnVertices(&lambda[0],AdjointEnum); - for(int i=0;iSetValues(NUMVERTICES,vertexpidlist,gradient_g,INS_VAL); -} -/*}}}*/ -void Tria::GradjAdotBulancethickness2(Vector* gradient,int control_index){/*{{{*/ - - /*Intermediaries*/ - IssmDouble lambda,potential,weight,Jdet; - IssmDouble grade_g[NUMVERTICES]; - IssmDouble basis[NUMVERTICES]; - IssmDouble xyz_list[NUMVERTICES][3]; - int vertexpidlist[NUMVERTICES]; - - /*Recover inputs*/ - Input* lambda_input = inputs->GetInput(AdjointEnum); _assert_(lambda_input); - Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GradientIndexing(&vertexpidlist[0],control_index); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(&basis[0],gauss); - weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum); - lambda_input->GetInputValue(&lambda,gauss); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for(int i=0;iweight*basis[i]*lambda; - } - } - gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ -void Tria::GradjVxBalancedthickness(Vector* gradient,int control_index){/*{{{*/ - - /*Intermediaries*/ - int i; - int vertexpidlist[NUMVERTICES]; - IssmDouble thickness,Jdet; - IssmDouble basis[3]; - IssmDouble Dlambda[2],dp[2]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble grade_g[NUMVERTICES] = {0.0}; - GaussTria *gauss = NULL; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GradientIndexing(&vertexpidlist[0],control_index); - - /*Retrieve all inputs we will be needing: */ - Input* adjoint_input=inputs->GetInput(AdjointEnum); _assert_(adjoint_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); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis, gauss); - - adjoint_input->GetInputDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); - thickness_input->GetInputValue(&thickness, gauss); - thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss); - - for(i=0;iweight*basis[i]; - } - - gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ -void Tria::GradjVyBalancedthickness(Vector* gradient,int control_index){/*{{{*/ - - /*Intermediaries*/ - int i; - int vertexpidlist[NUMVERTICES]; - IssmDouble thickness,Jdet; - IssmDouble basis[3]; - IssmDouble Dlambda[2],dp[2]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble grade_g[NUMVERTICES] = {0.0}; - GaussTria *gauss = NULL; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GradientIndexing(&vertexpidlist[0],control_index); - - /*Retrieve all inputs we will be needing: */ - Input* adjoint_input=inputs->GetInput(AdjointEnum); _assert_(adjoint_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); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis, gauss); - - adjoint_input->GetInputDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); - thickness_input->GetInputValue(&thickness, gauss); - thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss); - - for(i=0;iweight*basis[i]; - } - gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ -void Tria::GradjThicknessBalancethicknessSoft(Vector* gradient,int control_index){/*{{{*/ - - /*Intermediaries */ - int i,resp; - int vertexpidlist[NUMVERTICES]; - IssmDouble Jdet; - IssmDouble thickness,thicknessobs,weight; - int num_responses; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[3]; - IssmDouble dbasis[NDOF2][NUMVERTICES]; - IssmDouble dH[2]; - IssmDouble vx,vy,vel; - IssmDouble dvx[2],dvy[2]; - IssmDouble dhdt,basal_melting,surface_mass_balance; - GaussTria *gauss = NULL; - int *responses = NULL; - IssmDouble grade_g[NUMVERTICES] = {0.0}; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GradientIndexing(&vertexpidlist[0],control_index); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); - this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum); - 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); - 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); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis, gauss); - GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); - - thickness_input->GetInputValue(&thickness, gauss); - thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss); - thicknessobs_input->GetInputValue(&thicknessobs, gauss); - - /*Loop over all requested responses*/ - for(resp=0;respGetInputValue(&weight,gauss,responses[resp]); - - switch(responses[resp]){ - - case ThicknessAbsMisfitEnum: - for(i=0;iweight*basis[i]; - break; - case ThicknessAbsGradientEnum: - for(i=0;iweight; - for(i=0;iweight; - break; - case ThicknessAlongGradientEnum: - 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); - for(i=0;iweight; - break; - case ThicknessAcrossGradientEnum: - 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); - for(i=0;iweight; - break; - case BalancethicknessMisfitEnum: - 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); - for(i=0;iweight*( - (vx*dH[0]+vy*dH[1] + thickness*(dvx[0]+dvy[1]))*(vx*dbasis[0][i]+ vy*dbasis[1][i] + basis[i]*(dvx[0]+dvy[1])) - -(surface_mass_balance-basal_melting-dhdt)*(vx*dbasis[0][i]+ vy*dbasis[1][i] + basis[i]*(dvx[0]+dvy[1])) - ); - } - break; - default: - _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); - } - } - } - - gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); - - /*Clean up and return*/ - delete gauss; - xDelete(responses); -} -/*}}}*/ 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 18061) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18062) @@ -136,17 +136,4 @@ #endif - void Gradj(Vector* gradient,int control_type,int control_index); - void GradjBGradient(Vector* gradient,int control_index); - void GradjDGradient(Vector* gradient,int control_index); - void GradjBSSA(Vector* gradient,int control_index); - void GradjDSSA(Vector* gradient,int control_index); - void GradjDragSSA(Vector* gradient,int control_index); - void GradjDragFS(Vector* gradient,int control_index); - void GradjDragGradient(Vector* gradient,int control_index); - void GradjDhDtBalancedthickness(Vector* gradient,int control_index); - void GradjVxBalancedthickness(Vector* gradient,int control_index); - void GradjVyBalancedthickness(Vector* gradient,int control_index); - void GradjThicknessBalancethicknessSoft(Vector* gradient,int control_index); - void GradjAdotBulancethickness2(Vector* gradient,int control_index); void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data); void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);