Changeset 15741 for issm/trunk-jpl/src/c/classes/Elements

Timestamp:

08/07/13 11:25:05 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: keep the exponent as an integer may execute faster (compiler dependent)

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Penta.cpp (modified) (3 diffs)
• PentaRef.cpp (modified) (1 diff)
• Tria.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -3938 +3938 @@
                         /*Build K: */
                         GetElementSizes(&hx,&hy,&hz);
-                        vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14;
-                        h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.));
+                        vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
+                        h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
                         K[0][0]=h/(2*vel)*vx*vx;  K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz;
                         K[1][0]=h/(2*vel)*vy*vx;  K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz;
@@ -4167 +4167 @@
                         /*Build K: */
                         GetElementSizes(&hx,&hy,&hz);
-                        vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14;
-                        h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.));
+                        vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
+                        h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
 
                         K[0][0]=h/(2*vel)*fabs(vx*vx);  K[0][1]=h/(2*vel)*fabs(vx*vy); K[0][2]=h/(2*vel)*fabs(vx*vz);
@@ -8496 +8496 @@
                 /*Deal with lower surface*/
                 if(IsOnBed()){
-                        constant_part=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness;
+                        constant_part=-1.58*1.e-10*rho_ice*gravity*thickness;
                         ub=constant_part*slope[0];
                         vb=constant_part*slope[1];

issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp

@@ -1874 +1874 @@
         /*Jdet = (Area of the trapezoid)/(Area trapezoid ref) with AreaRef = 4*/
         /*Area of a trabezoid = altitude * (base1 + base2)/2 */
-        *Jdet= pow(pow(x2-x1,2.) + pow(y2-y1,2.),0.5) * (z4-z1 + z3-z2)/8;
-        if(*Jdet<0) _error_("negative jacobian determinant!");
+        *Jdet= sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)) * (z4-z1 + z3-z2)/8.;
+        if(*Jdet<0.) _error_("negative jacobian determinant!");
 
 }

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -3321 +3321 @@
                 constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.));
 
-                ub=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness*slope[0];
-                vb=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness*slope[1];
+                ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
+                vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
 
                 pe->values[2*i+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/connectivity;
@@ -4429 +4429 @@
 
                 /*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*(dp[0]*dp[0] + dp[1]*dp[1])*Jdet*gauss->weight;
         }