diff --git a/common/TwoPhase.h b/common/TwoPhase.h index 20bc7dcc..1fe0e926 100644 --- a/common/TwoPhase.h +++ b/common/TwoPhase.h @@ -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 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; +} + diff --git a/common/pmmc.h b/common/pmmc.h index 434f9233..3581fc14 100644 --- a/common/pmmc.h +++ b/common/pmmc.h @@ -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,66 +2271,68 @@ else if ( A(i,j+1,k+1) == 0){ // n = number of vertices in this grid cell - 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; - } - } + // 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; - } - } + // 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