Changeset 13528

Timestamp:

10/03/12 20:58:05 (12 years ago)

Author:

Mathieu Morlighem

Message:

BUG: reverted back to original regularization

File:

• issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (modified) (3 diffs)

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

@@ -3659 +3659 @@
                         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.));
-                        vz=0;
-                        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);
-                        K[2][0]=h/(2*vel)*(vz*vx);  K[2][1]=h/(2*vel)*(vz*vy); K[2][2]=h/(2*vel)*(vz*vz);
+
+                        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);
+                        K[1][0]=h/(2*vel)*fabs(vy*vx);  K[1][1]=h/(2*vel)*fabs(vy*vy); K[1][2]=h/(2*vel)*fabs(vy*vz);
+                        K[2][0]=h/(2*vel)*fabs(vz*vx);  K[2][1]=h/(2*vel)*fabs(vz*vy); K[2][2]=h/(2*vel)*fabs(vz*vz);
 
                         D_scalar_stab=gauss->weight*Jdet;
@@ -4319 +4319 @@
         IssmDouble B_average,s_average;
         int*   doflist=NULL;
-        IssmDouble pressure[numdof];
+        //IssmDouble pressure[numdof];
 
         /*Get dof list: */
@@ -4327 +4327 @@
         for(i=0;i<numdof;i++){
                 values[i]=solution[doflist[i]];
-                GetInputListOnVertices(&pressure[0],PressureEnum);
-                if(values[i]>matpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]);
-                if(values[i]<matpar->TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50;
+                //GetInputListOnVertices(&pressure[0],PressureEnum);
+                //if(values[i]>matpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]);
+                //if(values[i]<matpar->TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50;
 
                 /*Check solution*/