From 66fc8556088e68dda83bd4260ccf435c3379a2fe Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 19 Jan 2022 16:18:29 -0500 Subject: [PATCH] convention for inside / outside membrane link direction --- analysis/Membrane.cpp | 148 ++++++++++++++++++++++++++++++----------- analysis/Membrane.h | 1 + tests/TestMembrane.cpp | 6 +- 3 files changed, 116 insertions(+), 39 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index 4618380f..99e05234 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -552,9 +552,11 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra /* allocate memory */ membraneLinks = new int [2*mlink]; membraneDist = new double [2*mlink]; - + membraneTag = new int [mlink]; + /* construct the membrane*/ mlink = 0; + int insideMem = 0; int outsideMem = 0; for (k=1;k Dm, DoubleArray &Distance, IntArra neighbor=Map(i+1,j,k); dist=Distance(i+1,j,k); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 1*Np; - membraneLinks[2*mlink+1] = neighbor + 2*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 1*Np; + membraneLinks[outsideMem] = neighbor + 2*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i,j+1,k); dist=Distance(i,j+1,k); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 3*Np; - membraneLinks[2*mlink+1] = neighbor + 4*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 3*Np; + membraneLinks[outsideMem] = neighbor + 4*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i,j,k+1); dist=Distance(i,j,k+1); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 5*Np; - membraneLinks[2*mlink+1] = neighbor + 6*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 5*Np; + membraneLinks[outsideMem] = neighbor + 6*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i+1,j+1,k); dist=Distance(i+1,j+1,k); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 7*Np; - membraneLinks[2*mlink+1] = neighbor+8*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 7*Np; + membraneLinks[outsideMem] = neighbor+8*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i+1,j-1,k); dist=Distance(i+1,j-1,k); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 9*Np; - membraneLinks[2*mlink+1] = neighbor + 10*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 9*Np; + membraneLinks[outsideMem] = neighbor + 10*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i+1,j,k+1); dist=Distance(i+1,j,k+1); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 11*Np; - membraneLinks[2*mlink+1] = neighbor + 12*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 11*Np; + membraneLinks[outsideMem] = neighbor + 12*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i+1,j,k-1); dist=Distance(i+1,j,k-1); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 13*Np; - membraneLinks[2*mlink+1] = neighbor + 14*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 13*Np; + membraneLinks[outsideMem] = neighbor + 14*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i,j+1,k+1); dist=Distance(i,j+1,k+1); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 15*Np; - membraneLinks[2*mlink+1] =neighbor + 16*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 15*Np; + membraneLinks[outsideMem] =neighbor + 16*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i,j+1,k-1); dist=Distance(i,j+1,k-1); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 17*Np; - membraneLinks[2*mlink+1] = neighbor + 18*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 17*Np; + membraneLinks[outsideMem] = neighbor + 18*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } } diff --git a/analysis/Membrane.h b/analysis/Membrane.h index 3a842b75..563af5b0 100644 --- a/analysis/Membrane.h +++ b/analysis/Membrane.h @@ -26,6 +26,7 @@ public: int *neighborList; // modified neighborlist int *membraneLinks; // D3Q19 links that cross membrane + int *membraneTag; // label each link in the membrane double *membraneDist; // distance to membrane for each linked site /** diff --git a/tests/TestMembrane.cpp b/tests/TestMembrane.cpp index 319992b0..048be117 100644 --- a/tests/TestMembrane.cpp +++ b/tests/TestMembrane.cpp @@ -118,7 +118,6 @@ int main(int argc, char **argv) dist[q*Np + idx] = 0.0; } } - /* create a membrane data structure */ Membrane M(Dm, neighborList, Np); @@ -135,7 +134,7 @@ int main(int argc, char **argv) for (int mlink=0; mlink 0.f){ Dm->id[n] = 127; } + if (sum < 0.f){ + Dm->id[n] = 64; + } } } }