diff --git a/applications/solvers/molecularDynamics/mdFoam/Make/files b/applications/solvers/molecularDynamics/mdFoam/Make/files new file mode 100755 index 0000000000..6b867ec0a8 --- /dev/null +++ b/applications/solvers/molecularDynamics/mdFoam/Make/files @@ -0,0 +1,3 @@ +mdFoam.C + +EXE = $(FOAM_APPBIN)/mdFoam diff --git a/applications/solvers/molecularDynamics/mdFoam/Make/options b/applications/solvers/molecularDynamics/mdFoam/Make/options new file mode 100755 index 0000000000..14f2ed2f93 --- /dev/null +++ b/applications/solvers/molecularDynamics/mdFoam/Make/options @@ -0,0 +1,16 @@ +EXE_INC = \ + -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \ + -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \ + -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lmeshTools \ + -lfiniteVolume \ + -llagrangian \ + -lmolecule \ + -lpotential \ + -lmolecularMeasurements + diff --git a/applications/solvers/molecularDynamics/mdFoam/mdFoam.C b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C new file mode 100644 index 0000000000..e971bfaa52 --- /dev/null +++ b/applications/solvers/molecularDynamics/mdFoam/mdFoam.C @@ -0,0 +1,62 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Application + mdFoam + +Description + molecular dynamics solver for fluid dynamics + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +//#include "md.h" +#include "interactionLists.H" + +int main(int argc, char *argv[]) +{ + +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" + + interactionLists il(mesh, (1e-9*1e-9), false); + + Info<< labelListList(il.dil()) << endl; + + Info << "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + runTime++; + + Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info << "End\n" << endl; + + return(0); +} diff --git a/src/lagrangian/molecularDynamics/molecule/Make/files b/src/lagrangian/molecularDynamics/molecule/Make/files index dc12adad01..c5e58f9cc7 100755 --- a/src/lagrangian/molecularDynamics/molecule/Make/files +++ b/src/lagrangian/molecularDynamics/molecule/Make/files @@ -8,7 +8,7 @@ directInteractionList = $(interactionLists)/directInteractionList $(referralLists)/receivingReferralList.C $(referralLists)/sendingReferralList.C // $(referredCellList)/referredCellList.C -// $(referredCell)/referredCell.C +$(referredCell)/referredCell.C $(referredMolecule)/referredMolecule.C $(directInteractionList)/directInteractionList.C $(interactionLists)/interactionLists.C diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.C b/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.C index d4f26eafad..5c00ddd0fa 100644 --- a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.C +++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.C @@ -25,37 +25,29 @@ License \*---------------------------------------------------------------------------*/ #include "directInteractionList.H" - +#include "interactionLists.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // void Foam::directInteractionList::buildDirectInteractionList ( - scalar rCutMax, bool pointPointListBuild ) { Info<< nl << "Building list of direct interaction neighbours" << endl; - scalar rCutMaxSqr = rCutMax*rCutMax; - - const polyMesh& mesh(interactionLists_.mesh()); + const polyMesh& mesh(il_.mesh()); List > directInteractionList(mesh.nCells()); if (pointPointListBuild) { - Info<< "Point-Point direct interaction list build." << endl; + Info<< tab << "Point-Point direct interaction list build." << endl; label pointJIndex; forAll (mesh.points(), pointIIndex) { - const point& ptI - ( - mesh.points()[pointIIndex] - ); - for ( pointJIndex = pointIIndex; @@ -63,12 +55,7 @@ void Foam::directInteractionList::buildDirectInteractionList ++pointJIndex ) { - const point& ptJ - ( - mesh.points()[pointJIndex] - ); - - if (magSqr(ptI - ptJ) <= rCutMaxSqr) + if (il_.testPointPointDistance(pointIIndex, pointJIndex)) { const labelList& ptICells ( @@ -111,13 +98,13 @@ void Foam::directInteractionList::buildDirectInteractionList } else { - Info<< "Point-Face, Edge-Edge direct interaction list build." << endl; + Info<< tab << "Point-Face, Edge-Edge direct interaction list build." << endl; forAll (mesh.points(), p) { forAll(mesh.faces(), f) { - if(testPointFaceDistance(p, f)) + if(il_.testPointFaceDistance(p, f)) { const labelList& pCells(mesh.pointCells()[p]); @@ -153,11 +140,7 @@ void Foam::directInteractionList::buildDirectInteractionList if (cellN > cellI) { - if - ( - findIndex(directInteractionList[cellI], cellN) - == -1 - ) + if (findIndex(directInteractionList[cellI], cellN) == -1) { directInteractionList[cellI].append(cellN); } @@ -165,11 +148,7 @@ void Foam::directInteractionList::buildDirectInteractionList if (cellI > cellN) { - if - ( - findIndex(directInteractionList[cellN], cellI) - == -1 - ) + if (findIndex(directInteractionList[cellN], cellI) == -1) { directInteractionList[cellN].append(cellI); } @@ -195,7 +174,7 @@ void Foam::directInteractionList::buildDirectInteractionList { const edge& eJ(mesh.edges()[edgeJIndex]); - if (testEdgeEdgeDistance(eI, eJ)) + if (il_.testEdgeEdgeDistance(eI, eJ)) { const labelList& eICells(mesh.edgeCells()[edgeIIndex]); @@ -211,11 +190,7 @@ void Foam::directInteractionList::buildDirectInteractionList if (cellJ > cellI) { - if - ( - findIndex(directInteractionList[cellI], cellJ) - == -1 - ) + if (findIndex(directInteractionList[cellI], cellJ) == -1) { directInteractionList[cellI].append(cellJ); } @@ -223,11 +198,7 @@ void Foam::directInteractionList::buildDirectInteractionList if (cellI > cellJ) { - if - ( - findIndex(directInteractionList[cellJ], cellI) - == -1 - ) + if (findIndex(directInteractionList[cellJ], cellI) == -1) { directInteractionList[cellJ].append(cellI); } @@ -246,6 +217,13 @@ void Foam::directInteractionList::buildDirectInteractionList directInteractionList[transDIL].shrink() ); } + + // sorting DILs + + forAll((*this), dIL) + { + sort((*this)[dIL]); + } } @@ -253,38 +231,29 @@ void Foam::directInteractionList::buildDirectInteractionList Foam::directInteractionList::directInteractionList ( - const interactionLists& interactionLists, - scalar rCutMax, + const interactionLists& il, bool pointPointListBuild ) : - labelListList(interactionLists.mesh().nCells()), - interactionLists_(interactionLists) + labelListList(il.mesh().nCells()), + il_(il) { - buildDirectInteractionList(rCutMax, pointPointListBuild); + buildDirectInteractionList(pointPointListBuild); } Foam::directInteractionList::directInteractionList ( - const interactionLists& interactionLists + const interactionLists& il ) : - labelListList(interactionLists.mesh().nCells()), - interactionLists_(interactionLists) + labelListList(il.mesh().nCells()), + il_(il) { Info<< "Read directInteractionList from disk not implemented" << endl; } -// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * // - -Foam::autoPtr Foam::directInteractionList::New() -{ - return autoPtr(new directInteractionList); -} - - // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // Foam::directInteractionList::~directInteractionList() diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.H index 88dfe594e4..520670ea9b 100644 --- a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.H +++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.H @@ -57,13 +57,12 @@ class directInteractionList { // Private data - const interactionLists& interactionLists_; + const interactionLists& il_; // Private Member Functions void buildDirectInteractionList ( - scalar rCutMax, bool pointPointListBuild ); @@ -81,8 +80,7 @@ public: directInteractionList ( const interactionLists& interactionLists, - scalar rCutMax, - bool pointPointListBuild = false + bool pointPointListBuild ); //- Construct from file @@ -100,7 +98,7 @@ public: // Access - inline const interactionLists& interactionLists() const; + inline const interactionLists& il() const; // Check diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionListI.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionListI.H index 0acf54431a..176a2468b5 100644 --- a/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionListI.H +++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionListI.H @@ -29,9 +29,9 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // inline const Foam::interactionLists& - Foam::interactionLists::interactionLists() const + Foam::directInteractionList::il() const { - return interactionLists_; + return il_; } diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.C b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.C index 257bcf2aec..547630ef25 100644 --- a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.C +++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.C @@ -26,20 +26,24 @@ License #include "interactionLists.H" +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +Foam::scalar Foam::interactionLists::transTol = 1e-12; + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::interactionLists::interactionLists -interactionLists ( const polyMesh& mesh, - scalar rCutMax, + scalar rCutMaxSqr, bool pointPointListBuild ) : mesh_(mesh), - dil_(*this, rCutMax, pointPointListBuild) + rCutMaxSqr_(rCutMaxSqr), + dil_(*this, pointPointListBuild) {} @@ -58,5 +62,277 @@ Foam::interactionLists::~interactionLists() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::interactionLists::testPointPointDistance +( + const label ptI, + const label ptJ +) const +{ + return (magSqr(mesh_.points()[ptI] - mesh_.points()[ptJ]) <= rCutMaxSqr_); +} + + +bool Foam::interactionLists::testEdgeEdgeDistance +( + const edge& eI, + const edge& eJ +) const +{ + const vector& eJs(mesh_.points()[eJ.start()]); + const vector& eJe(mesh_.points()[eJ.end()]); + + return testEdgeEdgeDistance(eI, eJs, eJe); +} + + +bool Foam::interactionLists::testPointFaceDistance +( + const label p, + const label faceNo +) const +{ + const vector& pointPosition(mesh_.points()[p]); + + return testPointFaceDistance(pointPosition, faceNo); +} + + +bool Foam::interactionLists::testPointFaceDistance +( + const label p, + const referredCell& refCell +) const +{ + const vector& pointPosition(mesh_.points()[p]); + + forAll (refCell.faces(), rCF) + { + if + ( + testPointFaceDistance + ( + pointPosition, + refCell.faces()[rCF], + refCell.vertexPositions(), + refCell.faceCentres()[rCF], + refCell.faceAreas()[rCF] + ) + ) + { + return true; + } + } + + return false; +} + + +bool Foam::interactionLists::testPointFaceDistance +( + const vectorList& pointsToTest, + const label faceNo +) const +{ + forAll(pointsToTest, pTT) + { + const vector& p(pointsToTest[pTT]); + + // if any point in the list is in range of the face + // then the rest do not need to be tested and + // true can be returned + + if (testPointFaceDistance(p, faceNo)) + { + return true; + } + } + + return false; +} + + +bool Foam::interactionLists::testPointFaceDistance +( + const vector& p, + const label faceNo +) const +{ + const face& faceToTest(mesh_.faces()[faceNo]); + + const vector& faceC(mesh_.faceCentres()[faceNo]); + + const vector& faceA(mesh_.faceAreas()[faceNo]); + + const vectorList& points(mesh_.points()); + + return testPointFaceDistance + ( + p, + faceToTest, + points, + faceC, + faceA + ); +} + +bool Foam::interactionLists::testPointFaceDistance +( + const vector& p, + const labelList& faceToTest, + const vectorList& points, + const vector& faceC, + const vector& faceA +) const +{ + vector faceN(faceA/mag(faceA)); + + scalar perpDist((p - faceC) & faceN); + + if (magSqr(perpDist) > rCutMaxSqr_) + { + return false; + } + + vector pointOnPlane = (p - faceN * perpDist); + + if (magSqr(faceC - pointOnPlane) < rCutMaxSqr_*1e-8) + { + // If pointOnPlane is very close to the face centre + // then defining the local axes will be inaccurate + // and it is very likely that pointOnPlane will be + // inside the face, so return true if the points + // are in range to be safe + + return (magSqr(pointOnPlane - p) <= rCutMaxSqr_); + } + + vector xAxis = (faceC - pointOnPlane)/mag(faceC - pointOnPlane); + + vector yAxis = + ((faceC - pointOnPlane) ^ faceN) + /mag((faceC - pointOnPlane) ^ faceN); + + List local2DVertices(faceToTest.size()); + + forAll(faceToTest, fTT) + { + const vector& V(points[faceToTest[fTT]]); + + if (magSqr(V-p) <= rCutMaxSqr_) + { + return true; + } + + local2DVertices[fTT] = vector2D + ( + ((V - pointOnPlane) & xAxis), + ((V - pointOnPlane) & yAxis) + ); + } + + scalar localFaceCx((faceC - pointOnPlane) & xAxis); + + scalar la_valid = -1; + + forAll(local2DVertices, fV) + { + const vector2D& va(local2DVertices[fV]); + + const vector2D& vb + ( + local2DVertices[(fV + 1) % local2DVertices.size()] + ); + + if (mag(vb.y()-va.y()) > SMALL) + { + scalar la = + ( + va.x() - va.y()*((vb.x() - va.x())/(vb.y() - va.y())) + ) + /localFaceCx; + + scalar lv = -va.y()/(vb.y() - va.y()); + + + if (la >= 0 && la <= 1 && lv >= 0 && lv <= 1) + { + la_valid = la; + + break; + } + } + } + + if (la_valid < 0) + { + // perpendicular point inside face, nearest point is pointOnPlane; + return (magSqr(pointOnPlane-p) <= rCutMaxSqr_); + } + else + { + // perpendicular point outside face, nearest point is + // on edge that generated la_valid; + return + ( + magSqr(pointOnPlane + la_valid*(faceC - pointOnPlane) - p) + <= rCutMaxSqr_ + ); + } + + // if the algorithm hasn't returned anything by now then something has + // gone wrong. + + FatalErrorIn("interactionLists.C") << nl + << "point " << p << " to face " << faceToTest + << " comparison did not find a nearest point" + << " to be inside or outside face." + << abort(FatalError); + + return false; +} + + +bool Foam::interactionLists::testEdgeEdgeDistance +( + const edge& eI, + const vector& eJs, + const vector& eJe +) const +{ + vector a(eI.vec(mesh_.points())); + vector b(eJe - eJs); + + const vector& eIs(mesh_.points()[eI.start()]); + + vector c(eJs - eIs); + + vector crossab = a ^ b; + scalar magCrossSqr = magSqr(crossab); + + if (magCrossSqr > VSMALL) + { + // If the edges are parallel then a point-face + // search will pick them up + + scalar s = ((c ^ b) & crossab)/magCrossSqr; + scalar t = ((c ^ a) & crossab)/magCrossSqr; + + // Check for end points outside of range 0..1 + // If the closest point is outside this range + // a point-face search will have found it. + + return + ( + s >= 0 + && s <= 1 + && t >= 0 + && t <= 1 + && magSqr(eIs + a*s - eJs - b*t) <= rCutMaxSqr_ + ); + } + + return false; +} + // ************************************************************************* // diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H index 08f7d4d2a3..6d95e8b0d2 100644 --- a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H +++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H @@ -39,6 +39,7 @@ SourceFiles #include "polyMesh.H" #include "directInteractionList.H" +#include "referredCell.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -55,8 +56,24 @@ class interactionLists const polyMesh& mesh_; + scalar rCutMaxSqr_; + directInteractionList dil_; + // referredCellList ril_; + + // labelList realCellsWithinRCutMaxOfAnyReferringPatch_; + + // labelList realFacesWithinRCutMaxOfAnyReferringPatch_; + + // labelList realEdgesWithinRCutMaxOfAnyReferringPatch_; + + // labelList realPointsWithinRCutMaxOfAnyReferringPatch_; + + // List cellSendingReferralLists_; + + // List cellReceivingReferralLists_; + // Private Member Functions //- Disallow default bitwise copy construct @@ -68,13 +85,18 @@ class interactionLists public: + // Static data members + + //- Tolerance for checking that faces on a patch segment + static scalar transTol; + // Constructors //- Construct and create all information from the mesh interactionLists ( const polyMesh& mesh, - scalar rCutMax, + scalar rCutMaxSqr, bool pointPointListBuild = false ); @@ -88,11 +110,64 @@ public: // Member Functions + bool testPointPointDistance + ( + const label ptI, + const label ptJ + ) const; + + bool testPointFaceDistance + ( + const label p, + const label faceNo + ) const; + + bool testPointFaceDistance + ( + const label p, + const referredCell& refCell + ) const; + + bool testPointFaceDistance + ( + const vectorList& pointsToTest, + const label faceNo + ) const; + + bool testPointFaceDistance + ( + const vector& p, + const label faceNo + ) const; + + bool testPointFaceDistance + ( + const vector& p, + const labelList& faceToTest, + const vectorList& points, + const vector& faceC, + const vector& faceA + ) const; + + bool testEdgeEdgeDistance + ( + const edge& eI, + const edge& eJ + ) const; + + bool testEdgeEdgeDistance + ( + const edge& eI, + const vector& eJs, + const vector& eJe + ) const; + + // Access inline const polyMesh& mesh() const; - inline const directInteractionList& dil(); + inline const directInteractionList& dil() const; }; diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H index b9a2e326a4..ea8264632f 100644 --- a/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H +++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H @@ -26,31 +26,18 @@ License // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // +const Foam::polyMesh& Foam::interactionLists::mesh() const +{ + return mesh_; +} -// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // +const Foam::directInteractionList& Foam::interactionLists::dil() const +{ + return dil_; +} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.C b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.C index 98922d88ac..d971ec1487 100644 --- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.C +++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.C @@ -25,7 +25,7 @@ License \*----------------------------------------------------------------------------*/ #include "referredCell.H" -#include "moleculeCloud.H" +#include "interactionLists.H" namespace Foam { @@ -378,7 +378,7 @@ bool referredCell::duplicate(const referredCell& refCellDupl) const ( sourceProc_ == refCellDupl.sourceProc() && sourceCell_ == refCellDupl.sourceCell() - && mag(offset_ - refCellDupl.offset()) < moleculeCloud::transTol + && mag(offset_ - refCellDupl.offset()) < interactionLists::transTol ); } @@ -389,7 +389,7 @@ bool referredCell::duplicate(const label procNo,const label nCells) const ( sourceProc_ == procNo && sourceCell_ < nCells - && mag(offset_) < moleculeCloud::transTol + && mag(offset_) < interactionLists::transTol ); } diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/md.H b/src/lagrangian/molecularDynamics/molecule/mdTools/md.H index 5c018694bd..39ddd97b38 100644 --- a/src/lagrangian/molecularDynamics/molecule/mdTools/md.H +++ b/src/lagrangian/molecularDynamics/molecule/mdTools/md.H @@ -2,7 +2,6 @@ #define md_H # include "moleculeCloud.H" - # include "correlationFunction.H" # include "distribution.H" # include "reducedUnits.H"