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

This commit will probably cause merge conflicts :(

Browse Source
This commit is contained in:
James E McClure 2015-03-20 12:25:23 -04:00
parent ba66da2d0f
commit 77b2170e81
4 changed files with 206 additions and 463 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

115
common/TwoPhase.h
View File

@ -208,9 +208,11 @@ public:
void ComputeDelPhi();
void ColorToSignedDistance(double Beta, double *ColorData, double *DistData);
void ComputeLocal();
void ComputeLocalBlob();
void Reduce();
void NonDimensionalize(double D, double viscosity, double IFT);
void PrintAll(int timestep);
int GetCubeLabel(int i, int j, int k);
};
@ -342,7 +344,6 @@ void TwoPhase::UpdateMeshValues(){
//...........................................................................
}
void TwoPhase::ComputeLocal(){
int i,j,k,n;
double delphi;
@ -440,6 +441,104 @@ void TwoPhase::ComputeLocal(){
}
}
void TwoPhase::ComputeLocalBlob(){
int i,j,k,n,label;
double delphi;
int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};
for (int c=0;c<ncubes;c++){
// Get cube from the list
i = cubeList(0,c);
j = cubeList(1,c);
k = cubeList(2,c);
label=GetCubeLabel(i,j,k);
n_nw_pts=n_ns_pts=n_ws_pts=n_nws_pts=n_local_sol_pts=n_local_nws_pts=0;
n_nw_tris=n_ns_tris=n_ws_tris=n_nws_seg=n_local_sol_tris=0;
//...........................................................................
//...........................................................................
// Compute volume averages
for (int p=0;p<8;p++){
if ( SDs(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
// 1-D index for this cube corner
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
// compute the norm of the gradient of the phase indicator field
// Compute the non-wetting phase volume contribution
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
nwp_volume += 0.125;
// volume the excludes the interfacial region
if (DelPhi.data[n] < 1e-4){
vol_n += 0.125;
// pressure
pan += 0.125*Press.data[n];
// velocity
van(0) += 0.125*Vel_x.data[n];
van(1) += 0.125*Vel_y.data[n];
van(2) += 0.125*Vel_z.data[n];
}
}
else{
wp_volume += 0.125;
if (DelPhi.data[n] < 1e-4){
// volume the excludes the interfacial region
vol_w += 0.125;
// pressure
paw += 0.125*Press.data[n];
// velocity
vaw(0) += 0.125*Vel_x.data[n];
vaw(1) += 0.125*Vel_y.data[n];
vaw(2) += 0.125*Vel_z.data[n];
}
}
}
}
//...........................................................................
// Construct the interfaces and common curve
pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue,
nw_pts, nw_tris, Values, ns_pts, ns_tris, ws_pts, ws_tris,
local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
i, j, k, Nx, Ny, Nz);
// Integrate the contact angle
efawns += pmmc_CubeContactAngle(CubeValues,Values,SDn_x,SDn_y,SDn_z,SDs_x,SDs_y,SDs_z,
local_nws_pts,i,j,k,n_local_nws_pts);
// Integrate the mean curvature
Jwn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
Kwn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
// Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels)
pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues,
i,j,k,n_nw_pts,n_nw_tris,trimdist,trawn,trJwn);
pmmc_CubeTrimSurfaceInterpInverseValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues,
i,j,k,n_nw_pts,n_nw_tris,trimdist,dummy,trRwn);
// Compute the normal speed of the interface
pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris,
NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z,SDs_x,SDs_y,SDs_z,
local_nws_pts,i,j,k,n_local_nws_pts);
pmmc_CurveCurvature(SDn, SDs, KNwns_values, KGwns_values, KNwns, KGwns,
nws_pts, n_nws_pts, i, j, k);
As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris);
// Compute the surface orientation and the interfacial area
awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris);
ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris);
aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris);
lwns += pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
//...........................................................................
}
}
void TwoPhase::Reduce(){
int i;
double iVol_global=1.0/Volume;
@ -552,3 +651,17 @@ void TwoPhase::PrintAll(int timestep){
}
}
inline int TwoPhase::GetCubeLabel(int i, int j, int k){
int label;
label=LocalBlobID(i,j,k);
label=max(label,LocalBlobID(i+1,j,k));
label=max(label,LocalBlobID(i,j+1,k));
label=max(label,LocalBlobID(i+1,j+1,k));
label=max(label,LocalBlobID(i,j,k+1));
label=max(label,LocalBlobID(i+1,j,k+1));
label=max(label,LocalBlobID(i,j+1,k+1));
label=max(label,LocalBlobID(i+1,j+1,k+1));
return label;
}

408
common/pmmc.h
View File

@ -2066,42 +2066,7 @@ inline void MC( DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, i
}
}
}
// if point A(i,j,k) is in the solid phase
// else if ( A(i,j,k) == 0 ){
// pt.x = i;
// pt.y = j;
// pt.z = k;
// val = EXTRAP(A, v, i+1,j,k, 1,pt);
// // If extrapolated value gives a vertex
// if ( (A(i+1,j,k)- v)*(val-v) < 0 ){
// P.x = i + (val-v)/(val-A(i+1,j,k));
// P.y = j;
// P.z = k;
// if ( INTERP(solid(i,j,k), solid(i+1,j,k), P.x-i) > 0 ){
// nw_pts(n_nw_pts++) = P;
// N++;
// }
// }
// }
// // if point A(i+1,j,k) is in the solid phase
// else if ( A(i+1,j,k) == 0 ){
// pt.x = i+1;
// pt.y = j;
// pt.z = k;
// val = EXTRAP(A, v, i,j,k, 4,pt);
// // If extrapolated value gives a vertex
// if ( (A(i,j,k)- v)*(val-v) < 0 ){
// P.x = i + (A(i,j,k)-v)/(A(i,j,k)-val);
// P.y = j;
// P.z = k;
// if ( INTERP(solid(i,j,k), solid(i+1,j,k), P.x-i) > 0 ){
// nw_pts(n_nw_pts++) = P;
// N++;
// }
// }
// }
// }
// 2
if ((A(i+1,j,k)-v)*(A(i+1,j+1,k)-v) < 0)
{
if ( A(i+1,j,k) != 0 && A(i+1,j+1,k) != 0 ){
@ -2118,42 +2083,7 @@ inline void MC( DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, i
}
}
}
// else if ( A(i+1,j,k) == 0 ){
// pt.x = i+1;
// pt.y = j;
// pt.z = k;
// val = EXTRAP(A, v, i+1,j+1,k, 2,pt);
// // If extrapolated value gives a vertex
// if ( (A(i+1,j+1,k)- v)*(val-v) < 0 ){
// P.x = i+1;
// P.y = j + (val-v)/(val-A(i+1,j+1,k));
// P.z = k;
// if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
// INTERP(solid(i+1,j,k), solid(i+1,j+1,k), P.y-j) > 0 ){ // P is a new vertex (not counted twice)
// nw_pts(n_nw_pts++) = P;
// N++;
// }
// }
// }
// else if ( A(i+1,j+1,k) == 0){
// pt.x = i+1;
// pt.y = j+1;
// pt.z = k;
// val = EXTRAP(A, v, i+1,j,k, 5,pt);
// // If extrapolated value gives a vertex
// if ( (A(i+1,j,k)- v)*(val-v) < 0 ){
// P.x = i+1;
// P.y = j + (A(i+1,j,k)-v)/(A(i+1,j,k)-val);
// P.z = k;
// if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
// INTERP(solid(i+1,j,k), solid(i+1,j+1,k), P.y-j) > 0 ){
// nw_pts(n_nw_pts++) = P;
// N++;
// }
// }
// }
// }
//3
if ((A(i+1,j+1,k)-v)*(A(i,j+1,k)-v) < 0 )
{
if ( A(i+1,j+1,k) != 0 && A(i,j+1,k) != 0 ){
@ -2170,42 +2100,6 @@ inline void MC( DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, i
}
}
}
/* else if ( A(i,j+1,k) == 0 ){
pt.x = i;
pt.y = j+1;
pt.z = k;
val = EXTRAP(A, v, i+1,j+1,k, 1,pt);
// If extrapolated value gives a vertex
if ( (A(i+1,j+1,k)- v)*(val-v) < 0 ){
P.x = i + (val-v) / (val-A(i+1,j+1,k));
P.y = j+1;
P.z = k;
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j+1,k), solid(i+1,j+1,k), P.x-i) > 0 ){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
else if ( A(i+1,j+1,k) == 0){
pt.x = i+1;
pt.y = j+1;
pt.z = k;
val = EXTRAP(A, v, i,j+1,k, 4, pt);
// If extrapolated value gives a vertex
if ( (A(i,j+1,k)- v)*(val-v) < 0 ){
P.x = i + (A(i,j+1,k)-v) / (A(i,j+1,k)-val);
P.y = j+1;
P.z = k;
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j+1,k), solid(i+1,j+1,k), P.x-i) > 0 ){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
}
*/
//4
if ((A(i,j+1,k)-v)*(A(i,j,k)-v) < 0 )
@ -2224,41 +2118,7 @@ else if ( A(i+1,j+1,k) == 0){
}
}
}
/* else if ( A(i,j+1,k) == 0 ){
pt.x = i;
pt.y = j+1;
pt.z = k;
val = EXTRAP(A, v, i,j,k, 5,pt);
// If extrapolated value gives a vertex
if ( (A(i,j,k)- v)*(val-v) < 0 ){
P.x = i;
P.y = j + (A(i,j,k)-v) / (A(i,j,k)-val);
P.z = k;
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j,k), solid(i,j+1,k), P.y-j) > 0 ){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
}
}
}
else if ( A(i,j,k) == 0){
pt.x = i;
pt.y = j;
pt.z = k;
val = EXTRAP(A, v, i,j+1,k, 2,pt);
// If extrapolated value gives a vertex
if ( (A(i,j+1,k)- v)*(val-v) < 0 ){
P.x = i;
P.y = j + (val-v) / (val-A(i,j+1,k));
P.z = k;
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j,k), solid(i,j+1,k), P.y-j) > 0){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
}
*/
//5
if ((A(i,j,k)-v)*(A(i,j,k+1)-v) < 0 )
{
@ -2276,42 +2136,7 @@ else if ( A(i,j,k) == 0){
}
}
}
/* else if ( A(i,j,k) == 0 ){
pt.x = i;
pt.y = j;
pt.z = k;
val = EXTRAP(A, v, i,j,k+1, 3,pt);
// If extrapolated value gives a vertex
if ( (A(i,j,k+1)- v)*(val-v) < 0 ){
P.x = i;
P.y = j;
P.z = k + (val-v) / (val-A(i,j,k+1));
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j,k), solid(i,j,k+1), P.z-k) > 0){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
else if ( A(i,j,k+1) == 0){
pt.x = i;
pt.y = j;
pt.z = k+1;
val = EXTRAP(A, v, i,j,k, 6,pt);
// If extrapolated value gives a vertex
if ( (A(i,j,k)- v)*(val-v) < 0 ){
P.x = i;
P.y = j;
P.z = k + (A(i,j,k)-v) / (A(i,j,k)-val);
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j,k), solid(i,j,k+1), P.z-k) > 0){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
}
*/
//6
if ((A(i+1,j,k)-v)*(A(i+1,j,k+1)-v) < 0 )
{
@ -2329,43 +2154,7 @@ else if ( A(i,j,k+1) == 0){
}
}
}
/* else if ( A(i+1,j,k) == 0 ){
pt.x = i+1;
pt.y = j;
pt.z = k;
val = EXTRAP(A, v, i+1,j,k+1, 3,pt);
// If extrapolated value gives a vertex
if ( (A(i+1,j,k+1)- v)*(val-v) < 0 ){
P.x = i+1;
P.y = j;
P.z = k + (val-v) / (val-A(i+1,j,k+1));
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i+1,j,k), solid(i+1,j,k+1), P.z-k) > 0){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
else if ( A(i+1,j,k+1) == 0){
pt.x = i+1;
pt.y = j;
pt.z = k+1;
val = EXTRAP(A, v, i+1,j,k, 6,pt);
// If extrapolated value gives a vertex
if ( (A(i+1,j,k)- v)*(val-v) < 0 ){
P.x = i+1;
P.y = j;
P.z = k + (A(i+1,j,k)-v) / (A(i+1,j,k)-val);
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i+1,j,k), solid(i+1,j,k+1), P.z-k) > 0){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
}
*/
//7
if ((A(i+1,j+1,k)-v)*(A(i+1,j+1,k+1)-v) < 0 )
{
@ -2383,42 +2172,7 @@ else if ( A(i+1,j,k+1) == 0){
}
}
}
/* else if ( A(i+1,j+1,k) == 0 ){
pt.x = i+1;
pt.y = j+1;
pt.z = k;
val = EXTRAP(A, v, i+1,j+1,k+1, 3,pt);
// If extrapolated value gives a vertex
if ( (A(i+1,j+1,k+1)- v)*(val-v) < 0 ){
P.x = i+1;
P.y = j+1;
P.z = k + (val-v) / (val-A(i+1,j+1,k+1));
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i+1,j+1,k), solid(i+1,j+1,k+1), P.z-k) > 0 ){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
else if ( A(i+1,j+1,k+1) == 0){
pt.x = i+1;
pt.y = j+1;
pt.z = k+1;
val = EXTRAP(A, v, i+1,j+1,k, 6,pt);
// If extrapolated value gives a vertex
if ( (A(i+1,j+1,k)- v)*(val-v) < 0 ){
P.x = i+1;
P.y = j+1;
P.z = k + (A(i+1,j+1,k)-v) / (A(i+1,j+1,k)-val);
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i+1,j+1,k), solid(i+1,j+1,k+1), P.z-k) > 0 ){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
}
*/
//8
if ((A(i,j+1,k)-v)*(A(i,j+1,k+1)-v) < 0 )
{
@ -2436,42 +2190,7 @@ else if ( A(i+1,j+1,k+1) == 0){
}
}
}
/* else if ( A(i,j+1,k) == 0 ){
pt.x = i;
pt.y = j+1;
pt.z = k;
val = EXTRAP(A, v, i,j+1,k+1, 3,pt);
// If extrapolated value gives a vertex
if ( (A(i,j+1,k+1)- v)*(val-v) < 0 ){
P.x = i;
P.y = j+1;
P.z = k + (val-v) / (val-A(i,j+1,k+1));
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j+1,k), solid(i,j+1,k+1), P.z-k) > 0 ){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
else if ( A(i,j+1,k+1) == 0){
pt.x = i;
pt.y = j+1;
pt.z = k+1;
val = EXTRAP(A, v, i,j+1,k, 6,pt);
// If extrapolated value gives a vertex
if ( (A(i,j+1,k)- v)*(val-v) < 0 ){
P.x = i;
P.y = j+1;
P.z = k + (A(i,j+1,k)-v) / (A(i,j+1,k)-val);
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j+1,k), solid(i,j+1,k+1), P.z-k) > 0 ){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
}
*/
//9
if ((A(i,j,k+1)-v)*(A(i+1,j,k+1)-v) < 0 )
{
@ -2489,42 +2208,7 @@ else if ( A(i,j+1,k+1) == 0){
}
}
}
/* else if ( A(i,j,k+1) == 0 ){
pt.x = i;
pt.y = j;
pt.z = k+1;
val = EXTRAP(A, v, i+1,j,k+1, 1,pt);
// If extrapolated value gives a vertex
if ( (A(i+1,j,k+1)- v)*(val-v) < 0 ){
P.x = i + (val-v) / (val-A(i+1,j,k+1));
P.y = j;
P.z = k+1;
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j,k+1), solid(i+1,j,k+1), P.x-i) > 0 ){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
else if ( A(i+1,j,k+1) == 0){
pt.x = i+1;
pt.y = j;
pt.z = k+1;
val = EXTRAP(A, v, i,j,k+1, 4,pt);
// If extrapolated value gives a vertex
if ( (A(i,j,k+1)- v)*(val-v) < 0 ){
P.x = i + (A(i,j,k+1)-v) / (A(i,j,k+1)-val);
P.y = j;
P.z = k+1;
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j,k+1), solid(i+1,j,k+1), P.x-i) > 0){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
}
*/
//10
if ((A(i+1,j,k+1)-v)*(A(i+1,j+1,k+1)-v) < 0 )
{
@ -2542,42 +2226,7 @@ else if ( A(i+1,j,k+1) == 0){
}
}
}
/* else if ( A(i+1,j,k+1) == 0 ){
pt.x = i+1;
pt.y = j;
pt.z = k+1;
val = EXTRAP(A, v, i+1,j+1,k+1, 2,pt);
// If extrapolated value gives a vertex
if ( (A(i+1,j+1,k+1)- v)*(val-v) < 0 ){
P.x = i+1;
P.y = j + (val-v) / (val-A(i+1,j+1,k+1));
P.z = k+1;
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i+1,j,k+1), solid(i+1,j+1,k+1), P.y-j) > 0){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
else if ( A(i+1,j+1,k+1) == 0){
pt.x = i+1;
pt.y = j+1;
pt.z = k+1;
val = EXTRAP(A, v, i+1,j,k+1, 5,pt);
// If extrapolated value gives a vertex
if ( (A(i+1,j,k+1)- v)*(val-v) < 0 ){
P.x = i+1;
P.y = j + (A(i+1,j,k+1)-v) / (A(i+1,j,k+1)-val);
P.z = k+1;
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i+1,j,k+1), solid(i+1,j+1,k+1), P.y-j) > 0){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
}
*/
//11
if ((A(i+1,j+1,k+1)-v)*(A(i,j+1,k+1)-v) < 0 )
{
@ -2595,42 +2244,7 @@ else if ( A(i+1,j+1,k+1) == 0){
}
}
}
/* else if ( A(i+1,j+1,k+1) == 0 ){
pt.x = i+1;
pt.y = j+1;
pt.z = k+1;
val = EXTRAP(A, v, i,j+1,k+1, 4,pt);
// If extrapolated value gives a vertex
if ( (A(i,j+1,k+1)- v)*(val-v) < 0 ){
P.x = i+(A(i,j+1,k+1)-v) / (A(i,j+1,k+1)-val);
P.y = j+1;
P.z = k+1;
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j+1,k+1), solid(i+1,j+1,k+1), P.x-i) > 0){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
else if ( A(i,j+1,k+1) == 0){
pt.x = i;
pt.y = j+1;
pt.z = k+1;
val = EXTRAP(A, v, i+1,j+1,k+1, 1,pt);
// If extrapolated value gives a vertex
if ( (A(i+1,j+1,k+1)- v)*(val-v) < 0 ){
P.x = i+(val-v) / (val-A(i+1,j+1,k+1));
P.y = j+1;
P.z = k+1;
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1 &&
INTERP(solid(i,j+1,k+1), solid(i+1,j+1,k+1), P.x-i) > 0){ // P is a new vertex (not counted twice)
nw_pts(n_nw_pts++) = P;
N++;
}
}
}
}
*/
//12
if ((A(i,j+1,k+1)-v)*(A(i,j,k+1)-v) < 0 )
{
@ -2657,6 +2271,8 @@ else if ( A(i,j+1,k+1) == 0){
// n = number of vertices in this grid cell
// Assemble the triangles as long as points are found
if (N > 0){
for (m = n_nw_pts-N; m < n_nw_pts-2; m++) {
for (o = m+2; o < n_nw_pts-1; o++) {
if (ShareSide(nw_pts(m), nw_pts(o)) == 1) {
@ -2716,7 +2332,7 @@ else if ( A(i,j+1,k+1) == 0){
nw_tris(2,n_nw_tris) = p;
n_nw_tris++;
}
}
}
//-------------------------------------------------------------------------------
inline void EDGE(DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k, int &m, int &n, int &o,

1
tests/BlobAnalysis.cpp
View File

@ -12,7 +12,6 @@
using namespace std;
inline void ReadCheckpoint(char *FILENAME, double *cDen, double *cDistEven, double *cDistOdd, int N)
{
int q,n;

29
tests/lbpm_BlobAnalysis.cpp
View File

@ -55,9 +55,6 @@ int main(int argc, char **argv)
int BC; // type of boundary condition applied: 0-periodic, 1-pressure/velocity
int nblobs_global; // number of blobs in the global system
int *CubeList;
if (rank==0){
//.......................................................................
// Reading the domain information file
@ -166,10 +163,13 @@ int main(int argc, char **argv)
Averages.Vel_z.data[n]=vz;
}
// Count the number of cubes for each blob
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=1;i<Nx;i++){
int label;
DoubleArray BlobAverages(50,n_blobs_global);
for (k=1;k<Nz-1;k++){
for (j=1;j<Ny-1;j++){
for (i=1;i<Nx-1;i++){
// Assign the label for the cube
label = Averages.GetCubeLabel(i,j,k);
}
}
@ -188,3 +188,18 @@ int main(int argc, char **argv)
MPI_Finalize();
// ****************************************************
}
inline int GetCubeLabel(int i, int j, int k){
int label;
label=BlobLabel(i,j,k);
label=max(label,BlobLabel(i+1,j,k));
label=max(label,BlobLabel(i,j+1,k));
label=max(label,BlobLabel(i+1,j+1,k));
label=max(label,BlobLabel(i,j,k+1));
label=max(label,BlobLabel(i+1,j,k+1));
label=max(label,BlobLabel(i,j+1,k+1));
label=max(label,BlobLabel(i+1,j+1,k+1));
return label;
}