Changeset 232

Timestamp:

05/05/09 08:51:00 (16 years ago)

Author:

Mathieu Morlighem

Message:

Fixed Du

File:

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

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

@@ -1369 +1369 @@
         double obs_vx_list[numgrids];
         double obs_vy_list[numgrids];
-        double absolutex[numgrids];
-        double absolutey[numgrids];
-        double relativex[numgrids];
-        double relativey[numgrids];
-        double logarithmicx[numgrids];
-        double logarithmicy[numgrids];
+        double absolutex_list[numgrids];
+        double absolutey_list[numgrids];
+        double relativex_list[numgrids];
+        double relativey_list[numgrids];
+        double logarithmicx_list[numgrids];
+        double logarithmicy_list[numgrids];
 
         /* gaussian points: */
@@ -1388 +1388 @@
         double  velocity_x,velocity_y,velocity_mag;
         double  obs_velocity_x,obs_velocity_y,obs_velocity_mag;
+        double  absolutex,absolutey,relativex,relativey,logarithmicx,logarithmicy;
 
         /*element vector : */
@@ -1432 +1433 @@
                 /*We are using an absolute misfit: */
                 for (i=0;i<numgrids;i++){
-                        absolutex[i]=vx_list[i]-obs_vx_list[i];
-                        absolutey[i]=vy_list[i]-obs_vy_list[i];
+                        absolutex_list[i]=vx_list[i]-obs_vx_list[i];
+                        absolutey_list[i]=vy_list[i]-obs_vy_list[i];
                 }
         }
@@ -1443 +1444 @@
                         if(obs_vx_list[i]==0)scalex=0;
                         if(obs_vy_list[i]==0)scaley=0;
-                        relativex[i]=scalex*(vx_list[i]-obs_vx_list[i]);
-                        relativey[i]=scaley*(vy_list[i]-obs_vy_list[i]);
+                        relativex_list[i]=scalex*(vx_list[i]-obs_vx_list[i]);
+                        relativey_list[i]=scaley*(vy_list[i]-obs_vy_list[i]);
                 }
         }
@@ -1453 +1454 @@
                         obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
                         scale=8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
-                        logarithmicx[i]=scale*vx_list[i];
-                        logarithmicy[i]=scale*vy_list[i];
+                        logarithmicx_list[i]=scale*vx_list[i];
+                        logarithmicy_list[i]=scale*vy_list[i];
                 }
         }
@@ -1509 +1510 @@
                 if(fit==0){
                         /*We are using an absolute misfit: */
+
+                        /*Compute absolute(x/y) at gaussian point: */
+                        GetParameterValue(&absolutex, &absolutex_list[0],gauss_l1l2l3);
+                        GetParameterValue(&absolutey, &absolutey_list[0],gauss_l1l2l3);
+
+                        /*compute Du*/
                         for (i=0;i<numgrids;i++){
-                                due_g_gaussian[i*NDOF2+0]=absolutex[i]*Jdet*gauss_weight*l1l2l3[i]; 
-                                due_g_gaussian[i*NDOF2+1]=absolutey[i]*Jdet*gauss_weight*l1l2l3[i]; 
+                                due_g_gaussian[i*NDOF2+0]=absolutex*Jdet*gauss_weight*l1l2l3[i]; 
+                                due_g_gaussian[i*NDOF2+1]=absolutey*Jdet*gauss_weight*l1l2l3[i]; 
                         }
                 }
                 else if(fit==1){
                         /*We are using a 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;
+
+                        /*Compute relative(x/y) at gaussian point: */
+                        GetParameterValue(&relativex, &relativex_list[0],gauss_l1l2l3);
+                        GetParameterValue(&relativey, &relativey_list[0],gauss_l1l2l3);
+
+                        /*compute Du*/
                         for (i=0;i<numgrids;i++){
-                                due_g_gaussian[i*NDOF2+0]=relativex[i]*Jdet*gauss_weight*l1l2l3[i]; 
-                                due_g_gaussian[i*NDOF2+1]=relativey[i]*Jdet*gauss_weight*l1l2l3[i]; 
+                                due_g_gaussian[i*NDOF2+0]=relativex*Jdet*gauss_weight*l1l2l3[i]; 
+                                due_g_gaussian[i*NDOF2+1]=relativey*Jdet*gauss_weight*l1l2l3[i]; 
                         }
                 }
                 else if(fit==2){
                         /*We are using a 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.
-                        scale=8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
+
+                        /*Compute logarithmic(x/y) at gaussian point: */
+                        GetParameterValue(&logarithmicx, &logarithmicx_list[0],gauss_l1l2l3);
+                        GetParameterValue(&logarithmicy, &logarithmicy_list[0],gauss_l1l2l3);
+
+                        /*compute Du*/
                         for (i=0;i<numgrids;i++){
-                                due_g_gaussian[i*NDOF2+0]=logarithmicx[i]*Jdet*gauss_weight*l1l2l3[i]; 
-                                due_g_gaussian[i*NDOF2+1]=logarithmicy[i]*Jdet*gauss_weight*l1l2l3[i]; 
+                                due_g_gaussian[i*NDOF2+0]=logarithmicx*Jdet*gauss_weight*l1l2l3[i]; 
+                                due_g_gaussian[i*NDOF2+1]=logarithmicy*Jdet*gauss_weight*l1l2l3[i]; 
                         }
                 }
@@ -1816 +1828 @@
 
         /*element vector at the gaussian points: */
-        double  grade_g[numgrids];
+        double  grade_g[numdofs];
         double  grade_g_gaussian[numgrids];
 
@@ -1912 +1924 @@
 
                 /*Add grade_g_gaussian to grade_g: */
-                for( i=0; i<numdofs;i++)grade_g[i]+=grade_g_gaussian[i];
+                for( i=0; i<numgrids;i++) grade_g[i*numberofdofspernode]+=grade_g_gaussian[i];
         }