Changeset 419

Timestamp:

05/14/09 08:50:30 (16 years ago)

Author:

Mathieu Morlighem

Message:

better way to compute Misfit integral

File:

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

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

@@ -1577 +1577 @@
 
         /* parameters: */
-        double  velocity_x,velocity_y,velocity_mag;
-        double  obs_velocity_x,obs_velocity_y,obs_velocity_mag;
+        double  obs_velocity_mag,velocity_mag;
         double  absolutex,absolutey,relativex,relativey,logarithmicx,logarithmicy;
 
@@ -2095 +2094 @@
 #define __FUNCT__ "Tria::Misfit"
 double Tria::Misfit(double* velocity,double* obs_velocity,void* vinputs,int analysis_type){
+
         int i;
         
@@ -2113 +2113 @@
         double obs_vx_list[numgrids];
         double obs_vy_list[numgrids];
+        double absolute_list[numgrids];
+        double relative_list[numgrids];
+        double logarithmic_list[numgrids];
 
         /* gaussian points: */
@@ -2124 +2127 @@
 
         /* parameters: */
-        double  velocity_x,velocity_y,velocity_mag;
-        double  obs_velocity_x,obs_velocity_y,obs_velocity_mag;
+        double  velocity_mag,obs_velocity_mag;
+        double  absolute,relative,logarithmic;
 
         /* Jacobian: */
@@ -2155 +2158 @@
         }
         
+        /*Compute Misfit at the 3 nodes (integration of the linearized function)*/
+        if(fit==0){
+                /*We are using an absolute misfit: */
+                for (i=0;i<numgrids;i++){
+                        absolute_list[i]=0.5*(pow((vx_list[i]-obs_vx_list[i]),2)+pow((vy_list[i]-obs_vy_list[i]),2));
+                }
+        }
+        else if(fit==1){
+                /*We are using a relative misfit: */
+                for (i=0;i<numgrids;i++){
+                        scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
+                        scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
+                        if(obs_vx_list[i]==0)scalex=0;
+                        if(obs_vy_list[i]==0)scaley=0;
+                        relative_list[i]=0.5*(scalex*pow((vx_list[i]-obs_vx_list[i]),2)+scaley*pow((vy_list[i]-obs_vy_list[i]),2))*Jdet*gauss_weight;
+                }
+        }
+        else if(fit==2){
+                /*We are using a logarithmic misfit: */
+                for (i=0;i<numgrids;i++){
+                        velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
+                        obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
+                        logarithmic_list[i]=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
+                }
+        }
+        else{
+                /*Not supported yet! : */
+                throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit));
+        }
+
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
@@ -2172 +2205 @@
                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
 
-                /*Compute velocities at gaussian point: */
-                GetParameterValue(&velocity_x, &vx_list[0],gauss_l1l2l3);
-                GetParameterValue(&velocity_y, &vy_list[0],gauss_l1l2l3);
-                #ifdef _DEBUG_ 
-                        printf("Velocity: %g %g\n", velocity_x,velocity_y);
-                #endif
-        
-                /*Compute obs_velocities at gaussian point: */
-                GetParameterValue(&obs_velocity_x, &obs_vx_list[0],gauss_l1l2l3);
-                GetParameterValue(&obs_velocity_y, &obs_vy_list[0],gauss_l1l2l3);
-                #ifdef _DEBUG_ 
-                        printf("Observed velocity: %g %g\n", obs_velocity_x,obs_velocity_y);
-                #endif
-
                 /* Get Jacobian determinant: */
                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
@@ -2194 +2213 @@
                 /*Differents misfits are allowed: */
                 if(fit==0){
-                        /*Absolute misfit: */
-                        Jelem+=.5*(pow((velocity_x-obs_velocity_x),2)+pow((velocity_y-obs_velocity_y),2))*Jdet*gauss_weight;
+                        /*Compute absolute misfit at gaussian point: */
+                        GetParameterValue(&absolute, &absolute_list[0],gauss_l1l2l3);
+
+                        /*compute Misfit*/
+                        Jelem+=absolute*Jdet*gauss_weight;
                 }
                 else if(fit==1){
-                        /*Relative misfit: */
-                        scalex=pow(meanvel/(obs_velocity_x+epsvel),2);
-                        scaley=pow(meanvel/(obs_velocity_y+epsvel),2);
-                        if(obs_velocity_x==0)scalex=0;
-                        if(obs_velocity_y==0)scaley=0;
-                        Jelem+=.5*(scalex*pow((velocity_x-obs_velocity_x),2)+scaley*pow((velocity_y-obs_velocity_y),2))*Jdet*gauss_weight;
+                        /*Compute relative misfit at gaussian point: */
+                        GetParameterValue(&relative, &relative_list[0],gauss_l1l2l3);
+
+                        /*compute Misfit*/
+                        Jelem+=relative*Jdet*gauss_weight;
                 }       
                 else if(fit==2){
-                        /*Logarithmic misfit: */
-                        velocity_mag=sqrt(pow(velocity_x,2)+pow(velocity_y,2))+epsvel; //epsvel to avoid velocity being nil.
-                        obs_velocity_mag=sqrt(pow(obs_velocity_x,2)+pow(obs_velocity_y,2))+epsvel; //epsvel to avoid observed velocity being nil.
-                        Jelem+=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2)*Jdet*gauss_weight;
+                        /*Compute logarithmic misfit at gaussian point: */
+                        GetParameterValue(&logarithmic, &logarithmic_list[0],gauss_l1l2l3);
+
+                        /*compute Misfit*/
+                        Jelem+=logarithmic*Jdet*gauss_weight;
                 }
                 else throw ErrorException(__FUNCT__,exprintf("%s%i%s","fit type",fit," not supported yet!"));