Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp (revision 6130) +++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp (revision 6131) @@ -39,5 +39,5 @@ IoModelFetchData(&iomodel->thickness_obs,NULL,NULL,iomodel_handle,"thickness_obs"); IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights"); - if(iomodel->control_type==DhDtEnum){ + if(iomodel->control_type==DhDtEnum || iomodel->control_type==VxEnum || iomodel->control_type==VyEnum){ IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,EnumToModelField(iomodel->control_type)); //copy the control parameter in iomodel } Index: /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 6130) +++ /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 6131) @@ -26,5 +26,10 @@ /*What control type?*/ - if (iomodel->control_type!=DragCoefficientEnum && iomodel->control_type!=RheologyBbarEnum && iomodel->control_type!=DhDtEnum){ + if (iomodel->control_type!=DragCoefficientEnum && + iomodel->control_type!=RheologyBbarEnum && + iomodel->control_type!=DhDtEnum && + iomodel->control_type!=VxEnum && + iomodel->control_type!=VyEnum + ){ ISSMERROR("control_type %s not supported yet!",EnumToString(iomodel->control_type)); } Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6130) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6131) @@ -673,15 +673,17 @@ } else if (control_type==RheologyBbarEnum){ - Input* Bbar_input=matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(Bbar_input); - Bbar_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss); - //Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight; + //nothing } else if (control_type==DhDtEnum){ - Input* dhdt_input=inputs->GetInput(DhDtEnum); ISSMASSERT(dhdt_input); - dhdt_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss); - //Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight; + //nothing + } + else if (control_type==VxEnum){ + //nothing + } + else if (control_type==VyEnum){ + //nothing } else{ - ISSMERROR("%s%i","unsupported control type: ",control_type); + ISSMERROR("unsupported control type: %s",EnumToString(control_type)); } } @@ -882,4 +884,10 @@ GradjDhDt(gradient); break; + case VxEnum: + GradjVx(gradient); + break; + case VyEnum: + GradjVy(gradient); + break; default: ISSMERROR("%s%i","control type not supported yet: ",control_type); @@ -964,15 +972,10 @@ int i; - - /* node data: */ + int ig; + double xyz_list[NUMVERTICES][3]; int doflist1[NUMVERTICES]; double dh1dh3[NDOF2][NUMVERTICES]; - - /* gaussian points: */ - int ig; GaussTria *gauss=NULL; - - /* parameters: */ double dk[NDOF2]; double vx,vy; @@ -983,19 +986,9 @@ double drag; Friction* friction=NULL; - - /*element vector at the gaussian points: */ double grade_g[NUMVERTICES]={0.0}; double grade_g_gaussian[NUMVERTICES]; - - /* Jacobian: */ double Jdet; - - /*nodal functions: */ double l1l2l3[3]; - - /* strain rate: */ double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - - /*parameters: */ double cm_noisedmp; int analysis_type; @@ -1100,4 +1093,88 @@ } /*}}}*/ +/*FUNCTION Tria::GradjVx{{{1*/ +void Tria::GradjVx(Vec gradient){ + + /*Intermediaries*/ + int i,ig; + int doflist1[NUMVERTICES]; + double l1l2l3[3]; + double thickness,Jdet; + double Dlambda[2]; + double xyz_list[NUMVERTICES][3]; + GaussTria *gauss = NULL; + double grade_g[NUMVERTICES] = {0.0}; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GetDofList1(&doflist1[0]); + + /*Retrieve all inputs we will be needing: */ + Input* adjoint_input=inputs->GetInput(AdjointEnum); ISSMASSERT(adjoint_input); + Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); + thickness_input->GetParameterValue(&thickness, gauss); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(l1l2l3, gauss); + + for(i=0;iweight*l1l2l3[i]; + } + + VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Tria::GradjVy{{{1*/ +void Tria::GradjVy(Vec gradient){ + + /*Intermediaries*/ + int i,ig; + int doflist1[NUMVERTICES]; + double thickness,Jdet; + double l1l2l3[3]; + double Dlambda[2]; + double xyz_list[NUMVERTICES][3]; + GaussTria *gauss = NULL; + double grade_g[NUMVERTICES] = {0.0}; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GetDofList1(&doflist1[0]); + + /*Retrieve all inputs we will be needing: */ + Input* adjoint_input=inputs->GetInput(AdjointEnum); ISSMASSERT(adjoint_input); + Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); + thickness_input->GetParameterValue(&thickness, gauss); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(l1l2l3, gauss); + + for(i=0;iweight*l1l2l3[i]; + } + + VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ /*FUNCTION Tria::InputControlUpdate{{{1*/ void Tria::InputControlUpdate(double scalar,bool save_parameter){ @@ -1114,58 +1191,61 @@ /*Rheology*/ - if(control_type==RheologyBbarEnum){ - - /*First, get revert to previous parameter value (kept in ControlParameter input)*/ - matice->inputs->DuplicateInput(ControlParameterEnum,RheologyBbarEnum); - - /*Now, update using the gradient new = old + scalar * gradient*/ - //matice->inputs->AXPY(RheologyBbarEnum,scalar,GradientEnum); - // For now: Gradient is in element (TO BE CHANGED) and parameter in matice - Input* input_B =(Input*)matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(input_B); - Input* input_Grad=(Input*)this->inputs->GetInput(GradientEnum); ISSMASSERT(input_Grad); - input_B->AXPY(input_Grad,scalar); - - /*Constrain constrain input*/ - input=(Input*)matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(input); - input->Constrain(cm_min,cm_max); - - /*Finally: save input if requested*/ - if (save_parameter) matice->inputs->DuplicateInput(RheologyBbarEnum,ControlParameterEnum); - - } - else if(control_type==DragCoefficientEnum){ - - /*First, get revert to previous parameter value (kept in ControlParameter input)*/ - this->inputs->DuplicateInput(ControlParameterEnum,DragCoefficientEnum); - - /*Now, update using the gradient new = old + scalar * gradient*/ - this->inputs->AXPY(DragCoefficientEnum,scalar,GradientEnum); - - /*Constrain input*/ - input=(Input*)this->inputs->GetInput(DragCoefficientEnum); ISSMASSERT(input); - input->Constrain(cm_min,cm_max); - - /*Finally: save input if requested*/ - if (save_parameter) inputs->DuplicateInput(DragCoefficientEnum,ControlParameterEnum); - - } - else if(control_type==DhDtEnum){ - - /*First, get revert to previous parameter value (kept in ControlParameter input)*/ - this->inputs->DuplicateInput(ControlParameterEnum,DhDtEnum); - - /*Now, update using the gradient new = old + scalar * gradient*/ - this->inputs->AXPY(DhDtEnum,scalar,GradientEnum); - - /*Constrain input*/ - input=(Input*)this->inputs->GetInput(DhDtEnum); ISSMASSERT(input); - input->Constrain(cm_min,cm_max); - - /*Finally: save input if requested*/ - if (save_parameter) inputs->DuplicateInput(DhDtEnum,ControlParameterEnum); - - } - else{ - ISSMERROR("control type %s not implemented yet",EnumToString(control_type)); + switch(control_type){ + + case RheologyBbarEnum:{ + + /*First, get revert to previous parameter value (kept in ControlParameter input)*/ + matice->inputs->DuplicateInput(ControlParameterEnum,RheologyBbarEnum); + + /*Now, update using the gradient new = old + scalar * gradient*/ + //matice->inputs->AXPY(RheologyBbarEnum,scalar,GradientEnum); + // For now: Gradient is in element (TO BE CHANGED) and parameter in matice + Input* input_B =(Input*)matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(input_B); + Input* input_Grad=(Input*)this->inputs->GetInput(GradientEnum); ISSMASSERT(input_Grad); + input_B->AXPY(input_Grad,scalar); + + /*Constrain constrain input*/ + input=(Input*)matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(input); + input->Constrain(cm_min,cm_max); + + /*Finally: save input if requested*/ + if (save_parameter) matice->inputs->DuplicateInput(RheologyBbarEnum,ControlParameterEnum); + } + break; + + case DragCoefficientEnum:{ + this->inputs->DuplicateInput(ControlParameterEnum,DragCoefficientEnum); + this->inputs->AXPY(DragCoefficientEnum,scalar,GradientEnum); + input=(Input*)this->inputs->GetInput(DragCoefficientEnum); ISSMASSERT(input); + input->Constrain(cm_min,cm_max); + if (save_parameter) inputs->DuplicateInput(DragCoefficientEnum,ControlParameterEnum); + } + break; + case DhDtEnum:{ + this->inputs->DuplicateInput(ControlParameterEnum,DhDtEnum); + this->inputs->AXPY(DhDtEnum,scalar,GradientEnum); + input=(Input*)this->inputs->GetInput(DhDtEnum); ISSMASSERT(input); + input->Constrain(cm_min,cm_max); + if (save_parameter) inputs->DuplicateInput(DhDtEnum,ControlParameterEnum); + } + break; + case VxEnum:{ + this->inputs->DuplicateInput(ControlParameterEnum,VxEnum); + this->inputs->AXPY(VxEnum,scalar,GradientEnum); + input=(Input*)this->inputs->GetInput(VxEnum); ISSMASSERT(input); + input->Constrain(cm_min,cm_max); + if (save_parameter) inputs->DuplicateInput(VxEnum,ControlParameterEnum); + } + break; + case VyEnum:{ + this->inputs->DuplicateInput(ControlParameterEnum,VyEnum); + this->inputs->AXPY(VyEnum,scalar,GradientEnum); + input=(Input*)this->inputs->GetInput(VyEnum); ISSMASSERT(input); + input->Constrain(cm_min,cm_max); + if (save_parameter) inputs->DuplicateInput(VyEnum,ControlParameterEnum); + } + break; + default: + ISSMERROR("control type %s not implemented yet",EnumToString(control_type)); } @@ -3894,5 +3974,5 @@ * 2 [ obs obs ] * - * dJ 2 + * dJ * DU = - -- = (u - u ) * du obs Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 6130) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 6131) @@ -86,4 +86,6 @@ void GradjDrag(Vec gradient); void GradjDhDt(Vec gradient); + void GradjVx(Vec gradient); + void GradjVy(Vec gradient); void InputControlUpdate(double scalar,bool save_parameter); void InputArtificialNoise(int enum_type,double min, double max); Index: /issm/trunk/src/c/solutions/adjointbalancedthickness_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/adjointbalancedthickness_core.cpp (revision 6130) +++ /issm/trunk/src/c/solutions/adjointbalancedthickness_core.cpp (revision 6131) @@ -38,7 +38,7 @@ /*Save results*/ - if(solution_type==AdjointSolutionEnum && !control_analysis){ + //if(solution_type==AdjointSolutionEnum && !control_analysis){ if(verbose)_printf_("saving results:\n"); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointEnum); - } + //} }