OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Updated common/pmmc.h

Browse Source
This commit is contained in:
James E McClure 2015-11-14 13:05:37 -05:00
parent 6ad044fd36
commit e31ac627e2
1 changed files with 63 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
common/pmmc.h
View File

@ -4355,8 +4355,10 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
if (norm==0.0) norm=1.0;
// Compute the interface speed from time derivative and gradient (Level Set Equation)
zeta = -Pt.eval(P) / norm;
//temp = sqrt(temp)/norm; <--- what was I thinking with this? (James)
// Normalize the normal vector
x /= norm;
y /= norm;
z /= norm;
// Compute the average
AvgVel(0) += area*zeta*x;
@ -4420,6 +4422,65 @@ inline double geomavg_EulerCharacteristic(DTMutableList<Point> &Points, IntArray
if (!newside) nside-=1;
}
EulerChar = 1.0*(nvert - nside + nface);
return EulerChar;
}
inline double mink_EulerCharacteristic(DTMutableList<Point> &Points, IntArray &Triangles,
DoubleArray &CubeValues, int &npts, int &ntris, int &i, int &j, int &k){
// Compute the Euler characteristic for triangles in a cube
// Exclude edges and vertices shared with between multiple cubes
double EulerChar;
int nvert=npts;
int nside=2*nvert-3;
int nface=nvert-2;
//if (ntris != nface){
// nface = ntris;
// nside =
//}
//...........................................................
// Check that this point is not on a previously computed face
// Note direction that the marching cubes algorithm marches
// In parallel, other sub-domains fill in the lower boundary
for (int p=0; p<npts; p++){
Point PT = Points(p);
if (PT.x - double(i) < 1e-12) nvert-=1;
else if (PT.y - double(j) < 1e-12) nvert-=1;
else if (PT.z - double(k) < 1e-12) nvert-=1;
}
// Remove redundantly computed edges (shared by two cubes across a cube face)
for (int p=0; p<ntris; p++){
Point A = Points(Triangles(0,p));
Point B = Points(Triangles(1,p));
Point C = Points(Triangles(2,p));
// Check side A-B
bool newside = true;
if (A.x - double(i) < 1e-12 && B.x - double(i) < 1e-12) newside=false;
if (A.y - double(j) < 1e-12 && B.y - double(j) < 1e-12) newside=false;
if (A.z - double(k) < 1e-12 && B.z - double(k) < 1e-12) newside=false;
if (!newside) nside-=1;
// Check side A-C
newside = true;
if (A.x - double(i)< 1e-12 && C.x - double(i) < 1e-12) newside=false;
if (A.y - double(j)< 1e-12 && C.y - double(j) < 1e-12) newside=false;
if (A.z - double(k)< 1e-12 && C.z - double(k) < 1e-12) newside=false;
if (!newside) nside-=1;
// Check side B-C
newside = true;
if (B.x - double(i) < 1e-12 && C.x - double(i) < 1e-12) newside=false;
if (B.y - double(j) < 1e-12 && C.y - double(j) < 1e-12) newside=false;
if (B.z - double(k) < 1e-12 && C.z - double(k) < 1e-12) newside=false;
if (!newside) nside-=1;
}
EulerChar = 1.0*(nvert - nside + nface);
return EulerChar;
}