Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 17971) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 17972) @@ -195,4 +195,47 @@ } /*}}}*/ +void Seg::GradjDragGradient(Vector* gradient,int control_index){/*{{{*/ + + int i; + int vertexpidlist[NUMVERTICES]; + IssmDouble Jdet,weight; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dbasis[NDOF2][NUMVERTICES]; + IssmDouble dk[NDOF2]; + IssmDouble grade_g[NUMVERTICES]={0.0}; + GaussSeg *gauss=NULL; + + /*Retrieve all inputs we will be needing: */ + if(IsFloating())return; + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + GradientIndexing(&vertexpidlist[0],control_index); + Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input); + Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussSeg(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); + weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum); + + /*Build alpha_complement_list: */ + dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;iweight*dbasis[0][i]*dk[0]; + _assert_(!xIsNan(grade_g[i])); + } + } + gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ bool Seg::IsIcefront(void){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17971) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17972) @@ -177,5 +177,5 @@ void GradjDragSSA(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 GradjDragGradient(Vector* gradient,int control_index); void GradjDhDtBalancedthickness(Vector* gradient,int control_index){_error_("not implemented yet");}; void GradjVxBalancedthickness(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 17971) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17972) @@ -2887,6 +2887,4 @@ } break; - GradjDragSSA(gradient,control_index); - break; case MaterialsRheologyBbarEnum: GradjBSSA(gradient,control_index); @@ -2933,5 +2931,24 @@ break; case DragCoefficientAbsGradientEnum: - GradjDragGradient(gradient,control_index); + inputs->GetInputValue(&approximation,ApproximationEnum); + switch(approximation){ + case SSAApproximationEnum: + GradjDragGradient(gradient,control_index); + break; + case FSApproximationEnum:{ + if(IsFloating() || !IsOnBase()) return; + int index1,index2; + this->EdgeOnBaseIndices(&index1,&index2); + Seg* seg = SpawnSeg(index1,index2); + seg->GradjDragGradient(gradient,control_index); + seg->DeleteMaterials(); delete seg; + break; + } + case NoneApproximationEnum: + /*Gradient is 0*/ + break; + default: + _error_("approximation " << EnumToStringx(approximation) << " not supported yet"); + } break; case RheologyBbarAbsGradientEnum: Index: /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp (revision 17971) +++ /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp (revision 17972) @@ -10,16 +10,12 @@ void DragCoefficientAbsGradientx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ - /*Intermediary*/ - int i; - Element* element=NULL; - /*output: */ - IssmDouble J=0; + IssmDouble J=0.; IssmDouble J_sum; /*Compute Misfit: */ - for (i=0;iSize();i++){ - element=dynamic_cast(elements->GetObjectByOffset(i)); - J+=element->DragCoefficientAbsGradient(); + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + J+=DragCoefficientAbsGradient(element); } @@ -32,2 +28,71 @@ *pJ=J; } + +IssmDouble DragCoefficientAbsGradient(Element* element){ + + int domaintype,numcomponents; + bool surfaceintegral; + IssmDouble Jelem=0.; + IssmDouble misfit,Jdet; + IssmDouble dp[2],weight; + IssmDouble* xyz_list = NULL; + + /*Get basal element*/ + if(!element->IsOnBase()) return 0.; + + /*If on water, return 0: */ + if(!element->IsIceInElement()) return 0.; + + /*Get problem dimension*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DverticalEnum: + surfaceintegral = true; + numcomponents = 1; + break; + case Domain3DEnum: + surfaceintegral = true; + numcomponents = 2; + break; + case Domain2DhorizontalEnum: + surfaceintegral = false; + numcomponents = 2; + break; + default: _error_("not supported yet"); + } + + /*Spawn basal element*/ + Element* basalelement = element->SpawnBasalElement(); + + /* Get node coordinates*/ + basalelement->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + Input* weights_input=basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* drag_input =basalelement->GetInput(FrictionCoefficientEnum); _assert_(drag_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum); + drag_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); + + /*Compute Tikhonov regularization J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ + Jelem+=weight*.5*dp[0]*dp[0]*Jdet*gauss->weight; + if(numcomponents==2) Jelem+=weight*.5*dp[1]*dp[1]*Jdet*gauss->weight; + + } + + /*clean up and Return: */ + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + xDelete(xyz_list); + delete gauss; + return Jelem; +} Index: /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h (revision 17971) +++ /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h (revision 17972) @@ -10,4 +10,5 @@ /* local prototypes: */ void DragCoefficientAbsGradientx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); +IssmDouble DragCoefficientAbsGradient(Element* element); #endif Index: /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp (revision 17971) +++ /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp (revision 17972) @@ -10,15 +10,11 @@ void SurfaceAbsVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ - /*Intermediary*/ - int i; - Element* element=NULL; - /*output: */ - IssmDouble J=0; + IssmDouble J=0.; IssmDouble J_sum; /*Compute Misfit: */ - for (i=0;iSize();i++){ - element=dynamic_cast(elements->GetObjectByOffset(i)); + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); J+=SurfaceAbsVelMisfit(element); } @@ -40,5 +36,5 @@ IssmDouble misfit,Jdet; IssmDouble vx,vy,vxobs,vyobs,weight; - IssmDouble* xyz_list_top = NULL; + IssmDouble* xyz_list = NULL; /*Get basal element*/ @@ -66,23 +62,23 @@ } + /*Spawn surface element*/ + Element* topelement = element->SpawnTopElement(); + /* Get node coordinates*/ - if(surfaceintegral) element->GetVerticesCoordinatesTop(&xyz_list_top); - else element->GetVerticesCoordinates(&xyz_list_top); + topelement->GetVerticesCoordinates(&xyz_list); /*Retrieve all inputs we will be needing: */ - Input* weights_input=element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); - Input* vxobs_input =element->GetInput(InversionVxObsEnum); _assert_(vxobs_input); + Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* vx_input =topelement->GetInput(VxEnum); _assert_(vx_input); + Input* vxobs_input =topelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input); Input* vy_input = NULL; Input* vyobs_input = NULL; if(numcomponents==2){ - vy_input =element->GetInput(VyEnum); _assert_(vy_input); - vyobs_input =element->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + vy_input =topelement->GetInput(VyEnum); _assert_(vy_input); + vyobs_input =topelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input); } /* Start looping on the number of gaussian points: */ - Gauss* gauss=NULL; - if(surfaceintegral) gauss=element->NewGaussTop(2); - else gauss=element->NewGauss(2); + Gauss* gauss=topelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ @@ -90,6 +86,5 @@ /* Get Jacobian determinant: */ - if(surfaceintegral) element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss); - else element->JacobianDeterminant(&Jdet,xyz_list_top,gauss); + topelement->JacobianDeterminant(&Jdet,xyz_list,gauss); /*Get all parameters at gaussian point*/ @@ -117,4 +112,6 @@ /*clean up and Return: */ + if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;}; + xDelete(xyz_list); delete gauss; return Jelem;