Changeset 205

Timestamp:

05/04/09 10:55:32 (16 years ago)

Author:

Mathieu Morlighem

Message:

Linearize Du before the integration

File:

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

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

@@ -663 +663 @@
                 GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
 
-                if (velocity_param){
-                        DL_scalar=alpha2*gauss_weight*Jdet;
-                }
-                else{
-                        DL_scalar=0;
-                }
-                        
+                DL_scalar=alpha2*gauss_weight*Jdet;
                 for (i=0;i<2;i++){
                         DL[i][i]=DL_scalar;
@@ -922 +916 @@
         B_param=ParameterInputsRecover(inputs,"B");
         if(temperature_average_param){
+                temperature_average-0;
                 for(i=0;i<3;i++) temperature_average+=temperature_average_param[doflist[i*numberofdofspernode+0]];
                 temperature_average=temperature_average/3.0;
@@ -1374 +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];
 
         /* gaussian points: */
@@ -1399 +1400 @@
 
         /*relative and algorithmic fitting: */
-        double meanvel=0;
-        double epsvel=0;
-        double scalex=1;
-        double scaley=1;
+        double scalex=0;
+        double scaley=0;
         double scale=0;
         double* fit_param=NULL;
@@ -1427 +1426 @@
                 obs_vx_list[i]=obs_velocity[doflist[i*numberofdofspernode+0]];
                 obs_vy_list[i]=obs_velocity[doflist[i*numberofdofspernode+1]];
-
-        }
-
+        }
+
+        /*Get Du at the 3 nodes (integration of the linearized function)*/
+        if(fit==0){
+                /*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];
+                }
+        }
+        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;
+                        relativex[i]=scalex*(vx_list[i]-obs_vx_list[i]);
+                        relativey[i]=scaley*(vy_list[i]-obs_vy_list[i]);
+                }
+        }
+        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.
+                        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];
+                }
+        }
+        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): */
@@ -1479 +1510 @@
                         /*We are using an absolute misfit: */
                         for (i=0;i<numgrids;i++){
-                                due_g_gaussian[i*NDOF2+0]=(velocity_x-obs_velocity_x)*Jdet*gauss_weight*l1l2l3[i]; 
-                                due_g_gaussian[i*NDOF2+1]=(velocity_y-obs_velocity_y)*Jdet*gauss_weight*l1l2l3[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]; 
                         }
                 }
@@ -1490 +1521 @@
                         if(obs_velocity_y==0)scaley=0;
                         for (i=0;i<numgrids;i++){
-                                due_g_gaussian[i*NDOF2+0]=scalex*(velocity_x-obs_velocity_x)*Jdet*gauss_weight*l1l2l3[i]; 
-                                due_g_gaussian[i*NDOF2+1]=scaley*(velocity_y-obs_velocity_y)*Jdet*gauss_weight*l1l2l3[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]; 
                         }
                 }
                 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.
@@ -1501 +1531 @@
                         scale=8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
                         for (i=0;i<numgrids;i++){
-                                due_g_gaussian[i*NDOF2+0]=scale*velocity_x*Jdet*gauss_weight*l1l2l3[i]; 
-                                due_g_gaussian[i*NDOF2+1]=scale*velocity_y*Jdet*gauss_weight*l1l2l3[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]; 
                         }
                 }
@@ -1510 +1540 @@
                 }
                 
-                
                 /*Add due_g_gaussian vector to due_g: */
                 for( i=0; i<numdofs; i++){
@@ -1516 +1545 @@
                 }
         }
-        
         
         /*Add due_g to global vector du_g: */
@@ -2127 +2155 @@
                                         &Ke_gg_gaussian[0][0],0);
 
+                
                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
                 for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
@@ -2356 +2385 @@
                 /*Get bed slope: */
                 GetParameterDerivativeValue(&slope[0], &b[0],&xyz_list[0][0], gauss_l1l2l3);
-                dbdx=slope[0];
-                dbdy=slope[1];
+                dbdx=-slope[0];
+                dbdy=-slope[1];
 
                 /* Get Jacobian determinant: */