Changeset 8611

Timestamp:

06/10/11 17:57:23 (14 years ago)

Author:

Mathieu Morlighem

Message:

preparing thikonov regularization in misfits

Location:

issm/trunk/src

Files:

• c/objects/Elements/Penta.cpp (modified) (2 diffs)
• c/objects/Elements/Tria.cpp (modified) (5 diffs)
• m/model/ismodelselfconsistent.m (modified) (1 diff)

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -3849 +3849 @@
 double Penta::DragCoefficientAbsGradient(bool process_units,int weight_index){
 
-        _error_("Not implemented yet");
+        double J;
+        Tria*  tria=NULL;
+
+        /*If on water, on shelf or not on bed, skip: */
+        if(IsOnWater()|| IsOnShelf() || !IsOnBed()) return 0;
+
+        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria
+        J=tria->DragCoefficientAbsGradient(process_units,weight_index);
+        delete tria->matice; delete tria;
+        return J;
 }
 /*}}}*/
@@ -6751 +6760 @@
 double Penta::RheologyBbarAbsGradient(bool process_units,int weight_index){
 
-        _error_("Not implemented yet");
+        double J;
+        Tria*  tria=NULL;
+
+        /*If on water, on shelf or not on bed, skip: */
+        if(IsOnWater() || !IsOnBed()) return 0;
+
+        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria
+        J=tria->RheologyBbarAbsGradient(process_units,weight_index);
+        delete tria->matice; delete tria;
+        return J;
 }
 /*}}}*/

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -1812 +1812 @@
                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
                                         }
+                                        break;
+                                case DragCoefficientAbsGradientEnum:
+                                        /*Nothing in P vector*/
+                                        break;
+                                case ThicknessAbsGradientEnum:
+                                        /*Nothing in P vector*/
+                                        break;
+                                case RheologyBbarAbsGradientEnum:
+                                        /*Nothing in P vector*/
                                         break;
                                 default:
@@ -2396 +2405 @@
 double Tria::DragCoefficientAbsGradient(bool process_units,int weight_index){
 
-        _error_("Not implemented yet");
+        /* Intermediaries */
+        int        ig;
+        double     Jelem = 0;
+        double     weight;
+        double     Jdet;
+        double     xyz_list[NUMVERTICES][3];
+        double     dp[NDOF2];
+        GaussTria *gauss = NULL;
+
+        /*retrieve parameters and inputs*/
+
+        /*If on water, return 0: */
+        if(IsOnWater()) return 0;
+
+        /*Retrieve all inputs we will be needing: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        Input* weights_input=inputs->GetInput(WeightsEnum);         _assert_(weights_input);
+        Input* drag_input   =inputs->GetInput(DragCoefficientEnum); _assert_(drag_input);
+
+        /* Start looping on the number of gaussian points: */
+        gauss=new GaussTria(2);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+
+                /*Get all parameters at gaussian point*/
+                weights_input->GetParameterValue(&weight,gauss,weight_index);
+                drag_input->GetParameterDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
+
+                /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+                //Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
+        }
+
+        /*Clean up and return*/
+        delete gauss;
+        return Jelem;
 }
 /*}}}*/
@@ -4679 +4726 @@
 double Tria::RheologyBbarAbsGradient(bool process_units,int weight_index){
 
-        _error_("Not implemented yet");
+        /* Intermediaries */
+        int        ig;
+        double     Jelem = 0;
+        double     weight;
+        double     Jdet;
+        double     xyz_list[NUMVERTICES][3];
+        double     dp[NDOF2];
+        GaussTria *gauss = NULL;
+
+        /*retrieve parameters and inputs*/
+
+        /*If on water, return 0: */
+        if(IsOnWater()) return 0;
+
+        /*Retrieve all inputs we will be needing: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        Input* weights_input  =inputs->GetInput(WeightsEnum);              _assert_(weights_input);
+        Input* rheologyb_input=matice->inputs->GetInput(RheologyBbarEnum); _assert_(rheologyb_input);
+
+        /* Start looping on the number of gaussian points: */
+        gauss=new GaussTria(2);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+
+                /*Get all parameters at gaussian point*/
+                weights_input->GetParameterValue(&weight,gauss,weight_index);
+                rheologyb_input->GetParameterDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
+
+                /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+                //Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
+        }
+
+        /*Clean up and return*/
+        delete gauss;
+        return Jelem;
 }
 /*}}}*/
@@ -5211 +5296 @@
 double Tria::ThicknessAbsGradient(bool process_units,int weight_index){
 
-        _error_("Not implemented yet");
+        /* Intermediaries */
+        int        ig;
+        double     Jelem = 0;
+        double     weight;
+        double     Jdet;
+        double     xyz_list[NUMVERTICES][3];
+        double     dp[NDOF2];
+        GaussTria *gauss = NULL;
+
+        /*retrieve parameters and inputs*/
+
+        /*If on water, return 0: */
+        if(IsOnWater()) return 0;
+
+        /*Retrieve all inputs we will be needing: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        Input* weights_input  =inputs->GetInput(WeightsEnum);   _assert_(weights_input);
+        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+
+        /* Start looping on the number of gaussian points: */
+        gauss=new GaussTria(2);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+
+                /*Get all parameters at gaussian point*/
+                weights_input->GetParameterValue(&weight,gauss,weight_index);
+                thickness_input->GetParameterDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
+
+                /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+                //Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
+        }
+
+        /*Clean up and return*/
+        delete gauss;
+        return Jelem;
 }
 /*}}}*/
 /*FUNCTION Tria::ThicknessAbsMisfit {{{1*/
 double Tria::ThicknessAbsMisfit(bool process_units,int weight_index){
-
-        /* Constants */
-        const int    numdof=1*NUMVERTICES;
 
         /*Intermediaries*/
@@ -5233 +5353 @@
         if(IsOnWater())return 0;
 
-        /* Get node coordinates and dof list: */
+        /*Retrieve all inputs we will be needing: */
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-
-        /*Retrieve all inputs we will be needing: */
         Input* thickness_input   =inputs->GetInput(ThicknessEnum);   _assert_(thickness_input);
         Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);_assert_(thicknessobs_input);

issm/trunk/src/m/model/ismodelselfconsistent.m

@@ -199 +199 @@
         checkvalues(md,{'cm_responses'},...
                 [SurfaceAbsVelMisfitEnum SurfaceRelVelMisfitEnum SurfaceLogVelMisfitEnum SurfaceLogVxVyMisfitEnum SurfaceAverageVelMisfitEnum...
-                ThicknessAbsMisfitEnum]);
+                ThicknessAbsMisfitEnum DragCoefficientAbsGradientEnum RheologyBbarAbsGradientEnum]);
 
         %WEIGHTS