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

Added other minkowski measures (curvature, surface area)

Browse Source
This commit is contained in:
James E McClure 2015-09-04 16:13:32 -04:00
parent 941c3cb98f
commit 4c26a33038
3 changed files with 22 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
common/TwoPhase.cpp
View File

@ -367,14 +367,6 @@ void TwoPhase::ComputeLocal()
if (Dm.BoundaryCondition > 0 && Dm.kproc == 0) kmin=4;
if (Dm.BoundaryCondition > 0 && Dm.kproc == Dm.nprocz-1) kmax=Nz-4;
// Map solid to erode the fluid so that interfaces can be calculated accurately
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
SDs(i,j,k) += 1.0;
}
}
}
for (k=kmin; k<kmax; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
@ -472,7 +464,11 @@ void TwoPhase::ComputeLocal()
//...........................................................................
// Compute the integral curvature of the non-wetting phase
n_nw_pts=n_nw_tris=0;
geomavg_MarchingCubes(SDn,fluid_isovalue,i,j,k,nw_pts,n_nw_pts,nw_tris,n_nw_tris);
// Compute the non-wetting phase surface and associated area
An += geomavg_MarchingCubes(SDn,fluid_isovalue,i,j,k,nw_pts,n_nw_pts,nw_tris,n_nw_tris);
// Compute the integral of mean curvature
Jn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,
i,j,k,n_nw_pts,n_nw_tris);
// Compute Euler characteristic from integral of gaussian curvature
euler += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,
i,j,k,n_nw_pts,n_nw_tris);
@ -482,14 +478,6 @@ void TwoPhase::ComputeLocal()
}
}
// Map solid back
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
SDs(i,j,k) -= 1.0;
}
}
}
}

18
common/pmmc.h
View File

@ -1,7 +1,6 @@
#ifndef pmmc_INC
#define pmmc_INC
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
@ -1836,7 +1835,7 @@ inline void TRIM(DTMutableList<Point> &local_sol_pts, int &n_local_sol_pts, doub
}
}
}
inline void geomavg_MarchingCubes( DoubleArray &A, double &v, int &i, int &j, int &k,
inline double geomavg_MarchingCubes( DoubleArray &A, double &v, int &i, int &j, int &k,
DTMutableList<Point> &nw_pts, int &n_nw_pts, IntArray &nw_tris,
int &n_nw_tris)
{
@ -1847,6 +1846,7 @@ inline void geomavg_MarchingCubes( DoubleArray &A, double &v, int &i, int &j, in
int m;
int o;
int p;
int area;
// Go over each corner -- check to see if the corners are themselves vertices
//1
@ -2149,6 +2149,20 @@ inline void geomavg_MarchingCubes( DoubleArray &A, double &v, int &i, int &j, in
n_nw_tris++;
}
}
// Compute the Interfacial Area
double s1,s2,s3,s;
for (int r=n_nw_tris_beg;r<n_nw_tris;r++){
A = nw_pts(nw_tris(0,r));
B = nw_pts(nw_tris(1,r));
C = nw_pts(nw_tris(2,r));
// Compute length of sides (assume dx=dy=dz)
s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
area+=sqrt(s*(s-s1)*(s-s2)*(s-s3));
}
return area;
}
//-------------------------------------------------------------------------------
inline void MC( DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k,

2
tests/lbpm_color_simulator.cpp
View File

@ -576,7 +576,7 @@ int main(int argc, char **argv)
MPI_Barrier(comm);
//.......................................................................
// Once phase has been initialized, map solid to account for 'smeared' interface
//for (i=0; i<N; i++) Averages->SDs(i) -= (1.0); //
for (i=0; i<N; i++) Averages->SDs(i) -= (1.0); //
//.......................................................................
// Finalize setup for averaging domain
//Averages->SetupCubes(Dm);