ENH: molecularDynamics now using new InteractionLists.

This commit is contained in:
graham
2010-04-29 20:14:54 +01:00
parent 3e19b35b22
commit 23b5edd02d
28 changed files with 376 additions and 4829 deletions

View File

@ -39,6 +39,20 @@ int main(int argc, char *argv[])
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
Info<< "\nReading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
potential pot(mesh); potential pot(mesh);
moleculeCloud molecules(mesh, pot); moleculeCloud molecules(mesh, pot);

View File

@ -39,6 +39,20 @@ int main(int argc, char *argv[])
# include "createTime.H" # include "createTime.H"
# include "createMesh.H" # include "createMesh.H"
Info<< "\nReading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
potential pot(mesh); potential pot(mesh);
moleculeCloud molecules(mesh, pot); moleculeCloud molecules(mesh, pot);

View File

@ -1,24 +1,9 @@
interactionLists = interactionLists reducedUnits/reducedUnits.C
referredMolecule = $(interactionLists)/referredMolecule reducedUnits/reducedUnitsIO.C
referredCellList = $(interactionLists)/referredCellList
referredCell = $(interactionLists)/referredCell
directInteractionList = $(interactionLists)/directInteractionList
$(referredCellList)/referredCellList.C molecule/molecule.C
$(referredCell)/referredCell.C molecule/moleculeIO.C
$(referredMolecule)/referredMolecule.C
$(directInteractionList)/directInteractionList.C
$(interactionLists)/interactionLists.C
reducedUnits = reducedUnits moleculeCloud/moleculeCloud.C
$(reducedUnits)/reducedUnits.C
$(reducedUnits)/reducedUnitsIO.C
molecule = molecule
$(molecule)/molecule.C
$(molecule)/moleculeIO.C
moleculeCloud = moleculeCloud
$(moleculeCloud)/moleculeCloud.C
LIB = $(FOAM_LIBBIN)/libmolecule LIB = $(FOAM_LIBBIN)/libmolecule

View File

@ -1,11 +1,13 @@
EXE_INC = \ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \ -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \ -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude
EXE_LIBS = \ EXE_LIBS = \
-lfiniteVolume \ -lfiniteVolume \
-lmeshTools \
-llagrangian \ -llagrangian \
-lpotential \ -lpotential \
-lmolecularMeasurements -lmolecularMeasurements

View File

@ -1,346 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "directInteractionList.H"
#include "interactionLists.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::directInteractionList::buildDirectInteractionList
(
bool pointPointListBuild
)
{
Info<< nl << "Building list of direct interaction neighbours" << endl;
const polyMesh& mesh(il_.mesh());
List<DynamicList<label> > directInteractionList(mesh.nCells());
if (pointPointListBuild)
{
Info<< tab << "Point-Point direct interaction list build." << endl;
label pointJIndex;
forAll(mesh.points(), pointIIndex)
{
for
(
pointJIndex = pointIIndex;
pointJIndex != mesh.points().size();
++pointJIndex
)
{
if (il_.testPointPointDistance(pointIIndex, pointJIndex))
{
const labelList& ptICells
(
mesh.pointCells()[pointIIndex]
);
const labelList& ptJCells
(
mesh.pointCells()[pointJIndex]
);
forAll(ptICells, pIC)
{
const label cellI(ptICells[pIC]);
forAll(ptJCells, pJC)
{
const label cellJ(ptJCells[pJC]);
if (cellJ > cellI)
{
if
(
findIndex
(
directInteractionList[cellI],
cellJ
)
== -1
)
{
directInteractionList[cellI].append(cellJ);
}
}
if (cellI > cellJ)
{
if
(
findIndex
(
directInteractionList[cellJ],
cellI
)
==
-1
)
{
directInteractionList[cellJ].append(cellI);
}
}
}
}
}
}
}
}
else
{
Info<< tab << "Point-Face, Edge-Edge direct interaction list build."
<< endl;
forAll(mesh.points(), p)
{
forAll(mesh.faces(), f)
{
if (il_.testPointFaceDistance(p, f))
{
const labelList& pCells(mesh.pointCells()[p]);
const label cellO(mesh.faceOwner()[f]);
forAll(pCells, pC)
{
const label cellI(pCells[pC]);
// cells are not added to their own DIL
if (cellO > cellI)
{
if
(
findIndex
(
directInteractionList[cellI],
cellO
)
==
-1
)
{
directInteractionList[cellI].append(cellO);
}
}
if (cellI > cellO)
{
if
(
findIndex
(
directInteractionList[cellO],
cellI
)
==
-1
)
{
directInteractionList[cellO].append(cellI);
}
}
if (mesh.isInternalFace(f))
{
// boundary faces will not have neighbour
// information
const label cellN(mesh.faceNeighbour()[f]);
if (cellN > cellI)
{
if
(
findIndex
(
directInteractionList[cellI],
cellN
)
==
-1
)
{
directInteractionList[cellI].append(cellN);
}
}
if (cellI > cellN)
{
if
(
findIndex
(
directInteractionList[cellN],
cellI
)
==
-1
)
{
directInteractionList[cellN].append(cellI);
}
}
}
}
}
}
}
label edgeJIndex;
forAll(mesh.edges(), edgeIIndex)
{
const edge& eI(mesh.edges()[edgeIIndex]);
for
(
edgeJIndex = edgeIIndex + 1;
edgeJIndex != mesh.edges().size();
++edgeJIndex
)
{
const edge& eJ(mesh.edges()[edgeJIndex]);
if (il_.testEdgeEdgeDistance(eI, eJ))
{
const labelList& eICells(mesh.edgeCells()[edgeIIndex]);
const labelList& eJCells(mesh.edgeCells()[edgeJIndex]);
forAll(eICells, eIC)
{
const label cellI(eICells[eIC]);
forAll(eJCells, eJC)
{
const label cellJ(eJCells[eJC]);
if (cellJ > cellI)
{
if
(
findIndex
(
directInteractionList[cellI],
cellJ
)
==
-1
)
{
directInteractionList[cellI].append(cellJ);
}
}
if (cellI > cellJ)
{
if
(
findIndex
(
directInteractionList[cellJ],
cellI
)
==
-1
)
{
directInteractionList[cellJ].append(cellI);
}
}
}
}
}
}
}
}
forAll(directInteractionList, transDIL)
{
(*this)[transDIL].transfer
(
directInteractionList[transDIL].shrink()
);
}
// sorting DILs
forAll((*this), dIL)
{
sort((*this)[dIL]);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::directInteractionList::directInteractionList
(
const interactionLists& il,
bool pointPointListBuild
)
:
labelListList(il.mesh().nCells()),
il_(il)
{
if ((*this).size() > 1)
{
buildDirectInteractionList(pointPointListBuild);
}
else if ((*this).size() == 1)
{
Info<< nl
<< "Single cell mesh, no direct interaction lists required."
<< endl;
(*this)[0].setSize(0);
}
}
Foam::directInteractionList::directInteractionList
(
const interactionLists& il
)
:
labelListList(il.mesh().nCells()),
il_(il)
{
Info<< "Read directInteractionList from disk not implemented" << endl;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::directInteractionList::~directInteractionList()
{}
// ************************************************************************* //

View File

@ -1,124 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::directInteractionList
Description
SourceFiles
directInteractionListI.H
directInteractionList.C
\*---------------------------------------------------------------------------*/
#ifndef directInteractionList_H
#define directInteractionList_H
#include "polyMesh.H"
#include "List.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class interactionLists;
/*---------------------------------------------------------------------------*\
Class directInteractionList Declaration
\*---------------------------------------------------------------------------*/
class directInteractionList
:
public labelListList
{
// Private data
const interactionLists& il_;
// Private Member Functions
void buildDirectInteractionList
(
bool pointPointListBuild
);
//- Disallow default bitwise copy construct
directInteractionList(const directInteractionList&);
//- Disallow default bitwise assignment
void operator=(const directInteractionList&);
public:
// Constructors
//- Construct lists by searching the mesh
directInteractionList
(
const interactionLists& il,
bool pointPointListBuild
);
//- Construct from file
directInteractionList
(
const interactionLists& il
);
//- Destructor
~directInteractionList();
// Member Functions
// Access
inline const interactionLists& il() const;
// IOstream Operators
friend Istream& operator>>(Istream&, directInteractionList&);
friend Ostream& operator<<(Ostream&, const directInteractionList&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "directInteractionListI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,34 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::interactionLists& Foam::directInteractionList::il() const
{
return il_;
}
// ************************************************************************* //

View File

@ -1,691 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "interactionLists.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
Foam::scalar Foam::interactionLists::transTol = 1e-12;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::interactionLists::buildCellReferralLists()
{
Info<< nl << "Determining molecule referring schedule" << endl;
const referredCellList& refIntL(ril());
DynamicList<label> referralProcs;
// Run through all referredCells to build list of interacting processors
forAll(refIntL, rIL)
{
const referredCell& rC(refIntL[rIL]);
if (findIndex(referralProcs, rC.sourceProc()) == -1)
{
referralProcs.append(rC.sourceProc());
}
}
referralProcs.shrink();
// Pout << "referralProcs: " << nl << referralProcs << endl;
List<DynamicList<label> > cellSendingReferralLists(referralProcs.size());
List<DynamicList<DynamicList<label> > >
cellReceivingReferralLists(referralProcs.size());
// Run through all referredCells again building up send and receive info
forAll(refIntL, rIL)
{
const referredCell& rC(refIntL[rIL]);
label rPI = findIndex(referralProcs, rC.sourceProc());
DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
label existingSource = findIndex(sRL, rC.sourceCell());
// Check to see if this source cell has already been allocated to
// come to this processor. If not, add the source cell to the sending
// list and add the current referred cell to the receiving list.
// It shouldn't be possible for the sending and receiving lists to be
// different lengths, because their append operations happen at the
// same time.
if (existingSource == -1)
{
sRL.append(rC.sourceCell());
rRL.append
(
DynamicList<label> (labelList(1,rIL))
);
}
else
{
rRL[existingSource].append(rIL);
rRL[existingSource].shrink();
}
}
forAll(referralProcs, rPI)
{
DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
sRL.shrink();
rRL.shrink();
}
// It is assumed that cell exchange is reciprocal, if proc A has cells to
// send to proc B, then proc B must have some to send to proc A.
cellReceivingReferralLists_.setSize(referralProcs.size());
cellSendingReferralLists_.setSize(referralProcs.size());
forAll(referralProcs, rPI)
{
DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
labelListList translLL(rRL.size());
forAll(rRL, rRLI)
{
translLL[rRLI] = rRL[rRLI];
}
cellReceivingReferralLists_[rPI] = receivingReferralList
(
referralProcs[rPI],
translLL
);
}
// Send sendingReferralLists to each interacting processor.
forAll(referralProcs, rPI)
{
DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
if (referralProcs[rPI] != Pstream::myProcNo())
{
if (Pstream::parRun())
{
OPstream toInteractingProc
(
Pstream::blocking,
referralProcs[rPI]
);
toInteractingProc << sendingReferralList
(
Pstream::myProcNo(),
sRL
);
}
}
}
// Receive sendingReferralLists from each interacting processor.
forAll(referralProcs, rPI)
{
if (referralProcs[rPI] != Pstream::myProcNo())
{
if (Pstream::parRun())
{
IPstream fromInteractingProc
(
Pstream::blocking,
referralProcs[rPI]
);
fromInteractingProc >> cellSendingReferralLists_[rPI];
}
}
else
{
cellSendingReferralLists_[rPI] = sendingReferralList
(
Pstream::myProcNo(),
cellSendingReferralLists[rPI]
);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::interactionLists::interactionLists
(
const polyMesh& mesh,
scalar rCutMaxSqr,
bool pointPointListBuild
)
:
mesh_(mesh),
rCutMaxSqr_(rCutMaxSqr),
dil_(*this, pointPointListBuild),
ril_(*this, pointPointListBuild),
cellSendingReferralLists_(),
cellReceivingReferralLists_()
{
buildCellReferralLists();
}
Foam::interactionLists::interactionLists(const polyMesh& mesh)
:
mesh_(mesh),
dil_(*this),
ril_(*this)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
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<vector2D> 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::testPointFaceDistance"
"("
"const vector&, "
"const labelList&, "
"const vectorList&, "
"const vector&, "
"const vector&"
") const"
) << "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;
}
const Foam::labelList Foam::interactionLists::realCellsInRangeOfSegment
(
const labelList& segmentFaces,
const labelList& segmentEdges,
const labelList& segmentPoints
) const
{
DynamicList<label> realCellsFoundInRange;
forAll(segmentFaces, sF)
{
const label f = segmentFaces[sF];
forAll(mesh_.points(), p)
{
if (testPointFaceDistance(p, f))
{
const labelList& pCells(mesh_.pointCells()[p]);
forAll(pCells, pC)
{
const label cellI(pCells[pC]);
if (findIndex(realCellsFoundInRange, cellI) == -1)
{
realCellsFoundInRange.append(cellI);
}
}
}
}
}
forAll(segmentPoints, sP)
{
const label p = segmentPoints[sP];
forAll(mesh_.faces(), f)
{
if (testPointFaceDistance(p, f))
{
const label cellO(mesh_.faceOwner()[f]);
if (findIndex(realCellsFoundInRange, cellO) == -1)
{
realCellsFoundInRange.append(cellO);
}
if (mesh_.isInternalFace(f))
{
// boundary faces will not have neighbour information
const label cellN(mesh_.faceNeighbour()[f]);
if (findIndex(realCellsFoundInRange, cellN) == -1)
{
realCellsFoundInRange.append(cellN);
}
}
}
}
}
forAll(segmentEdges, sE)
{
const edge& eJ(mesh_.edges()[segmentEdges[sE]]);
forAll(mesh_.edges(), edgeIIndex)
{
const edge& eI(mesh_.edges()[edgeIIndex]);
if (testEdgeEdgeDistance(eI, eJ))
{
const labelList& eICells(mesh_.edgeCells()[edgeIIndex]);
forAll(eICells, eIC)
{
const label cellI(eICells[eIC]);
if (findIndex(realCellsFoundInRange, cellI) == -1)
{
realCellsFoundInRange.append(cellI);
}
}
}
}
}
return realCellsFoundInRange.shrink();
}
const Foam::labelList Foam::interactionLists::referredCellsInRangeOfSegment
(
const List<referredCell>& referredInteractionList,
const labelList& segmentFaces,
const labelList& segmentEdges,
const labelList& segmentPoints
) const
{
DynamicList<label> referredCellsFoundInRange;
forAll(segmentFaces, sF)
{
const label f = segmentFaces[sF];
forAll(referredInteractionList, rIL)
{
const vectorList& refCellPoints
= referredInteractionList[rIL].vertexPositions();
if (testPointFaceDistance(refCellPoints, f))
{
if (findIndex(referredCellsFoundInRange, rIL) == -1)
{
referredCellsFoundInRange.append(rIL);
}
}
}
}
forAll(segmentPoints, sP)
{
const label p = segmentPoints[sP];
forAll(referredInteractionList, rIL)
{
const referredCell& refCell(referredInteractionList[rIL]);
if (testPointFaceDistance(p, refCell))
{
if (findIndex(referredCellsFoundInRange, rIL) == -1)
{
referredCellsFoundInRange.append(rIL);
}
}
}
}
forAll(segmentEdges, sE)
{
const edge& eI(mesh_.edges()[segmentEdges[sE]]);
forAll(referredInteractionList, rIL)
{
const vectorList& refCellPoints
= referredInteractionList[rIL].vertexPositions();
const edgeList& refCellEdges
= referredInteractionList[rIL].edges();
forAll(refCellEdges, rCE)
{
const edge& eJ(refCellEdges[rCE]);
if
(
testEdgeEdgeDistance
(
eI,
refCellPoints[eJ.start()],
refCellPoints[eJ.end()]
)
)
{
if (findIndex(referredCellsFoundInRange, rIL) == -1)
{
referredCellsFoundInRange.append(rIL);
}
}
}
}
}
return referredCellsFoundInRange.shrink();
}
// ************************************************************************* //

View File

@ -1,214 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::interactionLists
Description
SourceFiles
interactionListsI.H
interactionLists.C
interactionListsIO.C
\*---------------------------------------------------------------------------*/
#ifndef interactionLists_H
#define interactionLists_H
#include "polyMesh.H"
#include "vector2D.H"
#include "directInteractionList.H"
#include "referredCell.H"
#include "referredCellList.H"
#include "sendingReferralList.H"
#include "receivingReferralList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class interactionLists Declaration
\*---------------------------------------------------------------------------*/
class interactionLists
{
// Private data
const polyMesh& mesh_;
scalar rCutMaxSqr_;
directInteractionList dil_;
referredCellList ril_;
List<sendingReferralList> cellSendingReferralLists_;
List<receivingReferralList> cellReceivingReferralLists_;
// Private Member Functions
//- Build referralLists which define how to send information
// to referredCells to source cells
void buildCellReferralLists();
//- Disallow default bitwise copy construct
interactionLists(const interactionLists&);
//- Disallow default bitwise assignment
void operator=(const 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 rCutMaxSqr,
bool pointPointListBuild = false
);
//- Construct from file
interactionLists(const polyMesh& mesh);
//- Destructor
~interactionLists();
// 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;
const labelList realCellsInRangeOfSegment
(
const labelList& segmentFaces,
const labelList& segmentEdges,
const labelList& segmentPoints
) const;
const labelList referredCellsInRangeOfSegment
(
const List<referredCell>& referredInteractionList,
const labelList& segmentFaces,
const labelList& segmentEdges,
const labelList& segmentPoints
) const;
// Access
inline const polyMesh& mesh() const;
inline const directInteractionList& dil() const;
inline const referredCellList& ril() const;
inline referredCellList& ril();
inline const List<sendingReferralList>&
cellSendingReferralLists() const;
inline const List<receivingReferralList>&
cellReceivingReferralLists() const;
inline label nInteractingProcs() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "interactionListsI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,77 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::polyMesh& Foam::interactionLists::mesh() const
{
return mesh_;
}
const Foam::directInteractionList& Foam::interactionLists::dil() const
{
return dil_;
}
inline const Foam::referredCellList& Foam::interactionLists::ril() const
{
return ril_;
}
inline Foam::referredCellList& Foam::interactionLists::ril()
{
return ril_;
}
inline const Foam::List<Foam::sendingReferralList>&
Foam::interactionLists::cellSendingReferralLists() const
{
return cellSendingReferralLists_;
}
inline const Foam::List<Foam::receivingReferralList>&
Foam::interactionLists::cellReceivingReferralLists() const
{
return cellReceivingReferralLists_;
}
inline Foam::label Foam::interactionLists::nInteractingProcs() const
{
return cellReceivingReferralLists_.size();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -1,448 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*----------------------------------------------------------------------------*/
#include "referredCell.H"
#include "interactionLists.H"
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void referredCell::setConstructionData
(
const polyMesh& mesh,
const label sourceCell
)
{
// Points
const labelList& points = mesh.cellPoints()[sourceCell];
vectorList sourceCellVertices(points.size());
forAll(sourceCellVertices, sCV)
{
sourceCellVertices[sCV] = mesh.points()[points[sCV]];
}
vertexPositions_ = referPositions(sourceCellVertices);
// Edges
const labelList& edges = mesh.cellEdges()[sourceCell];
edgeList sourceCellEdges(edges.size());
forAll(sourceCellEdges, sCE)
{
sourceCellEdges[sCE] = mesh.edges()[edges[sCE]];
}
locallyMapEdgeList(points, sourceCellEdges);
// Faces
labelList faces(mesh.cells()[sourceCell]);
vectorList sourceCellFaceCentres(faces.size());
vectorList sourceCellFaceAreas(faces.size());
labelListList sourceCellFaces(faces.size());
forAll(faces, f)
{
sourceCellFaces[f] = mesh.faces()[faces[f]];
sourceCellFaceCentres[f] = mesh.faceCentres()[faces[f]];
sourceCellFaceAreas[f] = mesh.faceAreas()[faces[f]];
}
locallyMapFaceList(points, sourceCellFaces);
faceCentres_ = referPositions(sourceCellFaceCentres);
faceAreas_ = rotateVectors(sourceCellFaceAreas);
}
void referredCell::locallyMapEdgeList
(
const labelList& points,
const edgeList& sourceCellEdges
)
{
edges_.setSize(sourceCellEdges.size());
forAll(sourceCellEdges, sCE)
{
const edge& e(sourceCellEdges[sCE]);
edges_[sCE].start() = findIndex(points, e.start());
edges_[sCE].end() = findIndex(points, e.end());
if
(
edges_[sCE].start() == -1
|| edges_[sCE].end() == -1
)
{
FatalErrorIn("Foam::referredCell::locallyMapEdgeList")
<< "edgeList and points labelList for "
<< "referred cell do not match: "
<< nl << "points: " << points
<< nl << "egdes: " << sourceCellEdges
<< abort(FatalError);
}
}
}
void referredCell::locallyMapFaceList
(
const labelList& points,
const labelListList& sourceCellFaces
)
{
faces_.setSize(sourceCellFaces.size());
forAll(sourceCellFaces, sCF)
{
const labelList& sourceCellFace(sourceCellFaces[sCF]);
labelList& localFace(faces_[sCF]);
localFace.setSize(sourceCellFace.size());
forAll(sourceCellFace, p)
{
localFace[p] = findIndex(points, sourceCellFace[p]);
if (localFace[p] == -1)
{
FatalErrorIn("Foam::referredCell::locallyMapEdgeList")
<< "edgeList and points labelList for "
<< "referred cell do not match: "
<< nl << "points: " << points
<< nl << "faces: " << sourceCellFaces
<< abort(FatalError);
}
}
}
}
vector referredCell::referPosition(const vector& positionToRefer)
{
return offset_ + (rotation_ & positionToRefer);
}
vectorList referredCell::referPositions(const vectorList& positionsToRefer)
{
return offset_ + (rotation_ & positionsToRefer);
}
vector referredCell::rotateVector(const vector& vectorToRotate)
{
return rotation_ & vectorToRotate;
}
vectorList referredCell::rotateVectors(const vectorList& vectorsToRotate)
{
return rotation_ & vectorsToRotate;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
referredCell::referredCell()
:
DynamicList<referredMolecule>(),
sourceProc_(-1),
sourceCell_(-1),
vertexPositions_(),
offset_(vector::zero),
rotation_(I)
{}
referredCell::referredCell
(
const polyMesh& mesh,
const label sourceProc,
const label sourceCell,
const vector& offset,
const tensor& rotation
)
:
DynamicList<referredMolecule>(),
sourceProc_(sourceProc),
sourceCell_(sourceCell),
offset_(offset),
rotation_(rotation)
{
setConstructionData(mesh, sourceCell);
}
referredCell::referredCell
(
const label sourceProc,
const label sourceCell,
const vectorList& vertexPositions,
const edgeList& localEdges,
const labelListList& localFaces,
const vectorList& faceCentres,
const vectorList& faceAreas,
const vector& offset,
const tensor& rotation
)
:
DynamicList<referredMolecule>(),
sourceProc_(sourceProc),
sourceCell_(sourceCell),
edges_(localEdges),
faces_(localFaces),
offset_(offset),
rotation_(rotation)
{
// Supplied vertexPositions, faceCentres and faceAreas are of the
// "original" cell, and need to be transformed to the referred
// locations on construction
vertexPositions_ = referPositions(vertexPositions);
faceCentres_ = referPositions(faceCentres);
faceAreas_ = rotateVectors(faceAreas);
}
referredCell::referredCell
(
const polyMesh& mesh,
const label sourceProc,
const label sourceCell,
const vector& cS,
const vector& cD,
const vector& nS,
const vector& nD
)
:
DynamicList<referredMolecule>(),
sourceProc_(sourceProc),
sourceCell_(sourceCell)
{
// It is assumed that the vectors originating from the faces being referred
// here are correct periodic faces - i.e. they have the same area etc.
vector nA = -nS/mag(nS);
vector nB = nD/mag(nD);
rotation_ = rotationTensor(nA, nB);
offset_ = cD - (rotation_ & cS);
// Allow sourceCell = -1 to create a dummy referredCell
// to obtain the transformation
if (sourceCell >= 0)
{
setConstructionData(mesh, sourceCell);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
referredCell::~referredCell()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
referredCell referredCell::reRefer
(
const vector& cS,
const vector& cD,
const vector& nS,
const vector& nD
)
{
vector nA = -nS/mag(nS);
vector nB = nD/mag(nD);
tensor newRotation = rotationTensor(nA, nB);
vector newOffset = cD - (newRotation & cS);
tensor reReferredRotation = newRotation & rotation_;
vector reReferredOffset = newOffset + (newRotation & offset_);
return referredCell
(
sourceProc_,
sourceCell_,
rotation_.T() & (vertexPositions_ - offset_),
edges_,
faces_,
rotation_.T() & (faceCentres_ - offset_),
rotation_.T() & (faceAreas_),
reReferredOffset,
reReferredRotation
);
}
vector referredCell::referPosition(const vector& positionToRefer) const
{
return offset_ + (rotation_ & positionToRefer);
}
vectorList referredCell::referPosition
(
const vectorList& positionsToRefer
) const
{
return offset_ + (rotation_ & positionsToRefer);
}
vector referredCell::rotateVector(const vector& vectorToRotate) const
{
return rotation_ & vectorToRotate;
}
vectorList referredCell::rotateVectors(const vectorList& vectorsToRotate) const
{
return rotation_ & vectorsToRotate;
}
void referredCell::referInMols(const List<referredMolecule>& incomingMols)
{
clear();
forAll(incomingMols, iM)
{
append
(
referredMolecule
(
incomingMols[iM].id(),
referPosition
(
incomingMols[iM].position()
),
referPosition
(
incomingMols[iM].sitePositions()
)
)
);
}
shrink();
}
bool referredCell::duplicate(const referredCell& refCellDupl) const
{
return
(
sourceProc_ == refCellDupl.sourceProc()
&& sourceCell_ == refCellDupl.sourceCell()
&& mag(offset_ - refCellDupl.offset()) < interactionLists::transTol
);
}
bool referredCell::duplicate(const label procNo,const label nCells) const
{
return
(
sourceProc_ == procNo
&& sourceCell_ < nCells
&& mag(offset_) < interactionLists::transTol
);
}
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
Istream& operator>>(Istream& is, referredCell& rC)
{
is >> rC.sourceProc_
>> rC.sourceCell_
>> rC.vertexPositions_
>> rC.edges_
>> rC.faces_
>> rC.faceCentres_
>> rC.faceAreas_
>> rC.offset_
>> rC.rotation_;
is.check("Istream& operator<<(Istream& f, const referredCell& rC");
return is;
}
Ostream& operator<<(Ostream& os, const referredCell& rC)
{
os << rC.sourceProc()
<< token::SPACE << rC.sourceCell()
<< token::SPACE << rC.vertexPositions()
<< token::SPACE << rC.edges()
<< token::SPACE << rC.faces()
<< token::SPACE << rC.faceCentres()
<< token::SPACE << rC.faceAreas()
<< token::SPACE << rC.offset()
<< token::SPACE << rC.rotation();
os.check("Ostream& operator<<(Ostream& f, const referredCell& rC");
return os;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -1,275 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::referredCell
Description
SourceFiles
referredCellI.H
referredCell.C
\*---------------------------------------------------------------------------*/
#ifndef referredCell_H
#define referredCell_H
#include "vector.H"
#include "vectorList.H"
#include "tensor.H"
#include "transform.H"
#include "DynamicList.H"
#include "labelList.H"
#include "edgeList.H"
#include "polyMesh.H"
#include "referredMolecule.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class referredCell Declaration
\*---------------------------------------------------------------------------*/
class referredCell
:
public DynamicList<referredMolecule>
{
// Private data
label sourceProc_;
label sourceCell_;
//- Referred vertex positions
vectorList vertexPositions_;
edgeList edges_;
labelListList faces_;
vectorList faceCentres_;
vectorList faceAreas_;
labelList realCellsForInteraction_;
vector offset_;
tensor rotation_;
// Private Member Functions
void setConstructionData
(
const polyMesh& mesh,
const label sourceCell
);
void locallyMapEdgeList
(
const labelList& points,
const edgeList& sourceCellEdges
);
void locallyMapFaceList
(
const labelList& points,
const labelListList& sourceCellFaces
);
vector referPosition(const vector& positionToRefer);
vectorList referPositions(const vectorList& positionsToRefer);
vector rotateVector(const vector& vectorToRotate);
vectorList rotateVectors(const vectorList& vectorsToRotate);
public:
// Constructors
//- Construct null
referredCell();
//- Construct from components with external edge information
referredCell
(
const polyMesh& mesh,
const label sourceProc,
const label sourceCell,
const vector& offset,
const tensor& rotation
);
//- Construct from components with existing local edge information
referredCell
(
const label sourceProc,
const label sourceCell,
const vectorList& vertexPositions,
const edgeList& localEdges,
const labelListList& localFaces,
const vectorList& faceCentres,
const vectorList& faceAreas,
const vector& offset,
const tensor& rotation
);
//- Construct from pair of face centers (c) and plain
// face normals (n) (no need to make unit vectors or
// reverse one direction)
// Order of vectors important (S = source, D = Destination).
// External edge information.
referredCell
(
const polyMesh& mesh,
const label sourceProc,
const label sourceCell,
const vector& cS,
const vector& cD,
const vector& nS,
const vector& nD
);
//- Destructor
virtual ~referredCell();
// Member Functions
//- Take this referredCell object that has already had its transform
// calculated and refer it on again, retaining same source info.
referredCell reRefer
(
const vector& cS,
const vector& cD,
const vector& nS,
const vector& nD
);
//- Use internal transformatation values to transform the given
// postion to its new location.
vector referPosition(const vector& positionToRefer) const;
//- Use internal transformatation values to transform the given
// list of postions to their new locations.
vectorList referPosition(const vectorList& positionsToRefer) const;
//- Use internal transformatation values to rotate the given vector
vector rotateVector(const vector& vectorToRotate) const;
//- Use internal transformatation values to rotate the given
// list of vectors
vectorList rotateVectors(const vectorList& vectorsToRotate) const;
//- referInMols clears, the stored list of referred mols and takes
// in a list of referred mols coming from a source processor,
// transformingtheir positions
void referInMols(const List<referredMolecule>& incomingMols);
//- duplicate() function to test whether a referred or real cell
// supplied by arguement is a duplicate of this referredCell.
// Can be used bi-directionally - i.e. can be called on an existing
// referred cell with a proposed referredCell as argument,
// or vice versa. Can only be called by a proposed referredCell with
// a real cell index as arguement to test to see if the proposed
// referredCell is a duplicate.
// A duplicate cell is defined as one which has the same source
// processor, source cell, and an equal offset. Real cells have zero
// offset by definition.
bool duplicate(const referredCell& refCellDupl) const;
bool duplicate(const label procNo, const label nCells) const;
// Access
inline label sourceProc() const;
inline label sourceCell() const;
inline const vector& offset() const;
inline const tensor& rotation() const;
inline const vectorList& vertexPositions() const;
inline const edgeList& edges() const;
inline const labelListList& faces() const;
inline const vectorList& faceCentres() const;
inline const vectorList& faceAreas() const;
inline labelList& realCells();
inline const labelList& realCellsForInteraction() const;
// Friend Operators
inline friend bool operator==
(
const referredCell& a,
const referredCell& b
);
inline friend bool operator!=
(
const referredCell& a,
const referredCell& b
);
// IOstream Operators
friend Istream& operator>>(Istream&, referredCell&);
friend Ostream& operator<<(Ostream&, const referredCell&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "referredCellI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,126 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline label referredCell::sourceProc() const
{
return sourceProc_;
}
inline label referredCell::sourceCell() const
{
return sourceCell_;
}
inline const vector& referredCell::offset() const
{
return offset_;
}
inline const tensor& referredCell::rotation() const
{
return rotation_;
}
inline const vectorList& referredCell::vertexPositions() const
{
return vertexPositions_;
}
inline const edgeList& referredCell::edges() const
{
return edges_;
}
inline const labelListList& referredCell::faces() const
{
return faces_;
}
inline const vectorList& referredCell::faceCentres() const
{
return faceCentres_;
}
inline const vectorList& referredCell::faceAreas() const
{
return faceAreas_;
}
inline labelList& referredCell::realCells()
{
return realCellsForInteraction_;
}
inline const labelList& referredCell::realCellsForInteraction() const
{
return realCellsForInteraction_;
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
inline bool operator==
(
const referredCell& a,
const referredCell& b
)
{
return const_cast<referredCell&>(a).duplicate
(
const_cast<const referredCell&>(b)
);
}
inline bool operator!=
(
const referredCell& a,
const referredCell& b
)
{
return !(a == b);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -1,109 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::referredCellList
Description
SourceFiles
referredCellListI.H
referredCellList.C
\*---------------------------------------------------------------------------*/
#ifndef referredCellList_H
#define referredCellList_H
#include "referredCell.H"
#include "molecule.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class interactionLists;
/*---------------------------------------------------------------------------*\
Class referredCellList Declaration
\*---------------------------------------------------------------------------*/
class referredCellList
:
public List<referredCell>
{
// Private data
const interactionLists& il_;
// Private Member Functions
void buildReferredCellList
(
bool pointPointListBuild
);
public:
// Constructors
//- Construct lists by searching the mesh
referredCellList
(
interactionLists& il,
bool pointPointListBuild
);
//- Construct from file
referredCellList (interactionLists& il);
//- Destructor
~referredCellList();
// Member Functions
void referMolecules(const List<DynamicList<molecule*> >& cellOccupancy);
inline const interactionLists& il() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "referredCellListI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,34 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline const Foam::interactionLists& Foam::referredCellList::il() const
{
return il_;
}
// ************************************************************************* //

View File

@ -1,85 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*----------------------------------------------------------------------------*/
#include "referredMolecule.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::referredMolecule::referredMolecule()
{}
Foam::referredMolecule::referredMolecule
(
const label id,
const vector& position,
const List<vector>& sitePositions
)
:
id_(id),
position_(position),
sitePositions_(sitePositions)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::referredMolecule::~referredMolecule()
{}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
Foam::Istream& Foam::operator>>
(
Istream& is,
referredMolecule& rM
)
{
is >> rM.id_ >> rM.position_ >> rM.sitePositions_;
is.check("Istream& operator<<(Istream& f, const referredMolecule& sRL");
return is;
}
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const referredMolecule& rM
)
{
os << rM.id()
<< token::SPACE << rM.position()
<< token::SPACE << rM.sitePositions();
os.check("Ostream& operator<<(Ostream& f, const referredMolecule& rM");
return os;
}
// ************************************************************************* //

View File

@ -1,135 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::referredMolecule
Description
SourceFiles
referredMoleculeI.H
referredMolecule.C
\*---------------------------------------------------------------------------*/
#ifndef referredMolecule_H
#define referredMolecule_H
#include "vector.H"
#include "List.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class referredMolecule Declaration
\*---------------------------------------------------------------------------*/
class referredMolecule
{
// Private data
label id_;
vector position_;
List<vector> sitePositions_;
public:
// Constructors
//- Construct null
referredMolecule();
//- Construct from components
referredMolecule
(
const label id,
const vector& position,
const List<vector>& sitePositions
);
//- Destructor
virtual ~referredMolecule();
// Member Functions
// Access
inline label id() const;
inline const vector& position() const;
inline const List<vector>& sitePositions() const;
// Friend Operators
inline friend bool operator==
(
const referredMolecule& a,
const referredMolecule& b
);
inline friend bool operator!=
(
const referredMolecule& a,
const referredMolecule& b
);
// IOstream Operators
friend Istream& operator>>
(
Istream&,
referredMolecule&
);
friend Ostream& operator<<
(
Ostream&,
const referredMolecule&
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "referredMoleculeI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,74 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline Foam::label Foam::referredMolecule::id() const
{
return id_;
}
inline const Foam::vector& Foam::referredMolecule::position() const
{
return position_;
}
inline const Foam::List<Foam::vector>&
Foam::referredMolecule::sitePositions() const
{
return sitePositions_;
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
inline bool Foam::operator==
(
const Foam::referredMolecule& a,
const Foam::referredMolecule& b
)
{
return
(
a.id() == b.id()
&& a.position() == b.position()
);
}
inline bool Foam::operator!=
(
const Foam::referredMolecule& a,
const Foam::referredMolecule& b
)
{
return !(a == b);
}
// ************************************************************************* //

View File

@ -209,8 +209,20 @@ bool Foam::molecule::move(molecule::trackData& td)
void Foam::molecule::transformProperties(const tensor& T) void Foam::molecule::transformProperties(const tensor& T)
{ {
Particle<molecule>::transformProperties(T);
Q_ = T & Q_; Q_ = T & Q_;
v_ = transform(T, v_);
a_ = transform(T, a_);
pi_ = Q_.T() & transform(T, Q_ & pi_);
tau_ = Q_.T() & transform(T, Q_ & tau_);
rf_ = transform(T, rf_);
sitePositions_ = position_ + (T & (sitePositions_ - position_)); sitePositions_ = position_ + (T & (sitePositions_ - position_));
siteForces_ = T & siteForces_; siteForces_ = T & siteForces_;
@ -219,10 +231,14 @@ void Foam::molecule::transformProperties(const tensor& T)
void Foam::molecule::transformProperties(const vector& separation) void Foam::molecule::transformProperties(const vector& separation)
{ {
Particle<molecule>::transformProperties(separation);
if (special_ == SPECIAL_TETHERED) if (special_ == SPECIAL_TETHERED)
{ {
specialPosition_ += separation; specialPosition_ += separation;
} }
sitePositions_ = sitePositions_ + separation;
} }

View File

@ -182,13 +182,6 @@ private:
// Private data // Private data
//- Be careful with the ordering of data.
// It has an impact on binary transfer:
// -# Put the largest data members 1st
// -# Pair up labels,
// -# Don't go scalar-label, scalar-label, because in 64bit mode,
// the labels will be padded by 4bytes.
tensor Q_; tensor Q_;
vector v_; vector v_;

View File

@ -242,11 +242,17 @@ void Foam::molecule::writeFields(const Cloud<molecule>& mC)
orientation2.write(); orientation2.write();
orientation3.write(); orientation3.write();
Info<< "writeFields " << mC.name() << endl;
if (isA<moleculeCloud>(mC))
{
const moleculeCloud& m = dynamic_cast<const moleculeCloud&>(mC); const moleculeCloud& m = dynamic_cast<const moleculeCloud&>(mC);
m.writeXYZ m.writeXYZ
( (
m.mesh().time().timePath() + "/lagrangian" + "/moleculeCloud.xmol" m.mesh().time().timePath()/cloud::prefix/"moleculeCloud.xmol"
); );
}
} }

View File

@ -123,23 +123,23 @@ void Foam::moleculeCloud::buildCellOccupancy()
{ {
cellOccupancy_[cO].shrink(); cellOccupancy_[cO].shrink();
} }
il_.ril().referMolecules(cellOccupancy());
} }
void Foam::moleculeCloud::calculatePairForce() void Foam::moleculeCloud::calculatePairForce()
{ {
iterator mol(this->begin()); PstreamBuffers pBufs(Pstream::nonBlocking);
// Start sending referred data
il_.sendReferredData(cellOccupancy(), pBufs);
molecule* molI = NULL;
molecule* molJ = NULL;
{ {
// Real-Real interactions // Real-Real interactions
molecule* molI = &mol(); const labelListList& dil = il_.dil();
molecule* molJ = &mol();
const directInteractionList& dil(il_.dil());
forAll(dil, d) forAll(dil, d)
{ {
@ -149,54 +149,62 @@ void Foam::moleculeCloud::calculatePairForce()
forAll(dil[d], interactingCells) forAll(dil[d], interactingCells)
{ {
List< molecule* > cellJ = List<molecule*> cellJ =
cellOccupancy_[dil[d][interactingCells]]; cellOccupancy_[dil[d][interactingCells]];
forAll(cellJ, cellJMols) forAll(cellJ, cellJMols)
{ {
molJ = cellJ[cellJMols]; molJ = cellJ[cellJMols];
evaluatePair(molI, molJ); evaluatePair(*molI, *molJ);
} }
} }
forAll(cellOccupancy_[d],cellIOtherMols) forAll(cellOccupancy_[d], cellIOtherMols)
{ {
molJ = cellOccupancy_[d][cellIOtherMols]; molJ = cellOccupancy_[d][cellIOtherMols];
if (molJ > molI) if (molJ > molI)
{ {
evaluatePair(molI, molJ); evaluatePair(*molI, *molJ);
} }
} }
} }
} }
} }
// Receive referred data
il_.receiveReferredData(pBufs);
{ {
// Real-Referred interactions // Real-Referred interactions
molecule* molK = &mol(); const labelListList& ril = il_.ril();
referredCellList& ril(il_.ril()); List<IDLList<molecule> >& referredMols = il_.referredParticles();
forAll(ril, r) forAll(ril, r)
{ {
const List<label>& realCells = ril[r].realCellsForInteraction(); const List<label>& realCells = ril[r];
forAll(ril[r], refMols) IDLList<molecule>& refMols = referredMols[r];
forAllIter
(
IDLList<molecule>,
refMols,
refMol
)
{ {
referredMolecule* molL(&(ril[r][refMols]));
forAll(realCells, rC) forAll(realCells, rC)
{ {
List<molecule*> cellK = cellOccupancy_[realCells[rC]]; List<molecule*> cellI = cellOccupancy_[realCells[rC]];
forAll(cellK, cellKMols) forAll(cellI, cellIMols)
{ {
molK = cellK[cellKMols]; molI = cellI[cellIMols];
evaluatePair(molK, molL); evaluatePair(*molI, refMol());
} }
} }
} }
@ -259,18 +267,14 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
buildCellOccupancy(); buildCellOccupancy();
// Real-Real interaction // Real-Real interaction
iterator mol(this->begin());
molecule* molI = NULL;
molecule* molJ = NULL;
{ {
mol = this->begin();
molecule* molI = &mol();
molecule* molJ = &mol();
DynamicList<molecule*> molsToDelete; DynamicList<molecule*> molsToDelete;
const directInteractionList& dil(il_.dil()); const labelListList& dil(il_.dil());
forAll(dil, d) forAll(dil, d)
{ {
@ -280,14 +284,14 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
forAll(dil[d], interactingCells) forAll(dil[d], interactingCells)
{ {
List< molecule* > cellJ = List<molecule*> cellJ =
cellOccupancy_[dil[d][interactingCells]]; cellOccupancy_[dil[d][interactingCells]];
forAll(cellJ, cellJMols) forAll(cellJ, cellJMols)
{ {
molJ = cellJ[cellJMols]; molJ = cellJ[cellJMols];
if (evaluatePotentialLimit(molI, molJ)) if (evaluatePotentialLimit(*molI, *molJ))
{ {
label idI = molI->id(); label idI = molI->id();
@ -314,13 +318,13 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
} }
} }
forAll(cellOccupancy_[d],cellIOtherMols) forAll(cellOccupancy_[d], cellIOtherMols)
{ {
molJ = cellOccupancy_[d][cellIOtherMols]; molJ = cellOccupancy_[d][cellIOtherMols];
if (molJ > molI) if (molJ > molI)
{ {
if (evaluatePotentialLimit(molI, molJ)) if (evaluatePotentialLimit(*molI, *molJ))
{ {
label idI = molI->id(); label idI = molI->id();
@ -355,79 +359,81 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
buildCellOccupancy(); buildCellOccupancy();
PstreamBuffers pBufs(Pstream::nonBlocking);
// Start sending referred data
il_.sendReferredData(cellOccupancy(), pBufs);
// Receive referred data
il_.receiveReferredData(pBufs);
// Real-Referred interaction // Real-Referred interaction
{ {
molecule* molK = &mol();
DynamicList<molecule*> molsToDelete; DynamicList<molecule*> molsToDelete;
referredCellList& ril(il_.ril()); const labelListList& ril(il_.ril());
List<IDLList<molecule> >& referredMols = il_.referredParticles();
forAll(ril, r) forAll(ril, r)
{ {
referredCell& refCell = ril[r]; IDLList<molecule>& refMols = referredMols[r];
forAll(refCell, refMols) forAllIter
(
IDLList<molecule>,
refMols,
refMol
)
{ {
referredMolecule* molL = &(refCell[refMols]); molJ = &refMol();
List <label> realCells = refCell.realCellsForInteraction(); const List<label>& realCells = ril[r];
forAll(realCells, rC) forAll(realCells, rC)
{ {
label cellK = realCells[rC]; label cellI = realCells[rC];
List<molecule*> cellKMols = cellOccupancy_[cellK]; List<molecule*> cellIMols = cellOccupancy_[cellI];
forAll(cellKMols, cKM) forAll(cellIMols, cIM)
{ {
molK = cellKMols[cKM]; molI = cellIMols[cIM];
if (evaluatePotentialLimit(molK, molL)) if (evaluatePotentialLimit(*molI, *molJ))
{ {
label idK = molK->id(); label idI = molI->id();
label idL = molL->id(); label idJ = molJ->id();
if if
( (
findIndex(pot_.removalOrder(), idK) findIndex(pot_.removalOrder(), idI)
< findIndex(pot_.removalOrder(), idL) < findIndex(pot_.removalOrder(), idJ)
) )
{ {
if (findIndex(molsToDelete, molK) == -1) if (findIndex(molsToDelete, molI) == -1)
{ {
molsToDelete.append(molK); molsToDelete.append(molI);
} }
} }
else if else if
( (
findIndex(pot_.removalOrder(), idK) findIndex(pot_.removalOrder(), idI)
== findIndex(pot_.removalOrder(), idL) == findIndex(pot_.removalOrder(), idJ)
) )
{ {
if // Remove one of the molecules
( // arbitrarily, assuring that a
Pstream::myProcNo() == refCell.sourceProc() // consistent decision is made for
&& cellK <= refCell.sourceCell() // both real-referred pairs.
)
{
if (findIndex(molsToDelete, molK) == -1)
{
molsToDelete.append(molK);
}
}
else if if (molI->origId() > molJ->origId())
(
Pstream::myProcNo() < refCell.sourceProc()
)
{ {
if (findIndex(molsToDelete, molK) == -1) if (findIndex(molsToDelete, molI) == -1)
{ {
molsToDelete.append(molK); molsToDelete.append(molI);
} }
} }
} }
@ -445,6 +451,12 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
buildCellOccupancy(); buildCellOccupancy();
// Start sending referred data
il_.sendReferredData(cellOccupancy(), pBufs);
// Receive referred data
il_.receiveReferredData(pBufs);
label molsRemoved = initialSize - this->size(); label molsRemoved = initialSize - this->size();
if (Pstream::parRun()) if (Pstream::parRun())
@ -722,10 +734,13 @@ void Foam::moleculeCloud::initialiseMolecules
+ latticePositions[p]; + latticePositions[p];
point globalPosition = point globalPosition =
anchor + (R & (latticeCellShape & latticePosition)); anchor
+ (R & (latticeCellShape & latticePosition));
partOfLayerInBounds = partOfLayerInBounds = mesh_.bounds().contains
mesh_.bounds().contains(globalPosition); (
globalPosition
);
label cell = mesh_.findCell(globalPosition); label cell = mesh_.findCell(globalPosition);
@ -787,12 +802,21 @@ void Foam::moleculeCloud::initialiseMolecules
) )
+ latticePositions[p]; + latticePositions[p];
point globalPosition = anchor point globalPosition =
+ (R anchor
&(latticeCellShape & latticePosition)); + (
R
& (
latticeCellShape
& latticePosition
)
);
partOfLayerInBounds = partOfLayerInBounds =
mesh_.bounds().contains(globalPosition); mesh_.bounds().contains
(
globalPosition
);
label cell = mesh_.findCell label cell = mesh_.findCell
( (
@ -848,12 +872,21 @@ void Foam::moleculeCloud::initialiseMolecules
) )
+ latticePositions[p]; + latticePositions[p];
point globalPosition = anchor point globalPosition =
+ (R anchor
&(latticeCellShape & latticePosition)); + (
R
& (
latticeCellShape
& latticePosition
)
);
partOfLayerInBounds = partOfLayerInBounds =
mesh_.bounds().contains(globalPosition); mesh_.bounds().contains
(
globalPosition
);
label cell = mesh_.findCell label cell = mesh_.findCell
( (
@ -1041,7 +1074,7 @@ Foam::moleculeCloud::moleculeCloud
mesh_(mesh), mesh_(mesh),
pot_(pot), pot_(pot),
cellOccupancy_(mesh_.nCells()), cellOccupancy_(mesh_.nCells()),
il_(mesh_, pot_.pairPotentials().rCutMaxSqr(), false), il_(mesh_, pot_.pairPotentials().rCutMax(), true),
constPropList_(), constPropList_(),
rndGen_(clock::getTime()) rndGen_(clock::getTime())
{ {
@ -1071,7 +1104,7 @@ Foam::moleculeCloud::moleculeCloud
Cloud<molecule>(mesh, "moleculeCloud", false), Cloud<molecule>(mesh, "moleculeCloud", false),
mesh_(mesh), mesh_(mesh),
pot_(pot), pot_(pot),
il_(mesh_), il_(mesh_, 0.0, false),
constPropList_(), constPropList_(),
rndGen_(clock::getTime()) rndGen_(clock::getTime())
{ {
@ -1137,7 +1170,7 @@ void Foam::moleculeCloud::applyConstraintsAndThermostats
) )
{ {
scalar temperatureCorrectionFactor = scalar temperatureCorrectionFactor =
sqrt(targetTemperature/measuredTemperature); sqrt(targetTemperature/max(VSMALL, measuredTemperature));
Info<< "----------------------------------------" << nl Info<< "----------------------------------------" << nl
<< "Temperature equilibration" << nl << "Temperature equilibration" << nl

View File

@ -40,7 +40,7 @@ SourceFiles
#include "molecule.H" #include "molecule.H"
#include "IOdictionary.H" #include "IOdictionary.H"
#include "potential.H" #include "potential.H"
#include "interactionLists.H" #include "InteractionLists.H"
#include "labelVector.H" #include "labelVector.H"
#include "Random.H" #include "Random.H"
#include "fileName.H" #include "fileName.H"
@ -69,7 +69,7 @@ private:
List<DynamicList<molecule*> > cellOccupancy_; List<DynamicList<molecule*> > cellOccupancy_;
interactionLists il_; InteractionLists<molecule> il_;
List<molecule::constantProperties> constPropList_; List<molecule::constantProperties> constPropList_;
@ -89,26 +89,14 @@ private:
inline void evaluatePair inline void evaluatePair
( (
molecule* molI, molecule& molI,
molecule* molJ molecule& molJ
);
inline void evaluatePair
(
molecule* molReal,
referredMolecule* molRef
); );
inline bool evaluatePotentialLimit inline bool evaluatePotentialLimit
( (
molecule* molI, molecule& molI,
molecule* molJ molecule& molJ
) const;
inline bool evaluatePotentialLimit
(
molecule* molReal,
referredMolecule* molRef
) const; ) const;
void calculateTetherForce(); void calculateTetherForce();
@ -197,7 +185,7 @@ public:
inline const List<DynamicList<molecule*> >& cellOccupancy() const; inline const List<DynamicList<molecule*> >& cellOccupancy() const;
inline const interactionLists& il() const; inline const InteractionLists<molecule>& il() const;
inline const List<molecule::constantProperties> constProps() const; inline const List<molecule::constantProperties> constProps() const;

View File

@ -31,17 +31,17 @@ using namespace Foam::constant;
inline void Foam::moleculeCloud::evaluatePair inline void Foam::moleculeCloud::evaluatePair
( (
molecule* molI, molecule& molI,
molecule* molJ molecule& molJ
) )
{ {
const pairPotentialList& pairPot = pot_.pairPotentials(); const pairPotentialList& pairPot = pot_.pairPotentials();
const pairPotential& electrostatic = pairPot.electrostatic(); const pairPotential& electrostatic = pairPot.electrostatic();
label idI = molI->id(); label idI = molI.id();
label idJ = molJ->id(); label idJ = molJ.id();
const molecule::constantProperties& constPropI(constProps(idI)); const molecule::constantProperties& constPropI(constProps(idI));
@ -70,7 +70,7 @@ inline void Foam::moleculeCloud::evaluatePair
if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ]) if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
{ {
vector rsIsJ = vector rsIsJ =
molI->sitePositions()[sI] - molJ->sitePositions()[sJ]; molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
scalar rsIsJMagSq = magSqr(rsIsJ); scalar rsIsJMagSq = magSqr(rsIsJ);
@ -82,38 +82,34 @@ inline void Foam::moleculeCloud::evaluatePair
(rsIsJ/rsIsJMag) (rsIsJ/rsIsJMag)
*pairPot.force(idsI, idsJ, rsIsJMag); *pairPot.force(idsI, idsJ, rsIsJMag);
molI->siteForces()[sI] += fsIsJ; molI.siteForces()[sI] += fsIsJ;
molJ->siteForces()[sJ] += -fsIsJ; molJ.siteForces()[sJ] += -fsIsJ;
scalar potentialEnergy scalar potentialEnergy
( (
pairPot.energy(idsI, idsJ, rsIsJMag) pairPot.energy(idsI, idsJ, rsIsJMag)
); );
molI->potentialEnergy() += 0.5*potentialEnergy; molI.potentialEnergy() += 0.5*potentialEnergy;
molJ->potentialEnergy() += 0.5*potentialEnergy; molJ.potentialEnergy() += 0.5*potentialEnergy;
vector rIJ = molI->position() - molJ->position(); vector rIJ = molI.position() - molJ.position();
tensor virialContribution = tensor virialContribution =
(rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq; (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
molI->rf() += virialContribution; molI.rf() += virialContribution;
molJ->rf() += virialContribution; molJ.rf() += virialContribution;
// molI->rf() += rsIsJ * fsIsJ;
// molJ->rf() += rsIsJ * fsIsJ;
} }
} }
if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ]) if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
{ {
vector rsIsJ = vector rsIsJ =
molI->sitePositions()[sI] - molJ->sitePositions()[sJ]; molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
scalar rsIsJMagSq = magSqr(rsIsJ); scalar rsIsJMagSq = magSqr(rsIsJ);
@ -129,149 +125,26 @@ inline void Foam::moleculeCloud::evaluatePair
(rsIsJ/rsIsJMag) (rsIsJ/rsIsJMag)
*chargeI*chargeJ*electrostatic.force(rsIsJMag); *chargeI*chargeJ*electrostatic.force(rsIsJMag);
molI->siteForces()[sI] += fsIsJ; molI.siteForces()[sI] += fsIsJ;
molJ->siteForces()[sJ] += -fsIsJ; molJ.siteForces()[sJ] += -fsIsJ;
scalar potentialEnergy = scalar potentialEnergy =
chargeI*chargeJ chargeI*chargeJ
*electrostatic.energy(rsIsJMag); *electrostatic.energy(rsIsJMag);
molI->potentialEnergy() += 0.5*potentialEnergy; molI.potentialEnergy() += 0.5*potentialEnergy;
molJ->potentialEnergy() += 0.5*potentialEnergy; molJ.potentialEnergy() += 0.5*potentialEnergy;
vector rIJ = molI->position() - molJ->position(); vector rIJ = molI.position() - molJ.position();
tensor virialContribution = tensor virialContribution =
(rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq; (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
molI->rf() += virialContribution; molI.rf() += virialContribution;
molJ->rf() += virialContribution; molJ.rf() += virialContribution;
// molI->rf() += rsIsJ * fsIsJ;
// molJ->rf() += rsIsJ * fsIsJ;
}
}
}
}
}
inline void Foam::moleculeCloud::evaluatePair
(
molecule* molReal,
referredMolecule* molRef
)
{
const pairPotentialList& pairPot = pot_.pairPotentials();
const pairPotential& electrostatic = pairPot.electrostatic();
label idReal = molReal->id();
label idRef = molRef->id();
const molecule::constantProperties& constPropReal(constProps(idReal));
const molecule::constantProperties& constPropRef(constProps(idRef));
List<label> siteIdsReal = constPropReal.siteIds();
List<label> siteIdsRef = constPropRef.siteIds();
List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
forAll(siteIdsReal, sReal)
{
label idsReal(siteIdsReal[sReal]);
forAll(siteIdsRef, sRef)
{
label idsRef(siteIdsRef[sRef]);
if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
{
vector rsRealsRef =
molReal->sitePositions()[sReal]
- molRef->sitePositions()[sRef];
scalar rsRealsRefMagSq = magSqr(rsRealsRef);
if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
{
scalar rsRealsRefMag = mag(rsRealsRef);
vector fsRealsRef =
(rsRealsRef/rsRealsRefMag)
*pairPot.force(idsReal, idsRef, rsRealsRefMag);
molReal->siteForces()[sReal] += fsRealsRef;
scalar potentialEnergy
(
pairPot.energy(idsReal, idsRef, rsRealsRefMag)
);
molReal->potentialEnergy() += 0.5*potentialEnergy;
vector rRealRef = molReal->position() - molRef->position();
molReal->rf() +=
(rsRealsRef*fsRealsRef)
*(rsRealsRef & rRealRef)
/rsRealsRefMagSq;
// molReal->rf() += rsRealsRef * fsRealsRef;
}
}
if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
{
vector rsRealsRef =
molReal->sitePositions()[sReal]
- molRef->sitePositions()[sRef];
scalar rsRealsRefMagSq = magSqr(rsRealsRef);
if (rsRealsRefMagSq <= electrostatic.rCutSqr())
{
scalar rsRealsRefMag = mag(rsRealsRef);
scalar chargeReal = constPropReal.siteCharges()[sReal];
scalar chargeRef = constPropRef.siteCharges()[sRef];
vector fsRealsRef =
(rsRealsRef/rsRealsRefMag)
*chargeReal*chargeRef
*electrostatic.force(rsRealsRefMag);
molReal->siteForces()[sReal] += fsRealsRef;
scalar potentialEnergy =
chargeReal*chargeRef
*electrostatic.energy(rsRealsRefMag);
molReal->potentialEnergy() += 0.5*potentialEnergy;
vector rRealRef = molReal->position() - molRef->position();
molReal->rf() +=
(rsRealsRef*fsRealsRef)
*(rsRealsRef & rRealRef)
/rsRealsRefMagSq;
// molReal->rf() += rsRealsRef * fsRealsRef;
} }
} }
} }
@ -281,17 +154,17 @@ inline void Foam::moleculeCloud::evaluatePair
inline bool Foam::moleculeCloud::evaluatePotentialLimit inline bool Foam::moleculeCloud::evaluatePotentialLimit
( (
molecule* molI, molecule& molI,
molecule* molJ molecule& molJ
) const ) const
{ {
const pairPotentialList& pairPot = pot_.pairPotentials(); const pairPotentialList& pairPot = pot_.pairPotentials();
const pairPotential& electrostatic = pairPot.electrostatic(); const pairPotential& electrostatic = pairPot.electrostatic();
label idI = molI->id(); label idI = molI.id();
label idJ = molJ->id(); label idJ = molJ.id();
const molecule::constantProperties& constPropI(constProps(idI)); const molecule::constantProperties& constPropI(constProps(idI));
@ -320,7 +193,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ]) if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
{ {
vector rsIsJ = vector rsIsJ =
molI->sitePositions()[sI] - molJ->sitePositions()[sJ]; molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
scalar rsIsJMagSq = magSqr(rsIsJ); scalar rsIsJMagSq = magSqr(rsIsJ);
@ -369,7 +242,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ]) if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
{ {
vector rsIsJ = vector rsIsJ =
molI->sitePositions()[sI] - molJ->sitePositions()[sJ]; molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
scalar rsIsJMagSq = magSqr(rsIsJ); scalar rsIsJMagSq = magSqr(rsIsJ);
@ -422,157 +295,6 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
} }
inline bool Foam::moleculeCloud::evaluatePotentialLimit
(
molecule* molReal,
referredMolecule* molRef
) const
{
const pairPotentialList& pairPot = pot_.pairPotentials();
const pairPotential& electrostatic = pairPot.electrostatic();
label idReal = molReal->id();
label idRef = molRef->id();
const molecule::constantProperties& constPropReal(constProps(idReal));
const molecule::constantProperties& constPropRef(constProps(idRef));
List<label> siteIdsReal = constPropReal.siteIds();
List<label> siteIdsRef = constPropRef.siteIds();
List<bool> pairPotentialSitesReal = constPropReal.pairPotentialSites();
List<bool> electrostaticSitesReal = constPropReal.electrostaticSites();
List<bool> pairPotentialSitesRef = constPropRef.pairPotentialSites();
List<bool> electrostaticSitesRef = constPropRef.electrostaticSites();
forAll(siteIdsReal, sReal)
{
label idsReal(siteIdsReal[sReal]);
forAll(siteIdsRef, sRef)
{
label idsRef(siteIdsRef[sRef]);
if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
{
vector rsRealsRef =
molReal->sitePositions()[sReal]
- molRef->sitePositions()[sRef];
scalar rsRealsRefMagSq = magSqr(rsRealsRef);
if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
{
scalar rsRealsRefMag = mag(rsRealsRef);
// Guard against pairPot.energy being evaluated
// if rRealRefMag < SMALL. A floating point exception will
// happen otherwise.
if (rsRealsRefMag < SMALL)
{
WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
<< "Molecule site pair closer than "
<< SMALL
<< ": mag separation = " << rsRealsRefMag
<< ". These may have been placed on top of each"
<< " other by a rounding error in mdInitialise in"
<< " parallel or a block filled with molecules"
<< " twice. Removing one of the molecules."
<< endl;
return true;
}
// Guard against pairPot.energy being evaluated if
// rRealRefMag < rMin. A tabulation lookup error will occur
// otherwise.
if (rsRealsRefMag < pairPot.rMin(idsReal, idsRef))
{
return true;
}
if
(
mag(pairPot.energy(idsReal, idsRef, rsRealsRefMag))
> pot_.potentialEnergyLimit()
)
{
return true;
};
}
}
if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
{
vector rsRealsRef =
molReal->sitePositions()[sReal]
- molRef->sitePositions()[sRef];
scalar rsRealsRefMagSq = magSqr(rsRealsRef);
if (pairPot.rCutMaxSqr(rsRealsRefMagSq))
{
scalar rsRealsRefMag = mag(rsRealsRef);
// Guard against pairPot.energy being evaluated
// if rRealRefMag < SMALL. A floating point exception will
// happen otherwise.
if (rsRealsRefMag < SMALL)
{
WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
<< "Molecule site pair closer than "
<< SMALL
<< ": mag separation = " << rsRealsRefMag
<< ". These may have been placed on top of each"
<< " other by a rounding error in mdInitialise in"
<< " parallel or a block filled with molecules"
<< " twice. Removing one of the molecules."
<< endl;
return true;
}
if (rsRealsRefMag < electrostatic.rMin())
{
return true;
}
scalar chargeReal = constPropReal.siteCharges()[sReal];
scalar chargeRef = constPropRef.siteCharges()[sRef];
if
(
mag
(
chargeReal
*chargeRef
*electrostatic.energy(rsRealsRefMag)
)
> pot_.potentialEnergyLimit()
)
{
return true;
};
}
}
}
}
return false;
}
inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
( (
scalar temperature, scalar temperature,
@ -638,7 +360,7 @@ inline const Foam::List<Foam::DynamicList<Foam::molecule*> >&
} }
inline const Foam::interactionLists& inline const Foam::InteractionLists<Foam::molecule>&
Foam::moleculeCloud::il() const Foam::moleculeCloud::il() const
{ {
return il_; return il_;

View File

@ -0,0 +1,37 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
periodicX
{
type cyclic;
}
periodicY
{
type cyclic;
}
periodicZ
{
type cyclic;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,37 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
periodicX
{
type cyclic;
}
periodicY
{
type cyclic;
}
periodicZ
{
type cyclic;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,55 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
sectionAEnd
{
type fixedValue;
value uniform (0 0 0);
}
sectionCEnd
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type fixedValue;
value uniform (0 0 0);
}
back
{
type fixedValue;
value uniform (0 0 0);
}
top
{
type fixedValue;
value uniform (0 0 0);
}
bottom
{
type fixedValue;
value uniform (0 0 0);
}
}
// ************************************************************************* //