Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 18059) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 18060) @@ -151,5 +151,168 @@ }/*}}}*/ 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){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h (revision 18059) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h (revision 18060) @@ -28,4 +28,7 @@ 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: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18059) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18060) @@ -923,5 +923,7 @@ /*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; @@ -1012,5 +1014,5 @@ /*Initialize some vectors*/ IssmDouble* dbasis = xNew(2*numvertices); - IssmDouble* ge = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); int* vertexpidlist = xNew(numvertices); @@ -1056,16 +1058,4 @@ }/*}}}*/ void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -void AdjointHorizAnalysis::GradientJDragHO(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -void AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -void AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector* gradient,int control_index){/*{{{*/ /*Intermediaries*/ @@ -1095,6 +1085,5 @@ /*Intermediaries*/ IssmDouble Jdet,weight; - IssmDouble thickness,dmudB; - IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3]; + IssmDouble dk[3]; IssmDouble *xyz_list= NULL; @@ -1103,6 +1092,6 @@ /*Initialize some vectors*/ - IssmDouble* basis = xNew(numvertices); - IssmDouble* ge = xNew(numvertices); + IssmDouble* dbasis = xNew(2*numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); int* vertexpidlist = xNew(numvertices); @@ -1110,10 +1099,94 @@ 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* 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){/*{{{*/ + + /*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: */ @@ -1122,4 +1195,225 @@ 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){/*{{{*/ + + /*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){/*{{{*/ + + /*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){/*{{{*/ + + /*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); @@ -1153,11 +1447,94 @@ }/*}}}*/ 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){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18059) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18060) @@ -234,4 +234,31 @@ /*Get viscosity*/ material->GetViscosity_B(&dmudB,eps_eff); + + /*Assign output pointer*/ + *pdmudB=dmudB; + +} +/*}}}*/ +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*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18059) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18060) @@ -61,4 +61,5 @@ 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); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18059) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18060) @@ -3000,4 +3000,5 @@ void Penta::GradjDragHO(Vector* gradient,int control_index){/*{{{*/ + printf("-------------- file: Penta.cpp line: %i\n",__LINE__); int i,j; int analysis_type; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18059) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18060) @@ -1457,8 +1457,22 @@ } /*}}}*/ +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); this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss); + +} +/*}}}*/ +void Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum); } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18059) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18060) @@ -191,7 +191,7 @@ 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); Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 18059) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 18060) @@ -4,5 +4,5 @@ /*Headers:*/ -/*{{{*//*{{{*/ +/*{{{*/ #ifdef HAVE_CONFIG_H #include Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 18059) +++ /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 18060) @@ -27,4 +27,5 @@ 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; Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18059) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18060) @@ -285,19 +285,8 @@ /*}}}*/ /*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] - - 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; +void Matice::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){ + + /*output: */ + IssmDouble dmudB; /*Intermediary: */ @@ -311,25 +300,42 @@ } - if (n==1.){ - /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */ - viscosity=(1.-D)/2.; - } - else{ - - /*if no strain rate, return maximum viscosity*/ - if(eps_eff==0.){ - viscosity = 1.e+14/2.; - } - - else{ - viscosity=(1.-D)/(2.*pow(eps_eff,(n-1.)/n)); - } - } - - /*Checks in debugging mode*/ - if(viscosity<=0) _error_("Negative viscosity"); + 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)); + } /*Return: */ - *pviscosity=viscosity; + *pdmudB=dmudB; +} +/*}}}*/ +/*FUNCTION Matice::GetViscosity_D {{{*/ +void Matice::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){ + + /*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)); + } + + /*Return: */ + *pdmudD=dmudD; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 18059) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 18060) @@ -56,4 +56,5 @@ 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); Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 18059) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 18060) @@ -77,4 +77,5 @@ 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");}; Index: /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 18059) +++ /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 18060) @@ -9,6 +9,6 @@ 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; @@ -22,8 +22,12 @@ 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); } @@ -31,19 +35,16 @@ /*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]; @@ -52,6 +53,6 @@ /*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*/ Index: /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h (revision 18059) +++ /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.h (revision 18060) @@ -7,4 +7,5 @@ #include "../../classes/classes.h" +#include "../../analyses/analyses.h" /* local prototypes: */