Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 17970) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 17971) @@ -332,5 +332,6 @@ /*Prepare coordinate system list*/ int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iFindParam(&dim,DomainDimensionEnum); + /*Fetch number of nodes and dof for this finite element*/ int vnumnodes = element->NumberofNodesVelocity(); int pnumnodes = element->NumberofNodesPressure(); - int vnumdof = vnumnodes*3; + int vnumdof = vnumnodes*dim; int pnumdof = pnumnodes*1; @@ -970,5 +972,6 @@ int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iFindParam(&num_controls,InversionNumControlParametersEnum); + + /*get gradient indices*/ + for(int i=0;ivertices[i]->Pid() + control_index; + } + +} +/*}}}*/ +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; +} +/*}}}*/ bool Seg::IsIcefront(void){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17970) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17971) @@ -169,5 +169,5 @@ IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");}; - void GradientIndexing(int* indexing,int control_index){_error_("not implemented yet");}; + void GradientIndexing(int* indexing,int control_index); void Gradj(Vector* gradient,int control_type,int control_index){_error_("not implemented yet");}; void GradjBGradient(Vector* gradient,int control_index){_error_("not implemented yet");}; @@ -176,5 +176,5 @@ 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 GradjDragFS(Vector* gradient,int control_index); void GradjDragGradient(Vector* gradient,int control_index){_error_("not implemented yet");}; void GradjDhDtBalancedthickness(Vector* gradient,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 17970) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17971) @@ -2863,4 +2863,7 @@ /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/ + int approximation; + Seg* seg=NULL; + /*If on water, grad = 0: */ if(!IsIceInElement()) return; @@ -2869,4 +2872,19 @@ 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; GradjDragSSA(gradient,control_index); break; @@ -3163,4 +3181,16 @@ 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; } /*}}}*/ Index: /issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp (revision 17970) +++ /issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp (revision 17971) @@ -38,6 +38,8 @@ if(VerboseSolution()) _printf0_(" saving results\n"); if (isFS){ - int outputs[4] = {AdjointxEnum,AdjointyEnum,AdjointzEnum,AdjointpEnum}; - femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4); + //int outputs[4] = {AdjointxEnum,AdjointyEnum,AdjointzEnum,AdjointpEnum}; + //femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4); + int outputs[3] = {AdjointxEnum,AdjointyEnum,AdjointpEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3); } else{