Index: ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18059) +++ ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18060) @@ -922,7 +922,9 @@ element->FindParam(&responses,NULL,InversionCostFunctionsEnum); /*Check that control_type is supported*/ - if(control_type!=MaterialsRheologyBbarEnum || control_type!=FrictionCoefficientEnum || control_type!=DamageDbarEnum){ + if(control_type!=MaterialsRheologyBbarEnum && + control_type!=FrictionCoefficientEnum && + control_type!=DamageDbarEnum){ _error_("Control "<* gradient,int control_index){/*{{{*/ - /*return if floating*/ + /*return if floating (gradient is 0)*/ if(element->IsFloating()) return; /*Intermediaries*/ @@ -1011,7 +1013,7 @@ /*Initialize some vectors*/ IssmDouble* dbasis = xNew(2*numvertices); - IssmDouble* ge = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); int* vertexpidlist = xNew(numvertices); /*Retrieve all inputs we will be needing: */ @@ -1055,16 +1057,309 @@ }/*}}}*/ void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); + + /*Intermediaries*/ + int domaintype,dim; + Element* basalelement; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + dim = 2; + break; + case Domain2DverticalEnum: + if(!element->IsOnBase()) return; + basalelement = element->SpawnBasalElement(); + dim = 1; + break; + case Domain3DEnum: + if(!element->IsOnBase()) return; + basalelement = element->SpawnBasalElement(); + dim = 2; + break; + default: _error_("mesh "<GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* dbasis = xNew(2*numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Retrieve all inputs we will be needing: */ + basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->GradientIndexing(&vertexpidlist[0],control_index); + Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); + Input* weights_input = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss); + weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum); + + /*Build alpha_complement_list: */ + rheologyb_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for(int i=0;iweight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]); + } + else{ + ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0]; + } + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(dbasis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); + + /*return if floating (gradient is 0)*/ + if(element->IsFloating()) return; + + /*Intermediaries*/ + int domaintype,dim; + Element* basalelement; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + dim = 2; + break; + case Domain2DverticalEnum: + if(!element->IsOnBase()) return; + basalelement = element->SpawnBasalElement(); + dim = 1; + break; + case Domain3DEnum: + if(!element->IsOnBase()) return; + basalelement = element->SpawnBasalElement(); + dim = 2; + break; + default: _error_("mesh "<GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Build friction element, needed later: */ + Friction* friction=new Friction(basalelement,dim); + + /*Retrieve all inputs we will be needing: */ + basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->GradientIndexing(&vertexpidlist[0],control_index); + Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); + Input* adjointx_input = basalelement->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input = basalelement->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* dragcoeff_input = basalelement->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(4); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + adjointx_input->GetInputValue(&lambda, gauss); + adjointy_input->GetInputValue(&mu, gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + dragcoeff_input->GetInputValue(&drag, gauss); + + friction->GetAlphaComplement(&dalpha2dk,gauss); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctionsP1(basis,gauss); + + /*Build gradient vector (actually -dJ/dD): */ + for(int i=0;iweight*basis[i]; + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; + delete friction; + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ void AdjointHorizAnalysis::GradientJDragHO(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); + + /*return if floating or not on bed (gradient is 0)*/ + if(element->IsFloating()) return; + if(!element->IsOnBase()) return; + + /*Intermediaries*/ + int dim=3; + IssmDouble Jdet,weight; + IssmDouble drag,dalpha2dk; + IssmDouble vx,vy,lambda,mu; + IssmDouble *xyz_list_base= NULL; + + /*Fetch number of vertices for this finite element*/ + int numvertices = element->GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Build friction element, needed later: */ + Friction* friction=new Friction(element,dim); + + /*Retrieve all inputs we will be needing: */ + element->GetVerticesCoordinatesBase(&xyz_list_base); + element->GradientIndexing(&vertexpidlist[0],control_index); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGaussBase(4); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + adjointx_input->GetInputValue(&lambda, gauss); + adjointy_input->GetInputValue(&mu, gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + dragcoeff_input->GetInputValue(&drag, gauss); + + friction->GetAlphaComplement(&dalpha2dk,gauss); + + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctionsP1(basis,gauss); + + /*Build gradient vector (actually -dJ/dD): */ + for(int i=0;iweight*basis[i]; + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list_base); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; + delete friction; }/*}}}*/ void AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); + + /*return if floating or not on bed (gradient is 0)*/ + if(element->IsFloating()) return; + if(!element->IsOnBase()) return; + + /*Intermediaries*/ + int dim=3; + IssmDouble Jdet,weight; + IssmDouble drag,dalpha2dk,normal[3]; + IssmDouble vx,vy,vz,lambda,mu,xi; + IssmDouble *xyz_list_base= NULL; + + /*Fetch number of vertices for this finite element*/ + int numvertices = element->GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Build friction element, needed later: */ + Friction* friction=new Friction(element,dim); + + /*Retrieve all inputs we will be needing: */ + element->GetVerticesCoordinatesBase(&xyz_list_base); + element->GradientIndexing(&vertexpidlist[0],control_index); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = element->GetInput(VzEnum); _assert_(vy_input); + Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* adjointz_input = element->GetInput(AdjointzEnum); _assert_(adjointz_input); + Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGaussBase(4); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + adjointx_input->GetInputValue(&lambda, gauss); + adjointy_input->GetInputValue(&mu, gauss); + adjointz_input->GetInputValue(&xi ,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vz_input->GetInputValue(&vz,gauss); + dragcoeff_input->GetInputValue(&drag, gauss); + + friction->GetAlphaComplement(&dalpha2dk,gauss); + element->NormalBase(&normal[0],xyz_list_base); + + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctionsP1(basis,gauss); + + /*Build gradient vector (actually -dJ/dk): */ + for(int i=0;iweight*basis[i]; + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list_base); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; + delete friction; }/*}}}*/ void AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector* gradient,int control_index){/*{{{*/ @@ -1103,25 +1398,24 @@ /*Initialize some vectors*/ IssmDouble* basis = xNew(numvertices); - IssmDouble* ge = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); int* vertexpidlist = xNew(numvertices); /*Retrieve all inputs we will be needing: */ basalelement->GetVerticesCoordinates(&xyz_list); basalelement->GradientIndexing(&vertexpidlist[0],control_index); - 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* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input); - Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input); - Input* rheologyb_input = element->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); + Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); + Input* adjointx_input = basalelement->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input = basalelement->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); /* Start looping on the number of gaussian points: */ Gauss* gauss=basalelement->NewGauss(4); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - thickness_input->GetInputValue(&thickness,gauss); vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); @@ -1152,13 +1446,96 @@ if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ void AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); + + /*WARNING: We use SSA as an estimate for now*/ + this->GradientJBbarSSA(element,gradient,control_index); }/*}}}*/ void AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); + /*WARNING: We use SSA as an estimate for now*/ + this->GradientJBbarSSA(element,gradient,control_index); }/*}}}*/ void AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); + + /*Intermediaries*/ + int domaintype,dim; + Element* basalelement; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + dim = 2; + break; + case Domain2DverticalEnum: + if(!element->IsOnBase()) return; + basalelement = element->SpawnBasalElement(); + dim = 1; + break; + case Domain3DEnum: + if(!element->IsOnBase()) return; + basalelement = element->SpawnBasalElement(); + dim = 2; + break; + default: _error_("mesh "<GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Retrieve all inputs we will be needing: */ + basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->GradientIndexing(&vertexpidlist[0],control_index); + Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); + Input* adjointx_input = basalelement->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input = basalelement->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(4); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + thickness_input->GetInputValue(&thickness,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss); + adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss); + + basalelement->dViscositydDSSA(&dmudD,dim,xyz_list,gauss,vx_input,vy_input); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctionsP1(basis,gauss); + + for(int i=0;iweight*basis[i]; + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ void AdjointHorizAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ int approximation; Index: ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 18059) +++ ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 18060) @@ -150,8 +150,171 @@ _error_("not implemented yet"); }/*}}}*/ void AdjointBalancethicknessAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ - _error_("Not implemented yet"); + /*The gradient of the cost function is calculated in 2 parts. + * + * dJ \partial J \partial lambda^T(KU-F) + * -- = ---------- + ------------------------ + * dk \partial k \parial k + * + * */ + + /*If on water, grad = 0: */ + if(!element->IsIceInElement()) return; + + /*Get list of cost functions*/ + int *responses = NULL; + int num_responses,resp; + element->FindParam(&num_responses,InversionNumCostFunctionsEnum); + element->FindParam(&responses,NULL,InversionCostFunctionsEnum); + + /*Check that control_type is supported*/ + if(control_type!=VxEnum && + control_type!=VyEnum && + control_type!=BalancethicknessThickeningRateEnum){ + _error_("Control "<(responses); + }/*}}}*/ +void AdjointBalancethicknessAnalysis::GradientJVx(Element* element,Vector* gradient,int control_index){/*{{{*/ + + /*Intermediaries*/ + IssmDouble Jdet,weight; + IssmDouble thickness,Dlambda[3],dp[3]; + IssmDouble *xyz_list= NULL; + + /*Fetch number of vertices for this finite element*/ + int numvertices = element->GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(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* adjoint_input = element->GetInput(AdjointEnum); _assert_(adjoint_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(4); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + adjoint_input->GetInputDerivativeValue(&Dlambda[0],xyz_list,gauss); + thickness_input->GetInputValue(&thickness, gauss); + thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsP1(basis,gauss); + + /*Build gradient vector (actually -dJ/dD): */ + for(int i=0;iweight*basis[i]; + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; +}/*}}}*/ +void AdjointBalancethicknessAnalysis::GradientJVy(Element* element,Vector* gradient,int control_index){/*{{{*/ + + /*Intermediaries*/ + IssmDouble Jdet,weight; + IssmDouble thickness,Dlambda[3],dp[3]; + IssmDouble *xyz_list= NULL; + + /*Fetch number of vertices for this finite element*/ + int numvertices = element->GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(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* adjoint_input = element->GetInput(AdjointEnum); _assert_(adjoint_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(4); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + adjoint_input->GetInputDerivativeValue(&Dlambda[0],xyz_list,gauss); + thickness_input->GetInputValue(&thickness, gauss); + thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsP1(basis,gauss); + + /*Build gradient vector (actually -dJ/dvy): */ + for(int i=0;iweight*basis[i]; + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; +}/*}}}*/ +void AdjointBalancethicknessAnalysis::GradientJDhDt(Element* element,Vector* gradient,int control_index){/*{{{*/ + + /*Fetch number of vertices for this finite element*/ + int numvertices = element->GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* ge = xNewZeroInit(numvertices); + IssmDouble* lambda = xNew(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Retrieve all inputs we will be needing: */ + element->GradientIndexing(&vertexpidlist[0],control_index); + element->GetInputListOnVertices(lambda,AdjointEnum); + for(int i=0;i(ge[i])); + } + gradient->SetValues(numvertices,vertexpidlist,ge,INS_VAL); + + /*Clean up and return*/ + xDelete(ge); + xDelete(lambda); + xDelete(vertexpidlist); +}/*}}}*/ void AdjointBalancethicknessAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ int domaintype; Index: ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h =================================================================== --- ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h (revision 18059) +++ ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h (revision 18060) @@ -27,6 +27,9 @@ ElementVector* CreatePVector(Element* element); void GetSolutionFromInputs(Vector* solution,Element* element); void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void GradientJVx(Element* element,Vector* gradient,int control_index); + void GradientJVy(Element* element,Vector* gradient,int control_index); + void GradientJDhDt(Element* element,Vector* gradient,int control_index); void InputUpdateFromSolution(IssmDouble* solution,Element* element); void UpdateConstraints(FemModel* femmodel); }; Index: ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 18059) +++ ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 18060) @@ -8,8 +8,8 @@ void Gradjx(Vector** pgradient,IssmDouble** pnorm_list, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){ - int i,j,numberofvertices; - int num_controls; + int numberofvertices; + int num_controls,analysisenum; IssmDouble norm_inf; IssmDouble *norm_list = NULL; int *control_type = NULL; @@ -21,38 +21,39 @@ parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); numberofvertices=vertices->NumberOfVertices(); + /*Get current analysis*/ + parameters->FindParam(&analysisenum,AnalysisTypeEnum); + Analysis* analysis = EnumToAnalysis(analysisenum); + /*Allocate gradient_list */ gradient_list = xNew*>(num_controls); norm_list = xNew(num_controls); - for(i=0;i(num_controls*numberofvertices); } gradient=new Vector(num_controls*numberofvertices); /*Compute all gradient_list*/ - for(i=0;iSize();j++){ - + for(int i=0;iSize();j++){ Element* element=(Element*)elements->GetObjectByOffset(j); - element->Gradj(gradient_list[i],control_type[i],i); + //element->Gradj(gradient_list[i],control_type[i],i); + analysis->GradientJ(gradient_list[i],element,control_type[i],i); } - gradient_list[i]->Assemble(); - norm_list[i]=gradient_list[i]->Norm(NORM_INF); } /*Add all gradient_list together*/ - for(i=0;iAXPY(gradient_list[i],1.0); delete gradient_list[i]; } /*Check that gradient is clean*/ norm_inf=gradient->Norm(NORM_INF); - if(norm_inf<=0) _error_("||∂J/∂α||∞ = 0 gradient norm is zero"); - if(xIsNan(norm_inf))_error_("||∂J/∂α||∞ = NaN gradient norm is NaN"); + if(norm_inf<=0) _error_("||dJ/dk|| = 0 gradient norm is zero"); + if(xIsNan(norm_inf))_error_("||dJ/dk|| = NaN gradient norm is NaN"); /*Clean-up and assign output pointer*/ if(pnorm_list){ Index: ../trunk-jpl/src/c/modules/Gradjx/Gradjx.h =================================================================== --- ../trunk-jpl/src/c/modules/Gradjx/Gradjx.h (revision 18059) +++ ../trunk-jpl/src/c/modules/Gradjx/Gradjx.h (revision 18060) @@ -6,6 +6,7 @@ #define _GRADJX_H #include "../../classes/classes.h" +#include "../../analyses/analyses.h" /* local prototypes: */ void Gradjx(Vector** pgrad_g,IssmDouble** pgrad_norm,Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters); Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18059) +++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18060) @@ -284,21 +284,10 @@ } /*}}}*/ /*FUNCTION Matice::GetViscosity_B {{{*/ -void Matice::GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){ - /*From a string tensor and a material object, return viscosity, using Glen's flow law. - (1-D) - viscosity= ----------------------- - 2 eps_eff ^[(n-1)/n] +void Matice::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){ - where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate - and n the flow law exponent. - - If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we - return 10^14, initial viscosity. - */ - /*output: */ - IssmDouble viscosity; + IssmDouble dmudB; /*Intermediary: */ IssmDouble D=0.,n; @@ -310,27 +299,44 @@ _assert_(D>=0. && D<1.); } - if (n==1.){ - /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */ - viscosity=(1.-D)/2.; + if(n==1.){ + /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */ + dmudB=(1.-D)/2.; } else{ + if(eps_eff==0.) dmudB = 0.; + else dmudB = (1.-D)/(2.*pow(eps_eff,(n-1.)/n)); + } - /*if no strain rate, return maximum viscosity*/ - if(eps_eff==0.){ - viscosity = 1.e+14/2.; - } + /*Return: */ + *pdmudB=dmudB; +} +/*}}}*/ +/*FUNCTION Matice::GetViscosity_D {{{*/ +void Matice::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){ - else{ - viscosity=(1.-D)/(2.*pow(eps_eff,(n-1.)/n)); - } + /*output: */ + IssmDouble dmudD; + + /*Intermediary: */ + IssmDouble n,B; + + /*Get B and n*/ + n=GetN(); _assert_(n>0.); + B=GetBbar(); + _assert_(this->isdamaged); + + if(n==1.){ + /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */ + dmudD=-B/2.; } + else{ + if(eps_eff==0.) dmudD = 0.; + else dmudD = -B/(2.*pow(eps_eff,(n-1.)/n)); + } - /*Checks in debugging mode*/ - if(viscosity<=0) _error_("Negative viscosity"); - /*Return: */ - *pviscosity=viscosity; + *pdmudD=dmudD; } /*}}}*/ /*FUNCTION Matice::GetViscosityBar {{{*/ Index: ../trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Material.h (revision 18059) +++ ../trunk-jpl/src/c/classes/Materials/Material.h (revision 18060) @@ -26,6 +26,7 @@ virtual void Configure(Elements* elements)=0; virtual void GetViscosity(IssmDouble* pviscosity,IssmDouble epseff)=0; virtual void GetViscosity_B(IssmDouble* pviscosity,IssmDouble epseff)=0; + virtual void GetViscosity_D(IssmDouble* pviscosity,IssmDouble epseff)=0; virtual void GetViscosityBar(IssmDouble* pviscosity,IssmDouble epseff)=0; virtual void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0; virtual void GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0; Index: ../trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 18059) +++ ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 18060) @@ -55,6 +55,7 @@ void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); void GetViscosity(IssmDouble* pviscosity, IssmDouble eps_eff); void GetViscosity_B(IssmDouble* pviscosity, IssmDouble eps_eff); + void GetViscosity_D(IssmDouble* pviscosity, IssmDouble eps_eff); void GetViscosityBar(IssmDouble* pviscosity, IssmDouble eps_eff); void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); void GetViscosityDComplement(IssmDouble*, IssmDouble*); Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 18059) +++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 18060) @@ -76,6 +76,7 @@ void Configure(Elements* elements); void GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");}; void GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");}; + void GetViscosity_D(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");}; void GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");}; void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");}; void GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");}; Index: ../trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 18059) +++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 18060) @@ -60,6 +60,7 @@ void DeepEcho(); void DeleteMaterials(void); void dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); IssmDouble Divergence(void); void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure); void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure); Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18059) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18060) @@ -1456,6 +1456,13 @@ } /*}}}*/ +void Tria::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetNodalFunctions(basis,(GaussTria*)gauss,P1Enum); + +} +/*}}}*/ void Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ _assert_(gauss->Enum()==GaussTriaEnum); @@ -1463,6 +1470,13 @@ } /*}}}*/ +void Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum); + +} +/*}}}*/ void Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ _assert_(gauss->Enum()==GaussTriaEnum); Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18059) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18060) @@ -190,9 +190,9 @@ Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");}; Gauss* NewGaussTop(int order); void NodalFunctions(IssmDouble* basis,Gauss* gauss); - void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; + void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss); void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); - void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; + void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); void NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss); Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18059) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18060) @@ -2999,6 +2999,7 @@ } /*}}}*/ 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]; Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18059) +++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18060) @@ -239,6 +239,33 @@ } /*}}}*/ +void Element::dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dmudB; + IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ + IssmDouble epsilon1d; /* epsilon=[exx]; */ + IssmDouble eps_eff; + + if(dim==2){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/ + this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]); + } + else{ + /* eps_eff^2 = 1/2 exx^2*/ + this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); + eps_eff = sqrt(epsilon1d*epsilon1d/2.); + } + + /*Get viscosity*/ + material->GetViscosity_D(&dmudB,eps_eff); + + /*Assign output pointer*/ + *pdmudB=dmudB; + +} +/*}}}*/ void Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/ matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure); }/*}}}*/ Index: ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 18059) +++ ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 18060) @@ -3,7 +3,7 @@ */ /*Headers:*/ -/*{{{*//*{{{*/ +/*{{{*/ #ifdef HAVE_CONFIG_H #include #else