Continuing to refactor Euler characteristic
This commit is contained in:
James E McClure
parent
3a4e564ca6
commit
df0a3f5635
@ -507,7 +507,7 @@ void TwoPhase::ComponentAverages()
|
||||
for (i=1; i<Nx-1; i++){
|
||||
LabelWP=GetCubeLabel(i,j,k,Label_WP);
|
||||
LabelNWP=GetCubeLabel(i,j,k,Label_NWP);
|
||||
|
||||
|
||||
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;
|
||||
|
||||
@ -592,12 +592,12 @@ void TwoPhase::ComponentAverages()
|
||||
*/
|
||||
if (n_nw_pts+n_ns_pts > 0){
|
||||
|
||||
/* double euler;
|
||||
/* double euler;
|
||||
euler = geomavg_EulerCharacteristic(nw_pts,nw_tris,n_nw_pts,n_nw_tris,i,j,k);
|
||||
euler += geomavg_EulerCharacteristic(ns_pts,ns_tris,n_ns_pts,n_ns_tris,i,j,k);
|
||||
// adjust for double-counted vertices and edges from the common curve
|
||||
if (n_nws_pts > 0) euler += 1.0*n_nws_pts;
|
||||
*/
|
||||
|
||||
int nvert = n_nw_pts+n_ns_pts-n_nws_pts;
|
||||
int nside = 2*nvert-3;
|
||||
int nface = nvert-2;
|
||||
@ -678,6 +678,8 @@ void TwoPhase::ComponentAverages()
|
||||
//ComponentAverages_NWP(NVERT,LabelNWP) += nvert;
|
||||
//ComponentAverages_NWP(NSIDE,LabelNWP) += nside;
|
||||
//ComponentAverages_NWP(NFACE,LabelNWP) += nface;
|
||||
*
|
||||
*/
|
||||
}
|
||||
|
||||
//...........................................................................
|
||||
@ -778,6 +780,54 @@ void TwoPhase::ComponentAverages()
|
||||
}
|
||||
//...........................................................................
|
||||
|
||||
// Compute the Euler characteristic
|
||||
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);
|
||||
if (n_nw_pts > 0 ){
|
||||
int nvert = n_nw_pts;
|
||||
int nside = 2*nvert-3;
|
||||
int nface = nvert-2;
|
||||
//...........................................................
|
||||
// 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<n_nw_pts; p++){
|
||||
Point PT = nw_pts(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 previously computed edges
|
||||
for (int p=0; p<n_nw_tris; p++){
|
||||
Point A = nw_pts(nw_tris(0,p));
|
||||
Point B = nw_pts(nw_tris(1,p));
|
||||
Point C = nw_pts(nw_tris(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;
|
||||
|
||||
}
|
||||
int euler=nvert-nside+nface; // euler characteristic for the cube
|
||||
ComponentAverages_NWP(EULER,LabelNWP) += 1.0*euler;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
322
common/pmmc.h
322
common/pmmc.h
@ -1837,6 +1837,320 @@ 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,
|
||||
DTMutableList<Point> &nw_pts, int &n_nw_pts, IntArray &nw_tris,
|
||||
int &n_nw_tris)
|
||||
{
|
||||
int N = 0; // n will be the number of vertices in this grid cell only
|
||||
Point P;
|
||||
Point pt;
|
||||
Point PlaceHolder;
|
||||
int m;
|
||||
int o;
|
||||
int p;
|
||||
|
||||
// Go over each corner -- check to see if the corners are themselves vertices
|
||||
//1
|
||||
if (A(i,j,k) == v){
|
||||
P.x = i;
|
||||
P.y = j;
|
||||
P.z = k;
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
//2
|
||||
if (A(i+1,j,k) == v){
|
||||
P.x = i+1;
|
||||
P.y = j;
|
||||
P.z = k;
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
//3
|
||||
if (A(i+1,j+1,k) == v){
|
||||
P.x = i+1;
|
||||
P.y = j+1;
|
||||
P.z = k;
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
//4
|
||||
if (A(i,j+1,k) == v){
|
||||
P.x = i;
|
||||
P.y = j+1;
|
||||
P.z = k;
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
//5
|
||||
if (A(i,j,k+1) == v){
|
||||
P.x = i;
|
||||
P.y = j;
|
||||
P.z = k+1;
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
//6
|
||||
if (A(i+1,j,k+1) == v){
|
||||
P.x = i+1;
|
||||
P.y = j;
|
||||
P.z = k+1;
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
//7
|
||||
if (A(i+1,j+1,k+1) == v){
|
||||
P.x = i+1;
|
||||
P.y = j+1;
|
||||
P.z = k+1;
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
//8
|
||||
if (A(i,j+1,k+1) == v){
|
||||
P.x = i;
|
||||
P.y = j+1;
|
||||
P.z = k+1;
|
||||
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
|
||||
// Go through each side, compute P for sides of box spiraling up
|
||||
|
||||
|
||||
// float val;
|
||||
if ((A(i,j,k)-v)*(A(i+1,j,k)-v) < 0)
|
||||
{
|
||||
// If both points are in the fluid region
|
||||
if (A(i,j,k) != 0 && A(i+1,j,k) != 0){
|
||||
P.x = i + (A(i,j,k)-v)/(A(i,j,k)-A(i+1,j,k));
|
||||
P.y = j;
|
||||
P.z = k;
|
||||
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
}
|
||||
|
||||
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 ){
|
||||
P.x = i+1;
|
||||
P.y = j + (A(i+1,j,k)-v)/(A(i+1,j,k)-A(i+1,j+1,k));
|
||||
P.z = k;
|
||||
|
||||
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){ // P is a new vertex (not counted twice)
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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 ){
|
||||
P.x = i + (A(i,j+1,k)-v) / (A(i,j+1,k)-A(i+1,j+1,k));
|
||||
P.y = j+1;
|
||||
P.z = k;
|
||||
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){ // 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 )
|
||||
{
|
||||
if (A(i,j+1,k) != 0 && A(i,j,k) != 0 ){
|
||||
P.x = i;
|
||||
P.y = j + (A(i,j,k)-v) / (A(i,j,k)-A(i,j+1,k));
|
||||
P.z = k;
|
||||
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){ // 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 )
|
||||
{
|
||||
if ( A(i,j,k) != 0 && A(i,j,k+1) != 0 ){
|
||||
P.x = i;
|
||||
P.y = j;
|
||||
P.z = k + (A(i,j,k)-v) / (A(i,j,k)-A(i,j,k+1));
|
||||
|
||||
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){ // 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 )
|
||||
{
|
||||
if ( A(i+1,j,k) != 0 && A(i+1,j,k+1) != 0 ){
|
||||
P.x = i+1;
|
||||
P.y = j;
|
||||
P.z = k + (A(i+1,j,k)-v) / (A(i+1,j,k)-A(i+1,j,k+1));
|
||||
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
|
||||
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 )
|
||||
{
|
||||
if ( A(i+1,j+1,k) != 0 && A(i+1,j+1,k+1) != 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)-A(i+1,j+1,k+1));
|
||||
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//8
|
||||
if ((A(i,j+1,k)-v)*(A(i,j+1,k+1)-v) < 0 )
|
||||
{
|
||||
if ( A(i,j+1,k) != 0 && A(i,j+1,k+1) != 0 ){
|
||||
P.x = i;
|
||||
P.y = j+1;
|
||||
P.z = k + (A(i,j+1,k)-v) / (A(i,j+1,k)-A(i,j+1,k+1));
|
||||
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//9
|
||||
if ((A(i,j,k+1)-v)*(A(i+1,j,k+1)-v) < 0 )
|
||||
{
|
||||
if ( A(i,j,k+1) != 0 && A(i+1,j,k+1) != 0 ){
|
||||
P.x = i + (A(i,j,k+1)-v) / (A(i,j,k+1)-A(i+1,j,k+1));
|
||||
P.y = j;
|
||||
P.z = k+1;
|
||||
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
|
||||
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 )
|
||||
{
|
||||
if ( A(i+1,j,k+1) != 0 && A(i+1,j+1,k+1) != 0 ){
|
||||
P.x = i+1;
|
||||
P.y = j + (A(i+1,j,k+1)-v) / (A(i+1,j,k+1)-A(i+1,j+1,k+1));
|
||||
P.z = k+1;
|
||||
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
|
||||
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 )
|
||||
{
|
||||
if ( A(i+1,j+1,k+1) != 0 && A(i,j+1,k+1) != 0 ){
|
||||
P.x = i+(A(i,j+1,k+1)-v) / (A(i,j+1,k+1)-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){
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//12
|
||||
if ((A(i,j+1,k+1)-v)*(A(i,j,k+1)-v) < 0 )
|
||||
{
|
||||
if ( A(i,j+1,k+1) != 0 && A(i,j,k+1) != 0 ){
|
||||
P.x = i;
|
||||
P.y = j + (A(i,j,k+1)-v) / (A(i,j,k+1)-A(i,j+1,k+1));
|
||||
P.z = k+1;
|
||||
if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
|
||||
nw_pts(n_nw_pts++) = P;
|
||||
N++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// 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) {
|
||||
PlaceHolder = nw_pts(m+1);
|
||||
nw_pts(m+1) = nw_pts(o);
|
||||
nw_pts(o) = PlaceHolder;
|
||||
}
|
||||
}
|
||||
|
||||
// make sure other neighbor of vertex 1 is in last spot
|
||||
if (m == n_nw_pts-N){
|
||||
for (p = m+2; p < n_nw_pts-1; p++){
|
||||
if (ShareSide(nw_pts(m), nw_pts(p)) == 1){
|
||||
PlaceHolder = nw_pts(n_nw_pts-1);
|
||||
nw_pts(n_nw_pts-1) = nw_pts(p);
|
||||
nw_pts(p) = PlaceHolder;
|
||||
}
|
||||
}
|
||||
}
|
||||
if ( ShareSide(nw_pts(n_nw_pts-2), nw_pts(n_nw_pts-3)) != 1 ){
|
||||
if (ShareSide( nw_pts(n_nw_pts-3), nw_pts(n_nw_pts-1)) == 1 &&
|
||||
ShareSide( nw_pts(n_nw_pts-N),nw_pts(n_nw_pts-2)) == 1 ){
|
||||
PlaceHolder = nw_pts(n_nw_pts-2);
|
||||
nw_pts(n_nw_pts-2) = nw_pts(n_nw_pts-1);
|
||||
nw_pts(n_nw_pts-1) = PlaceHolder;
|
||||
}
|
||||
}
|
||||
if ( ShareSide(nw_pts(n_nw_pts-1), nw_pts(n_nw_pts-2)) != 1 ){
|
||||
if (ShareSide( nw_pts(n_nw_pts-3), nw_pts(n_nw_pts-1)) == 1 &&
|
||||
ShareSide(nw_pts(n_nw_pts-4),nw_pts(n_nw_pts-2)) == 1 ){
|
||||
PlaceHolder = nw_pts(n_nw_pts-3);
|
||||
nw_pts(n_nw_pts-3) = nw_pts(n_nw_pts-2);
|
||||
nw_pts(n_nw_pts-2) = PlaceHolder;
|
||||
}
|
||||
if (ShareSide( nw_pts(n_nw_pts-N+1), nw_pts(n_nw_pts-3)) == 1 &&
|
||||
ShareSide(nw_pts(n_nw_pts-1),nw_pts(n_nw_pts-N+1)) == 1 ){
|
||||
PlaceHolder = nw_pts(n_nw_pts-2);
|
||||
nw_pts(n_nw_pts-2) = nw_pts(n_nw_pts-N+1);
|
||||
nw_pts(n_nw_pts-N+1) = PlaceHolder;
|
||||
}
|
||||
}
|
||||
if ( ShareSide(nw_pts(n_nw_pts-N), nw_pts(n_nw_pts-N+1)) != 1 ){
|
||||
if (ShareSide( nw_pts(n_nw_pts-N), nw_pts(n_nw_pts-2)) == 1 &&
|
||||
ShareSide(nw_pts(n_nw_pts-1), nw_pts(n_nw_pts-N+1)) == 1){
|
||||
PlaceHolder = nw_pts(n_nw_pts-1);
|
||||
nw_pts(n_nw_pts-1) = nw_pts(n_nw_pts-N);
|
||||
nw_pts(n_nw_pts-N) = PlaceHolder;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// * * * ESTABLISH TRIANGLE CONNECTIONS * * *
|
||||
|
||||
for (p=n_nw_pts-N+2; p<n_nw_pts; p++){
|
||||
nw_tris(0,n_nw_tris) = n_nw_pts-N;
|
||||
nw_tris(1,n_nw_tris) = p-1;
|
||||
nw_tris(2,n_nw_tris) = p;
|
||||
n_nw_tris++;
|
||||
}
|
||||
}
|
||||
}
|
||||
//-------------------------------------------------------------------------------
|
||||
inline void MC( DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k,
|
||||
DTMutableList<Point> &nw_pts, int &n_nw_pts, IntArray &nw_tris,
|
||||
@ -4037,14 +4351,14 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
|
||||
}
|
||||
//--------------------------------------------------------------------------------------------------------
|
||||
inline double geomavg_EulerCharacteristic(DTMutableList<Point> &Points, IntArray &Triangles,
|
||||
int npts, int ntris, int i, int j, int k){
|
||||
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*npts-3;
|
||||
int nface=npts-2;
|
||||
int nside=2*vert-3;
|
||||
int nface=nvert-2;
|
||||
//if (ntris != nface){
|
||||
// nface = ntris;
|
||||
// nside =
|
||||
@ -4087,7 +4401,7 @@ inline double geomavg_EulerCharacteristic(DTMutableList<Point> &Points, IntArray
|
||||
if (!newside) nside-=1;
|
||||
|
||||
}
|
||||
EulerChar = double(nvert - nside + nface);
|
||||
EulerChar = 1.0*(nvert - nside + nface);
|
||||
return EulerChar;
|
||||
}
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user