Changeset 17972

Timestamp:

05/09/14 14:03:34 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added DragCoefficientAbsGradientx in module

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Seg.cpp (modified) (1 diff)
• classes/Elements/Seg.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (2 diffs)
• modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp (modified) (2 diffs)
• modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h (modified) (1 diff)
• modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/classes/Elements/Seg.cpp ¶

@@ -195 +195 @@
 }
 /*}}}*/
+void       Seg::GradjDragGradient(Vector<IssmDouble>* 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();ig<gauss->end();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;i<NUMVERTICES;i++){
+                        grade_g[i]+=-weight*Jdet*gauss->weight*dbasis[0][i]*dk[0];
+                        _assert_(!xIsNan<IssmDouble>(grade_g[i]));
+                }
+        }
+        gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
+
+        /*Clean up and return*/
+        delete gauss;
+}
+/*}}}*/
 bool       Seg::IsIcefront(void){/*{{{*/
 

issm/trunk-jpl/src/c/classes/Elements/Seg.h ¶

@@ -177 +177 @@
                 void       GradjDragSSA(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
                 void       GradjDragFS(Vector<IssmDouble>* gradient,int control_index);
-                void       GradjDragGradient(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
+                void       GradjDragGradient(Vector<IssmDouble>* gradient,int control_index);
                 void       GradjDhDtBalancedthickness(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
                 void       GradjVxBalancedthickness(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp ¶

@@ -2887 +2887 @@
                         }
                         break;
-                        GradjDragSSA(gradient,control_index);
-                        break;
                 case MaterialsRheologyBbarEnum:
                         GradjBSSA(gradient,control_index);
@@ -2933 +2931 @@
                         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:

issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp ¶

@@ -10 +10 @@
 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;i<elements->Size();i++){
-                element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-                J+=element->DragCoefficientAbsGradient();
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+                J+=DragCoefficientAbsGradient(element);
         }
 
@@ -32 +28 @@
         *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();ig<gauss->end();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<IssmDouble>(xyz_list);
+        delete gauss;
+        return Jelem;
+}

issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h ¶

@@ -10 +10 @@
 /* local prototypes: */
 void DragCoefficientAbsGradientx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters);
+IssmDouble DragCoefficientAbsGradient(Element* element);
 
 #endif

issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp ¶

@@ -10 +10 @@
 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;i<elements->Size();i++){
-                element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
                 J+=SurfaceAbsVelMisfit(element);
         }
@@ -40 +36 @@
         IssmDouble misfit,Jdet;
         IssmDouble vx,vy,vxobs,vyobs,weight;
-        IssmDouble* xyz_list_top = NULL;
+        IssmDouble* xyz_list = NULL;
 
         /*Get basal element*/
@@ -66 +62 @@
         }
 
+        /*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();ig<gauss->end();ig++){
 
@@ -90 +86 @@
 
                 /* 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 +112 @@
 
         /*clean up and Return: */
+        if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
+        xDelete<IssmDouble>(xyz_list);
         delete gauss;
         return Jelem;