Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8646) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8647) @@ -184,5 +184,5 @@ CmResponseEnum, CmResponsesEnum, - CmNoiseDmpEnum, + CmNoiseDmpEnum, //TO BE DELETED (not used anymore) ConstantEnum, NumControlsEnum, Index: /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 8646) +++ /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 8647) @@ -43,5 +43,4 @@ parameters->AddObject(new DoubleParam(EpsCmEnum,iomodel->eps_cm)); parameters->AddObject(new DoubleParam(MeanVelEnum,iomodel->meanvel)); - parameters->AddObject(new DoubleParam(CmNoiseDmpEnum,iomodel->cm_noisedmp)); parameters->AddObject(new BoolParam(CmGradientEnum,iomodel->cm_gradient)); Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 8646) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 8647) @@ -50,5 +50,4 @@ virtual double RheologyBbarAbsGradient(bool process_units,int weight_index)=0; virtual double DragCoefficientAbsGradient(bool process_units,int weight_index)=0; - virtual double RegularizeInversion(void)=0; virtual double SurfaceArea(void)=0; virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8646) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8647) @@ -3986,4 +3986,10 @@ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. tria->GradjDragGradient(gradient,resp); + delete tria->matice; delete tria; + break; + case RheologyBbarAbsGradientEnum: + if (!IsOnBed()) return; + tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + tira->GradjBGradient(gradient,resp); delete tria->matice; delete tria; break; @@ -6724,61 +6730,4 @@ /*Affect value to the reduced matrix */ for(i=0;i<24;i++) *(Pe_reduced+i)=Pi[i]-Pright[i]; -} -/*}}}*/ -/*FUNCTION Penta::RegularizeInversion {{{1*/ -double Penta::RegularizeInversion(void){ - - double J; - Tria* tria=NULL; - - /*recover some inputs: */ - int approximation; - inputs->GetParameterValue(&approximation,ApproximationEnum); - - /*If on water, return 0: */ - if(IsOnWater())return 0; - - /*Bail out if this element if: - * -> Not MacAyeal and not on the surface - * -> MacAyeal (2d model) and not on bed) */ - if ((approximation!=MacAyealApproximationEnum && !IsOnSurface()) || (approximation==MacAyealApproximationEnum && !IsOnBed())){ - return 0; - } - else if (approximation==MacAyealApproximationEnum){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute RegularizeInversion*/ - - /*Depth Average B*/ - this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum); - - tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). - J=tria->RegularizeInversion(); - delete tria->matice; delete tria; - - /*delete B average*/ - this->matice->inputs->DeleteInput(RheologyBbarEnum); - - return J; - } - else{ - - /*Depth Average B and put it in inputs*/ - Penta* penta=GetBasalElement(); - penta->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum); - Input* B_input=penta->matice->inputs->GetInput(RheologyBbarEnum); - Input* B_copy=(Input*)B_input->copy(); - this->matice->inputs->AddInput((Input*)B_copy); - - tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). - J=tria->RegularizeInversion(); - delete tria->matice; delete tria; - - /*delete B average*/ - this->matice->inputs->DeleteInput(RheologyBbarEnum); - penta->matice->inputs->DeleteInput(RheologyBbarEnum); - - return J; - } } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 8646) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 8647) @@ -79,5 +79,4 @@ void Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); - double RegularizeInversion(void); void CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg,Vec df); void CreatePVector(Vec pg, Vec pf); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8646) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8647) @@ -2930,7 +2930,51 @@ GradjDragGradient(gradient,resp); break; + case RheologyBbarAbsGradientEnum: + GradjBGradient(gradient,resp); + break; default: _error_("response %s not supported yet",EnumToStringx(responses[resp])); } +} +/*}}}*/ +/*FUNCTION Tria::GradjBGradient{{{1*/ +void Tria::GradjBGradient(Vec gradient, int weight_index){ + + int i,ig; + int doflist1[NUMVERTICES]; + double Jdet,weight; + double xyz_list[NUMVERTICES][3]; + double dh1dh3[NDOF2][NUMVERTICES]; + double dk[NDOF2]; + double grade_g[NUMVERTICES]={0.0}; + GaussTria *gauss=NULL; + + /*Retrieve all inputs we will be needing: */ + if(IsOnShelf())return; + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GetDofList1(&doflist1[0]); + Input* rheologyb_input=matice->inputs->GetInput(RheologyBbarEnum); _assert_(rheologyb_input); + Input* weights_input=inputs->GetInput(WeightsEnum); _assert_(weights_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss); + weights_input->GetParameterValue(&weight,gauss,weight_index); + + /*Build alpha_complement_list: */ + rheologyb_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;iweight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); + } + VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); + + /*Clean up and return*/ + delete gauss; } /*}}}*/ @@ -2942,15 +2986,11 @@ int doflist[NUMVERTICES]; double vx,vy,lambda,mu,thickness,Jdet; - double cm_noisedmp,viscosity_complement; + double viscosity_complement; double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2]; double xyz_list[NUMVERTICES][3]; double basis[3],epsilon[3]; double dbasis[NDOF2][NUMVERTICES]; - double grad_g[NUMVERTICES]; double grad[NUMVERTICES]={0.0}; GaussTria *gauss = NULL; - - /*retrieve some parameters: */ - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); /* Get node coordinates and dof list: */ @@ -2987,10 +3027,7 @@ /*standard gradient dJ/dki*/ - for (i=0;iweight*basis[i]; - } - /*Add regularization term*/ - for (i=0;iweight*(dbasis[0][i]*dB[0]+dbasis[1][i]*dB[1]); - for(i=0;iweight*basis[i]; } @@ -3008,5 +3045,5 @@ int doflist1[NUMVERTICES]; double vx,vy,lambda,mu,alpha_complement,Jdet; - double bed,thickness,Neff,drag,cm_noisedmp; + double bed,thickness,Neff,drag; double xyz_list[NUMVERTICES][3]; double dh1dh3[NDOF2][NUMVERTICES]; @@ -3019,10 +3056,8 @@ GaussTria *gauss=NULL; + if(IsOnShelf())return; + /*retrive parameters: */ parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve some parameters ands return if iceshelf: */ - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - if(IsOnShelf())return; GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); GetDofList1(&doflist1[0]); @@ -3062,10 +3097,5 @@ /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ for (i=0;iweight*l1l2l3[i]; - - //noise dampening d/dki(1/2*(dk/dx)^2) - grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); } @@ -3092,5 +3122,5 @@ double bed,thickness,Neff; double lambda,mu,xi,Jdet,vx,vy,vz; - double alpha_complement,drag,cm_noisedmp; + double alpha_complement,drag; double surface_normal[3],bed_normal[3]; double xyz_list[NUMVERTICES][3]; @@ -3117,7 +3147,4 @@ Input* adjointz_input=inputs->GetInput(AdjointzEnum); _assert_(adjointz_input); - /*retrieve some parameters: */ - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - /*Get out if shelf*/ if(IsOnShelf())return; @@ -3179,7 +3206,4 @@ -xi *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2])) )*Jdet*gauss->weight*l1l2l3[i]; - - //Add regularization term - grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); } @@ -3261,5 +3285,5 @@ int i,ig; int doflist1[NUMVERTICES]; - double thickness,Jdet,cm_noisedmp; + double thickness,Jdet; double l1l2l3[3]; double dbasis[NDOF2][NUMVERTICES]; @@ -3274,5 +3298,4 @@ /*Retrieve all inputs we will be needing: */ - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); Input* adjoint_input=inputs->GetInput(AdjointEnum); _assert_(adjoint_input); Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); @@ -3293,7 +3316,4 @@ for(i=0;iweight*l1l2l3[i]; - - //noise dampening d/dki(1/2*(dk/dx)^2) - //for (i=0;iweight*(dbasis[0][i]*dp[0]+dbasis[1][i]*dp[1]); } @@ -3310,5 +3330,5 @@ int i,ig; int doflist1[NUMVERTICES]; - double thickness,Jdet,cm_noisedmp; + double thickness,Jdet; double l1l2l3[3]; double dbasis[NDOF2][NUMVERTICES]; @@ -3319,5 +3339,4 @@ /* Get node coordinates and dof list: */ - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); GetDofList1(&doflist1[0]); @@ -3342,9 +3361,5 @@ for(i=0;iweight*l1l2l3[i]; - - //noise dampening d/dki(1/2*(dk/dx)^2) - //for (i=0;iweight*(dbasis[0][i]*dp[0]+dbasis[1][i]*dp[1]); - } - + } VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); @@ -4729,76 +4744,4 @@ } /*}}}*/ -/*FUNCTION Tria::RegularizeInversion {{{1*/ -double Tria::RegularizeInversion(){ - - /*constants: */ - const int numdof=NDOF2*NUMVERTICES; - - /* Intermediaries */ - int i,j,ig; - int num_controls; - int *control_type=NULL; - double Jelem = 0; - double cm_noisedmp; - double Jdet; - double xyz_list[NUMVERTICES][3]; - double dk[NDOF2],dB[NDOF2]; - GaussTria *gauss = NULL; - - /*retrieve parameters and inputs*/ - this->parameters->FindParam(&num_controls,NumControlsEnum); - this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - - /*If on water, return 0: */ - if(IsOnWater()) return 0; - - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /*Add Tikhonov regularization term to misfit: - * - * WARNING: for now, the regularization is only taken into account by the gradient - * the misfit reflects the response only! - * - * */ - if (control_type[i]==DragCoefficientEnum){ - Input* drag_input=inputs->GetInput(DragCoefficientEnum); _assert_(drag_input); - drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - //Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight; - } - else if (control_type[i]==RheologyBbarEnum){ - //nothing - } - else if (control_type[i]==DhDtEnum){ - //nothing - } - else if (control_type[i]==VxEnum){ - //nothing - } - else if (control_type[i]==VyEnum){ - //nothing - } - else{ - _error_("unsupported control type: %s",EnumToStringx(control_type[i])); - } - } - - delete gauss; - } - - /*Clean up and return*/ - xfree((void**)&control_type); - return Jelem; -} -/*}}}*/ /*FUNCTION Tria::RheologyBbarAbsGradient{{{1*/ double Tria::RheologyBbarAbsGradient(bool process_units,int weight_index){ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 8646) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 8647) @@ -76,5 +76,4 @@ void Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); - double RegularizeInversion(void); void CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg,Vec df); void CreatePVector(Vec pg, Vec pf); Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 8646) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 8647) @@ -181,5 +181,4 @@ IoModelFetchData(&this->eps_cm,iomodel_handle,"eps_cm"); IoModelFetchData(&this->tolx,iomodel_handle,"tolx"); - IoModelFetchData(&this->cm_noisedmp,iomodel_handle,"cm_noisedmp"); IoModelFetchData(&this->cm_gradient,iomodel_handle,"cm_gradient"); IoModelFetchData(&this->eps_res,iomodel_handle,"eps_res"); @@ -354,5 +353,4 @@ this->tolx=0; this->maxiter=NULL; - this->cm_noisedmp=0; this->cm_min=NULL; this->cm_max=NULL; Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 8646) +++ /issm/trunk/src/c/objects/IoModel.h (revision 8647) @@ -125,10 +125,7 @@ double stokesreconditioning; int shelf_dampening; - double cm_noisedmp; double* cm_min; double* cm_max; int cm_gradient;; - - double cm_noisedampening; /*control methods: */