Changeset 8651

Timestamp:

06/16/11 15:51:23 (14 years ago)

Author:

Mathieu Morlighem

Message:

Added regularization for balance thickness

File:

• issm/trunk/src/c/objects/Elements/Tria.cpp (modified) (5 diffs)

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

@@ -1634 +1634 @@
                 thickness_input->GetParameterDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
                 thicknessobs_input->GetParameterValue(&thicknessobs, gauss);
-                weights_input->GetParameterValue(&weight, gauss,0);
 
                 /*Loop over all requested responses*/
@@ -1640 +1639 @@
 
                         case ThicknessAbsMisfitEnum:
+                                weights_input->GetParameterValue(&weight, gauss,resp);
                                 for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*l1l2l3[i];
-                                /*Regularization of the constraint: 2000000 79 N*/
-                                //for(i=0;i<numdof;i++) pe->values[i]+= - 1*100000*dH[0]*dbasis[0][i]*Jdet*gauss->weight;
-                                //for(i=0;i<numdof;i++) pe->values[i]+= - 1*100000*dH[1]*dbasis[1][i]*Jdet*gauss->weight;
+                                break;
+                        case ThicknessAbsGradientEnum:
+                                weights_input->GetParameterValue(&weight, gauss,resp);
+                                for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[0]*dbasis[0][i]*Jdet*gauss->weight;
+                                for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[1]*dbasis[1][i]*Jdet*gauss->weight;
                                 break;
                         default:
@@ -2918 +2920 @@
 
         for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
+                //FIXME: the control type should be checked somewhere (with respect to what variable are we taking the gradient!)
 
                 case ThicknessAbsMisfitEnum:
+                case ThicknessAbsGradientEnum:
                 case SurfaceAbsVelMisfitEnum:
                 case SurfaceRelVelMisfitEnum:
@@ -5349 +5353 @@
 
                 /*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;
+                Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
         }
 
@@ -5397 +5401 @@
                 /*compute ThicknessAbsMisfit*/
                 Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight;
-                //Jelem+=0.5*100000*(pow(dH[0],2.)+pow(dH[1],2.))*Jdet*gauss->weight;
         }