From 413a129dc7fb9b936dfe982ec2515302c20bfc7d Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Mon, 5 Aug 2019 10:13:44 -0400 Subject: [PATCH 1/3] separate pressure from interface in subphase --- analysis/SubPhase.cpp | 42 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 34 insertions(+), 8 deletions(-) diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 318c0df2..d993ced3 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -492,6 +492,10 @@ void SubPhase::Full(){ double vol_wc_bulk = 0.0; double vol_nd_bulk = 0.0; double vol_wd_bulk = 0.0; + double count_wc = 0.0; + double count_nc = 0.0; + double count_wd = 0.0; + double count_nd = 0.0; for (k=kmin; k 0.0){ + if (morph_n->label(i,j,k) > 0 ){ + count_nd += 1.0; + nd.p += Pressure(n); + } + else{ + count_nc += 1.0; + nc.p += Pressure(n); + } + } + else{ + // water region + if (morph_w->label(i,j,k) > 0 ){ + count_wd += 1.0; + wd.p += Pressure(n); + } + else{ + count_wc += 1.0; + wc.p += Pressure(n); + } + } + if ( phi > 0.0){ if (morph_n->label(i,j,k) > 0 ){ vol_nd_bulk += 1.0; nd.M += nA*rho_n; @@ -522,7 +547,6 @@ void SubPhase::Full(){ nd.Py += nA*rho_n*uy; nd.Pz += nA*rho_n*uz; nd.K += nA*rho_n*(ux*ux + uy*uy + uz*uz); - nd.p += Pressure(n); } else{ vol_nc_bulk += 1.0; @@ -531,7 +555,6 @@ void SubPhase::Full(){ nc.Py += nA*rho_n*uy; nc.Pz += nA*rho_n*uz; nc.K += nA*rho_n*(ux*ux + uy*uy + uz*uz); - nc.p += Pressure(n); } } else{ @@ -543,7 +566,6 @@ void SubPhase::Full(){ wd.Py += nB*rho_w*uy; wd.Pz += nB*rho_w*uz; wd.K += nB*rho_w*(ux*ux + uy*uy + uz*uz); - wd.p += Pressure(n); } else{ vol_wc_bulk += 1.0; @@ -552,40 +574,44 @@ void SubPhase::Full(){ wc.Py += nB*rho_w*uy; wc.Pz += nB*rho_w*uz; wc.K += nB*rho_w*(ux*ux + uy*uy + uz*uz); - wc.p += Pressure(n); } } } } } } + count_wc=sumReduce( Dm->Comm, count_wc); + count_nc=sumReduce( Dm->Comm, count_nc); + count_wd=sumReduce( Dm->Comm, count_wd); + count_nd=sumReduce( Dm->Comm, count_nd); + gnd.p=sumReduce( Dm->Comm, nd.p) / count_nd; + gwd.p=sumReduce( Dm->Comm, wd.p) / count_wd; + gnc.p=sumReduce( Dm->Comm, nc.p) / count_nc; + gwc.p=sumReduce( Dm->Comm, wc.p) / count_wc; + gnd.M=sumReduce( Dm->Comm, nd.M); gnd.Px=sumReduce( Dm->Comm, nd.Px); gnd.Py=sumReduce( Dm->Comm, nd.Py); gnd.Pz=sumReduce( Dm->Comm, nd.Pz); gnd.K=sumReduce( Dm->Comm, nd.K); - gnd.p=sumReduce( Dm->Comm, nd.p); gwd.M=sumReduce( Dm->Comm, wd.M); gwd.Px=sumReduce( Dm->Comm, wd.Px); gwd.Py=sumReduce( Dm->Comm, wd.Py); gwd.Pz=sumReduce( Dm->Comm, wd.Pz); gwd.K=sumReduce( Dm->Comm, wd.K); - gwd.p=sumReduce( Dm->Comm, wd.p); gnc.M=sumReduce( Dm->Comm, nc.M); gnc.Px=sumReduce( Dm->Comm, nc.Px); gnc.Py=sumReduce( Dm->Comm, nc.Py); gnc.Pz=sumReduce( Dm->Comm, nc.Pz); gnc.K=sumReduce( Dm->Comm, nc.K); - gnc.p=sumReduce( Dm->Comm, nc.p); gwc.M=sumReduce( Dm->Comm, wc.M); gwc.Px=sumReduce( Dm->Comm, wc.Px); gwc.Py=sumReduce( Dm->Comm, wc.Py); gwc.Pz=sumReduce( Dm->Comm, wc.Pz); gwc.K=sumReduce( Dm->Comm, wc.K); - gwc.p=sumReduce( Dm->Comm, wc.p); giwn.Mn=sumReduce( Dm->Comm, iwn.Mn); giwn.Pnx=sumReduce( Dm->Comm, iwn.Pnx); From 9514cffa2f389d17e5667ec3c777fbb302bf8007 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Mon, 5 Aug 2019 14:51:38 -0400 Subject: [PATCH 2/3] update dcel dotprod checks --- analysis/dcel.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/analysis/dcel.cpp b/analysis/dcel.cpp index 8267f09f..06a92216 100644 --- a/analysis/dcel.cpp +++ b/analysis/dcel.cpp @@ -312,6 +312,7 @@ double DECL::EdgeAngle(int edge) } if (dotprod > 1.f) dotprod=1.f; + if (dotprod < -1.f) dotprod=-1.f; angle = acos(dotprod); /* project onto plane of cube face also works W = U - dotprod*V; From e1a0194a69d20bb2adbc202233fa142908547dba Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 6 Aug 2019 10:28:00 -0400 Subject: [PATCH 3/3] added nan check to edge angle --- analysis/dcel.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/analysis/dcel.cpp b/analysis/dcel.cpp index 06a92216..4c7be292 100644 --- a/analysis/dcel.cpp +++ b/analysis/dcel.cpp @@ -345,6 +345,7 @@ double DECL::EdgeAngle(int edge) // concave angle = -angle; } + if (angle != angle) angle = 0.0; //printf("angle=%f,dot=%f (Edge=%i, twin=%i): P={%f, %f, %f}, Q={%f, %f, %f} U={%f, %f, %f}, V={%f, %f, %f}\n",angle,dotprod,edge,halfedge.twin(edge),P.x,P.y,P.z,Q.x,Q.y,Q.z,U.x,U.y,U.z,V.x,V.y,V.z); return angle; }