From 5c794e0bd008080d54f8adfcc54bb451118a74dc Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 15 Dec 2021 07:35:10 -0500 Subject: [PATCH] fix sign on geodesic curvature --- analysis/TwoPhase.cpp | 30 ++++++++++++++++++++++++++++-- analysis/TwoPhase.h | 2 ++ analysis/pmmc.h | 20 +++++++++++++------- tests/lbpm_TwoPhase_analysis.cpp | 10 ++++++++++ 4 files changed, 53 insertions(+), 9 deletions(-) diff --git a/analysis/TwoPhase.cpp b/analysis/TwoPhase.cpp index 20f4220b..ebebf222 100644 --- a/analysis/TwoPhase.cpp +++ b/analysis/TwoPhase.cpp @@ -361,6 +361,7 @@ void TwoPhase::Initialize() { wwndnw = 0.0; wwnsdnwn = 0.0; Jwnwwndnw = 0.0; + Xwn = Xns = Xws = 0.0; } void TwoPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB, @@ -780,6 +781,9 @@ void TwoPhase::ComputeStatic() { Jwn += pmmc_CubeSurfaceInterpValue( CubeValues, MeanCurvature, nw_pts, nw_tris, Values, i, j, k, n_nw_pts, n_nw_tris); + + Xwn += geomavg_EulerCharacteristic(nw_pts, nw_tris, n_nw_pts, + n_nw_tris, i, j, k); // Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels) pmmc_CubeTrimSurfaceInterpValues( @@ -812,6 +816,14 @@ void TwoPhase::ComputeStatic() { lwns += pmmc_CubeCurveLength(local_nws_pts, n_local_nws_pts); + + /* half contribution for vertices / edges at the common line + * each cube with contact line has a net of undercounting vertices + * each cube is undercounting edges due to internal counts + */ + Xwn += 0.25*n_local_nws_pts - 0.5; + Xws += 0.25*n_local_nws_pts - 0.5; + Xns += 0.25*n_local_nws_pts - 0.5; } // Solid interface averagees @@ -824,6 +836,12 @@ void TwoPhase::ComputeStatic() { n_ns_tris); aws += pmmc_CubeSurfaceOrientation(Gws, ws_pts, ws_tris, n_ws_tris); + + Xws += geomavg_EulerCharacteristic(ws_pts, ws_tris, n_ws_pts, + n_ws_tris, i, j, k); + + Xns += geomavg_EulerCharacteristic(ns_pts, ns_tris, n_ns_pts, + n_ns_tris, i, j, k); } //........................................................................... // Compute the integral curvature of the non-wetting phase @@ -848,9 +866,11 @@ void TwoPhase::ComputeStatic() { Kn += pmmc_CubeSurfaceInterpValue(CubeValues, GaussCurvature, nw_pts, nw_tris, Values, i, j, k, n_nw_pts, n_nw_tris); + euler += geomavg_EulerCharacteristic(nw_pts, nw_tris, n_nw_pts, n_nw_tris, i, j, k); + } } } @@ -1447,10 +1467,14 @@ void TwoPhase::Reduce() { trawn_global = Dm->Comm.sumReduce(trawn); trJwn_global = Dm->Comm.sumReduce(trJwn); trRwn_global = Dm->Comm.sumReduce(trRwn); - euler_global = Dm->Comm.sumReduce(euler); + Xwn_global = Dm->Comm.sumReduce(Xwn); + Xws_global = Dm->Comm.sumReduce(Xws); + Xns_global = Dm->Comm.sumReduce(Xns); An_global = Dm->Comm.sumReduce(An); Jn_global = Dm->Comm.sumReduce(Jn); Kn_global = Dm->Comm.sumReduce(Kn); + euler_global = Dm->Comm.sumReduce(euler); + Dm->Comm.barrier(); // Normalize the phase averages @@ -1533,7 +1557,7 @@ void TwoPhase::PrintStatic() { if (fseek(STATIC, 0, SEEK_SET) == fseek(STATIC, 0, SEEK_CUR)) { // If timelog is empty, write a short header to list the averages fprintf(STATIC, "sw awn ans aws Jwn Kwn lwns cwns KGws " - "KGwn "); // Scalar averages + "KGwn Xwn Xws Xns "); // Scalar averages fprintf(STATIC, "Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors fprintf(STATIC, "Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz "); @@ -1553,6 +1577,8 @@ void TwoPhase::PrintStatic() { fprintf(STATIC, "%.5g ", efawns_global); // average contact angle fprintf(STATIC, "%.5g %.5g ", KNwns_global, KGwns_global); // total curvature contributions of common line + fprintf(STATIC, "%.5g %.5g %.5g ", Xwn_global, Xns_global, + Xws_global); // Euler characteristic fprintf(STATIC, "%.5g %.5g %.5g %.5g %.5g %.5g ", Gwn_global(0), Gwn_global(1), Gwn_global(2), Gwn_global(3), Gwn_global(4), Gwn_global(5)); // orientation of wn interface diff --git a/analysis/TwoPhase.h b/analysis/TwoPhase.h index 535b763d..15247431 100644 --- a/analysis/TwoPhase.h +++ b/analysis/TwoPhase.h @@ -103,7 +103,9 @@ public: double lwns_global; double efawns, efawns_global; // averaged contact angle double euler, Kn, Jn, An; + double Xwn, Xns, Xws; double euler_global, Kn_global, Jn_global, An_global; + double Xwn_global, Xns_global, Xws_global; double rho_n, rho_w; double nu_n, nu_w; diff --git a/analysis/pmmc.h b/analysis/pmmc.h index 8dd964f5..5e07cbe8 100644 --- a/analysis/pmmc.h +++ b/analysis/pmmc.h @@ -4572,14 +4572,13 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, nwsx /= norm; nwsy /= norm; nwsz /= norm; - - /* JM does not remember why we should do this... - * if (nsx * nwnsx + nsy * nwnsy + nsz * nwnsz < 0.0) { - nwnsx = -nwnsx; - nwnsy = -nwnsy; - nwnsz = -nwnsz; + /* normal to ws interface boundary should point into fluid (same direction as gradient) */ + if (nwx * nwsx + nwy * nwsy + nwz * nwsz < 0.0) { + nwsx = -nwsx; + nwsy = -nwsy; + nwsz = -nwsz; } - */ + // common curve normal in the fluid surface tangent plane (rel. geodesic curvature) nwnx = twnsy * nwz - twnsz * nwy; nwny = twnsz * nwx - twnsx * nwz; @@ -4590,6 +4589,13 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, nwnx /= norm; nwny /= norm; nwnz /= norm; + /* normal to wn interface boundary should point into the solid */ + if (nsx * nwnx + nsy * nwny + nsz * nwnz > 0.0) { + nwnx = -nwnx; + nwny = -nwny; + nwnz = -nwnz; + } + if (length > 0.0) { // normal curvature component in the direction of the solid surface diff --git a/tests/lbpm_TwoPhase_analysis.cpp b/tests/lbpm_TwoPhase_analysis.cpp index e0d5ca05..902ef3cd 100644 --- a/tests/lbpm_TwoPhase_analysis.cpp +++ b/tests/lbpm_TwoPhase_analysis.cpp @@ -126,6 +126,16 @@ int main(int argc, char **argv) fillDouble.fill(Averages->SDs); if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); CalcDist(Averages->SDs,id,*Dm); + + /* move the solid surface by a voxel to improve contact angle measurement + for (k=0;kSDs(i,j,k) -= 1.0; + } + } + }*/ /* * * * * * * * * * * * * * * * * * * * * */ // Solve for the position of the fluid phase for (k=0;k