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

Added TriLinearPoly to common/pmmc.h

Browse Source
This commit is contained in:
James E McClure 2014-10-16 10:48:57 -04:00
parent 5b486429d0
commit 157a1bd6f1
1 changed files with 92 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

150
common/pmmc.h
View File

@ -301,7 +301,7 @@ char triTable[256][16] =
{0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}; {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}};
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------
class TriLinearPoly{ class TriLinPoly{
/* Compute a tri-linear polynomial within a given cube (i,j,k) x (i+1,j+1,k+1) /* Compute a tri-linear polynomial within a given cube (i,j,k) x (i+1,j+1,k+1)
* Values are provided at the corners in CubeValues * Values are provided at the corners in CubeValues
* x,y,z must be defined on [0,1] where the length of the cube edge is one * x,y,z must be defined on [0,1] where the length of the cube edge is one
@ -309,34 +309,40 @@ class TriLinearPoly{
int ic,jc,kc; int ic,jc,kc;
double a,b,c,d,e,f,g,h; double a,b,c,d,e,f,g,h;
double x,y,z; double x,y,z;
DoubleArray CubeValues;
public: public:
TriLinearPoly(){ DoubleArray Corners;
CubeValues.New(2,2,2);
TriLinPoly(){
Corners.New(2,2,2);
} }
// Assign the polynomial within a cube
void assign(DoubleArray &A, int i, int j, int k){ void assign(DoubleArray &A, int i, int j, int k){
// Save the cube indices
ic=i; jc=j; kc=k; ic=i; jc=j; kc=k;
CubeValues(0,0,0) = A(i,j,k); // local copy of the cube values
CubeValues(1,0,0) = A(i+1,j,k); Corners(0,0,0) = A(i,j,k);
CubeValues(0,1,0) = A(i,j+1,k); Corners(1,0,0) = A(i+1,j,k);
CubeValues(1,1,0) = A(i+1,j+1,k); Corners(0,1,0) = A(i,j+1,k);
CubeValues(0,0,1) = A(i,j,k+1); Corners(1,1,0) = A(i+1,j+1,k);
CubeValues(1,0,1) = A(i+1,j,k+1); Corners(0,0,1) = A(i,j,k+1);
CubeValues(0,1,1) = A(i,j+1,k+1); Corners(1,0,1) = A(i+1,j,k+1);
CubeValues(1,1,1) = A(i+1,j+1,k+1); Corners(0,1,1) = A(i,j+1,k+1);
Corners(1,1,1) = A(i+1,j+1,k+1);
a = CubeValues(0,0,0); // coefficients of the tri-linear approximation
b = CubeValues(1,0,0)-a; a = Corners(0,0,0);
c = CubeValues(0,1,0)-a; b = Corners(1,0,0)-a;
d = CubeValues(0,0,1)-a; c = Corners(0,1,0)-a;
e = CubeValues(1,1,0)-a-b-c; d = Corners(0,0,1)-a;
f = CubeValues(1,0,1)-a-b-d; e = Corners(1,1,0)-a-b-c;
g = CubeValues(0,1,1)-a-c-d; f = Corners(1,0,1)-a-b-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g; g = Corners(0,1,1)-a-c-d;
h = Corners(1,1,1)-a-b-c-d-e-f-g;
} }
// Evaluate the polynomial at a point
double eval(Point P){ double eval(Point P){
double returnValue; double returnValue;
x = P.x - 1.0*ic; x = P.x - 1.0*ic;
@ -345,6 +351,7 @@ public:
returnValue = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z; returnValue = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
return returnValue; return returnValue;
} }
}; };
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------
inline int ComputeBlob(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indicator, inline int ComputeBlob(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indicator,
@ -4243,50 +4250,76 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, DTMutableList<Po
double sxx,syy,szz,sxy,sxz,syz,sx,sy,sz; double sxx,syy,szz,sxy,sxz,syz,sx,sy,sz;
double Axx,Axy,Axz,Ayx,Ayy,Ayz,Azx,Azy,Azz; double Axx,Axy,Axz,Ayx,Ayy,Ayz,Azx,Azy,Azz;
double Tx[8],Ty[8],Tz[8]; // Tangent vector // double Tx[8],Ty[8],Tz[8]; // Tangent vector
double Nx[8],Ny[8],Nz[8]; // Principle normal // double Nx[8],Ny[8],Nz[8]; // Principle normal
double tx,ty,tz,nx,ny,nz,K; // tangent,normal and curvature double tx,ty,tz,nx,ny,nz,K; // tangent,normal and curvature
/* // Store the tangent and normal locally in the cube
DoubleArray Tx(2,2,2);
DoubleArray Ty(2,2,2);
DoubleArray Tz(2,2,2);
DoubleArray Nx(2,2,2);
DoubleArray Ny(2,2,2);
DoubleArray Nx(2,2,2);
*/
Point P; Point P;
TriLinPoly Tx,Ty,Tz,Nx,Ny,Nz;
for (int p=0; p<8; p++){ // Loop over the cube and compute the derivatives
i = ic; for (k=kc; k<kc+2; k++){
j = jc; for (j=jc; j<jc+2; j++){
k = kc; for (i=ic; i<ic+2; i++){
// Compute all of the derivatives using finite differences
// fluid phase indicator field // Compute all of the derivatives using finite differences
fx = 0.5*(f(i+1,j,k) - f(i-1,j,k)); // fluid phase indicator field
fy = 0.5*(f(i,j+1,k) - f(i,j-1,k)); fx = 0.5*(f(i+1,j,k) - f(i-1,j,k));
fz = 0.5*(f(i,j,k+1) - f(i,j,k-1)); fy = 0.5*(f(i,j+1,k) - f(i,j-1,k));
fxx = f(i+1,j,k) - 2.0*f(i,j,k) + f(i-1,j,k); fz = 0.5*(f(i,j,k+1) - f(i,j,k-1));
fyy = f(i,j+1,k) - 2.0*f(i,j,k) + f(i,j-1,k); fxx = f(i+1,j,k) - 2.0*f(i,j,k) + f(i-1,j,k);
fzz = f(i,j,k+1) - 2.0*f(i,j,k) + f(i,j,k-1); fyy = f(i,j+1,k) - 2.0*f(i,j,k) + f(i,j-1,k);
fxy = 0.25*(f(i+1,j+1,k) - f(i+1,j-1,k) - f(i-1,j+1,k) + f(i-1,j-1,k)); fzz = f(i,j,k+1) - 2.0*f(i,j,k) + f(i,j,k-1);
fxz = 0.25*(f(i+1,j,k+1) - f(i+1,j,k-1) - f(i-1,j,k+1) + f(i-1,j,k-1)); fxy = 0.25*(f(i+1,j+1,k) - f(i+1,j-1,k) - f(i-1,j+1,k) + f(i-1,j-1,k));
fyz = 0.25*(f(i,j+1,k+1) - f(i,j+1,k-1) - f(i,j-1,k+1) + f(i,j-1,k-1)); fxz = 0.25*(f(i+1,j,k+1) - f(i+1,j,k-1) - f(i-1,j,k+1) + f(i-1,j,k-1));
fyz = 0.25*(f(i,j+1,k+1) - f(i,j+1,k-1) - f(i,j-1,k+1) + f(i,j-1,k-1));
// solid distance function // solid distance function
sx = 0.5*(s(i+1,j,k) - s(i-1,j,k)); sx = 0.5*(s(i+1,j,k) - s(i-1,j,k));
sy = 0.5*(s(i,j+1,k) - s(i,j-1,k)); sy = 0.5*(s(i,j+1,k) - s(i,j-1,k));
sz = 0.5*(s(i,j,k+1) - s(i,j,k-1)); sz = 0.5*(s(i,j,k+1) - s(i,j,k-1));
sxx = s(i+1,j,k) - 2.0*s(i,j,k) + s(i-1,j,k); sxx = s(i+1,j,k) - 2.0*s(i,j,k) + s(i-1,j,k);
syy = s(i,j+1,k) - 2.0*s(i,j,k) + s(i,j-1,k); syy = s(i,j+1,k) - 2.0*s(i,j,k) + s(i,j-1,k);
szz = s(i,j,k+1) - 2.0*s(i,j,k) + s(i,j,k-1); szz = s(i,j,k+1) - 2.0*s(i,j,k) + s(i,j,k-1);
sxy = 0.25*(s(i+1,j+1,k) - s(i+1,j-1,k) - s(i-1,j+1,k) + s(i-1,j-1,k)); sxy = 0.25*(s(i+1,j+1,k) - s(i+1,j-1,k) - s(i-1,j+1,k) + s(i-1,j-1,k));
sxz = 0.25*(s(i+1,j,k+1) - s(i+1,j,k-1) - s(i-1,j,k+1) + s(i-1,j,k-1)); sxz = 0.25*(s(i+1,j,k+1) - s(i+1,j,k-1) - s(i-1,j,k+1) + s(i-1,j,k-1));
syz = 0.25*(s(i,j+1,k+1) - s(i,j+1,k-1) - s(i,j-1,k+1) + s(i,j-1,k-1)); syz = 0.25*(s(i,j+1,k+1) - s(i,j+1,k-1) - s(i,j-1,k+1) + s(i,j-1,k-1));
// Compute the Jacobean matrix for tangent vector // Compute the Jacobean matrix for tangent vector
Axx = sxy*fz + sy*fxz - sxz*fy - sz*fxy; Axx = sxy*fz + sy*fxz - sxz*fy - sz*fxy;
Axy = sxz*fx + sz*fxx - sxx*fz - sx*fxz; Axy = sxz*fx + sz*fxx - sxx*fz - sx*fxz;
Axz = sxx*fy + sx*fxy - sxy*fx - sy*fxx; Axz = sxx*fy + sx*fxy - sxy*fx - sy*fxx;
Ayx = syy*fz + sy*fyz - syz*fy - sz*fyy; Ayx = syy*fz + sy*fyz - syz*fy - sz*fyy;
Ayy = syz*fx + sz*fxy - sxy*fz - sx*fyz; Ayy = syz*fx + sz*fxy - sxy*fz - sx*fyz;
Ayz = sxy*fy + sx*fyy - syy*fx - sy*fxy; Ayz = sxy*fy + sx*fyy - syy*fx - sy*fxy;
Axx = syz*fz + sy*fzz - szz*fy - sz*fyz; Axx = syz*fz + sy*fzz - szz*fy - sz*fyz;
Axy = szz*fx + sz*fxz - sxz*fz - sx*fzz; Axy = szz*fx + sz*fxz - sxz*fz - sx*fzz;
Axz = sxz*fy + sx*fyz - syz*fx - sy*fxz; Axz = sxz*fy + sx*fyz - syz*fx - sy*fxz;
// Compute the tangent vector
Tx.Corners(ic-i,jc-j,kc-k) = sy*fz-sz*fy;
Ty.Corners(ic-i,jc-j,kc-k) = sz*fx-sx*fz;
Tz.Corners(ic-i,jc-j,kc-k) = sx*fy-sy*fx;
// Compute the normal
Nx.Corners(ic-i,jc-j,kc-k) = Tx.Corners(ic-i,jc-j,kc-k)*Axx + Ty.Corners(ic-i,jc-j,kc-k)*Ayx + Tz.Corners(ic-i,jc-j,kc-k)*Azx;
Ny.Corners(ic-i,jc-j,kc-k) = Tx.Corners(ic-i,jc-j,kc-k)*Axy + Ty.Corners(ic-i,jc-j,kc-k)*Ayy + Tz.Corners(ic-i,jc-j,kc-k)*Azy;
Nz.Corners(ic-i,jc-j,kc-k) = Tx.Corners(ic-i,jc-j,kc-k)*Axz + Ty.Corners(ic-i,jc-j,kc-k)*Ayz + Tz.Corners(ic-i,jc-j,kc-k)*Azz;
}
}
}
/* for (int p=0; p<8; p++){
// Compute the tangent vector // Compute the tangent vector
Tx[p] = sy*fz-sz*fy; Tx[p] = sy*fz-sz*fy;
Ty[p] = sz*fx-sx*fz; Ty[p] = sz*fx-sx*fz;
@ -4296,8 +4329,9 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, DTMutableList<Po
Nx[p] = Tx[p]*Axx + Ty[p]*Ayx + Tz[p]*Azx; Nx[p] = Tx[p]*Axx + Ty[p]*Ayx + Tz[p]*Azx;
Ny[p] = Tx[p]*Axy + Ty[p]*Ayy + Tz[p]*Azy; Ny[p] = Tx[p]*Axy + Ty[p]*Ayy + Tz[p]*Azy;
Nz[p] = Tx[p]*Axz + Ty[p]*Ayz + Tz[p]*Azz; Nz[p] = Tx[p]*Axz + Ty[p]*Ayz + Tz[p]*Azz;
}
}
*/
for (p=0; p<npts; p++){ for (p=0; p<npts; p++){
P = Points(p); P = Points(p);
x = P.x-1.0*i; x = P.x-1.0*i;