# [cig-commits] r14530 - in long/3D/Gale/trunk: . src/StGermain/Base/Foundation/src src/Underworld/plugins/EulerDeform

walter at geodynamics.org walter at geodynamics.org
Sun Mar 29 14:35:32 PDT 2009

```Author: walter
Date: 2009-03-29 14:35:32 -0700 (Sun, 29 Mar 2009)
New Revision: 14530

Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/StGermain/Base/Foundation/src/Numerics.c
long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c
Log:
r2619 at dante:  boo | 2009-03-29 14:34:56 -0700
Fixe some round-off bugs

Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2617
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2619

Modified: long/3D/Gale/trunk/src/StGermain/Base/Foundation/src/Numerics.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Base/Foundation/src/Numerics.c	2009-03-29 20:38:19 UTC (rev 14529)
+++ long/3D/Gale/trunk/src/StGermain/Base/Foundation/src/Numerics.c	2009-03-29 21:35:32 UTC (rev 14530)
@@ -35,13 +35,15 @@
#include "Numerics.h"

-const double	Num_Epsilon = 1e-12;
+const double    Num_Epsilon = 1e-10;

Bool Num_Approx( double var, double val ) {
-	return (var >= val - Num_Epsilon && var <= val + Num_Epsilon);
+	return (var >= val*(1.0-Num_Epsilon) && var <= val *(1.0+Num_Epsilon))
+          || (var >= val - Num_Epsilon && var <= val + Num_Epsilon);;
}

Bool Num_InRange( double var, double low, double upp ) {
-	return (var >= low - Num_Epsilon && var <= upp + Num_Epsilon);
+	return (var >= low *(1.0-Num_Epsilon) && var <= upp *(1.0+Num_Epsilon))
+          || (var >= low - Num_Epsilon && var <= upp + Num_Epsilon);
}

Modified: long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c	2009-03-29 20:38:19 UTC (rev 14529)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c	2009-03-29 21:35:32 UTC (rev 14530)
@@ -991,7 +991,7 @@
/* Interpolate the value. */
res = FieldVariable_InterpolateValueAt( field, newCrds[n_i], newVals + n_i * field->fieldComponentCount );
if( res == OTHER_PROC || res == OUTSIDE_GLOBAL ) {
-                  printf("Coordinate out of bounds %d %lf %lf %lf\n",
+                  printf("Coordinate out of bounds %d %g %g %g\n",
n_i,newCrds[n_i][0],newCrds[n_i][1],newCrds[n_i][2]);

assert( 0 );
@@ -1240,7 +1240,7 @@
#else
_EulerDeform_LineInterp( (const double**)crds, newCrd, 0, 1, &mesh->verts[centerInd][1] );
#endif
-                          mesh->verts[centerInd][1] -= 1e-15;
+                          mesh->verts[centerInd][1] *= (1.0-1.0e-10);
}
}
else if( grm->nDims == 3 ) {
@@ -1288,7 +1288,7 @@
memcpy( crds[3], oldCrds[ind], nDims * sizeof(double) );

if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->verts[centerInd][1] ) ) {
-					mesh->verts[centerInd][1] -= 1e-13;
+					mesh->verts[centerInd][1] *= (1.0-1.0e-10);
return;
}
}
@@ -1312,7 +1312,7 @@

ijk[2]--;
if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->verts[centerInd][1] ) ) {
-					mesh->verts[centerInd][1] -= 1e-13;
+					mesh->verts[centerInd][1] *= (1.0-1.0e-10);
return;
}
}
@@ -1336,7 +1336,7 @@

ijk[0]--;
if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->verts[centerInd][1] ) ) {
-					mesh->verts[centerInd][1] -= 1e-13;
+					mesh->verts[centerInd][1] *= (1.0-1.0e-10);
return;
}
}
@@ -1360,7 +1360,7 @@

ijk[0]--; ijk[2]--;
if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->verts[centerInd][1] ) ) {
-					mesh->verts[centerInd][1] -= 1e-13;
+					mesh->verts[centerInd][1] *= (1.0-1.0e-10);
return;
}
}
@@ -1422,7 +1422,7 @@
/* Calc barycenter. */
a1 = (newCrd[0] - leftCrd[0]) / (oldCrd[0] - leftCrd[0]);
a0 = 1.0 - a1;
-				mesh->nodeCoord[centerInd][1] = a0 * leftCrd[1] + a1 * oldCrd[1] + 1e-15;
+				mesh->nodeCoord[centerInd][1] = (a0 * leftCrd[1] + a1 * oldCrd[1])*(1.0+ 1e-10);
}
else {
XYZ		rightCrd;
@@ -1439,7 +1439,7 @@
/* Calc barycenter. */
a1 = (newCrd[0] - oldCrd[0]) / (rightCrd[0] - oldCrd[0]);
a0 = 1.0 - a1;
-				mesh->nodeCoord[centerInd][1] = a0 * oldCrd[1] + a1 * rightCrd[1] + 1e-15;
+				mesh->nodeCoord[centerInd][1] = (a0 * oldCrd[1] + a1 * rightCrd[1])*(1.0 + 1e-10);
}
}
else if( grm->nDims == 3 ) {
@@ -1473,7 +1473,7 @@
ijk[0]--; ijk[2]++;	GRM_Project( grm, ijk, &ind ); memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
ijk[0]++;		GRM_Project( grm, ijk, &ind ); memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
-					mesh->nodeCoord[centerInd][1] += 1e-15;
+                                  mesh->nodeCoord[centerInd][1] *=(1.0+1e-10);
return;
}
}
@@ -1485,7 +1485,7 @@
ijk[0]++;		GRM_Project( grm, ijk, &ind ); memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
ijk[2]--;
if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
-					mesh->nodeCoord[centerInd][1] += 1e-15;
+					mesh->nodeCoord[centerInd][1] *=(1.0+1e-10);
return;
}
}
@@ -1497,7 +1497,7 @@
ijk[0]++;		GRM_Project( grm, ijk, &ind ); memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
ijk[0]--;
if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
-					mesh->nodeCoord[centerInd][1] += 1e-15;
+					mesh->nodeCoord[centerInd][1] *=(1.0+1e-10);
return;
}
}
@@ -1509,7 +1509,7 @@
ijk[0]++;		GRM_Project( grm, ijk, &ind ); memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
ijk[0]--; ijk[2]--;
if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
-					mesh->nodeCoord[centerInd][1] += 1e-15;
+					mesh->nodeCoord[centerInd][1] *=(1.0+1e-10);
return;
}
}
@@ -1564,7 +1564,7 @@
/* Calc barycenter. */
a1 = (newCrd[1] - leftCrd[1]) / (oldCrd[1] - leftCrd[1]);
a0 = 1.0 - a1;
-				mesh->nodeCoord[centerInd][0] = a0 * leftCrd[0] + a1 * oldCrd[0] + 1e-15;
+				mesh->nodeCoord[centerInd][0] = (a0 * leftCrd[0] + a1 * oldCrd[0])*(1.0+1.0e-10);
}
else {
XYZ		rightCrd;
@@ -1581,7 +1581,7 @@
/* Calc barycenter. */
a1 = (newCrd[1] - oldCrd[1]) / (rightCrd[1] - oldCrd[1]);
a0 = 1.0 - a1;
-				mesh->nodeCoord[centerInd][0] = a0 * oldCrd[0] + a1 * rightCrd[0] + 1e-15;
+				mesh->nodeCoord[centerInd][0] = (a0 * oldCrd[0] + a1 * rightCrd[0])*(1.0+1.0e-10);
}
}
else if( grm->nDims == 3 ) {

```

More information about the CIG-COMMITS mailing list