Reading individual sites from moleculeProperties dictionary entries, WIP.

This commit is contained in:
graham 2008-10-07 20:14:10 +01:00
parent 4a6fc137c2
commit 69ab38bc8f
8 changed files with 187 additions and 38 deletions

View File

@ -31,7 +31,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
//#include "md.H"
#include "md.H"
#include "potential.H"
int main(int argc, char *argv[])
@ -43,6 +43,8 @@ int main(int argc, char *argv[])
potential pot(mesh);
moleculeCloud molecules(mesh, pot);
Info << "\nStarting time loop\n" << endl;
while (runTime.run())

View File

@ -53,7 +53,7 @@ int main(int argc, char *argv[])
# include "createAutoCorrelationFunctions.H"
label nAveragingSteps = 0;
Info << "\nStarting time loop\n" << endl;
while (runTime.run())
@ -87,6 +87,6 @@ int main(int argc, char *argv[])
# include "calculateTransportProperties.H"
Info << "End\n" << endl;
return(0);
}

View File

@ -80,16 +80,19 @@ public:
// Private data
List<vector> siteReferencePositions_;
List<vector> siteReferencePositions_;
List<scalar> siteCharges_;
List<scalar> siteCharges_;
List<label> siteIds_;
List<label> siteIds_;
diagTensor momentOfInertia_;
diagTensor momentOfInertia_;
scalar mass_;
scalar mass_;
// Private Member Functions
void checkSiteListSizes() const;
public:
@ -105,6 +108,9 @@ public:
scalar mass
);
//- Constructor
inline constantProperties(const dictionary& dict);
// Member functions
inline const List<vector>& siteReferencePositions() const;
@ -112,6 +118,7 @@ public:
inline const List<scalar>& siteCharges() const;
inline const List<label>& siteIds() const;
inline List<label>& siteIds();
inline const diagTensor& momentOfInertia() const;

View File

@ -33,9 +33,7 @@ inline Foam::molecule::constantProperties::constantProperties()
siteIds_(List<label>(0)),
momentOfInertia_(diagTensor(0, 0, 0)),
mass_(0)
{
}
{}
inline Foam::molecule::constantProperties::constantProperties
@ -53,16 +51,24 @@ inline Foam::molecule::constantProperties::constantProperties
momentOfInertia_(momentOfInertia),
mass_(mass)
{
if
(
siteIds_.size() != siteReferencePositions_.size()
|| siteIds_.size() != siteCharges_.size()
)
{
FatalErrorIn("molecule::constantProperties::constantProperties")
<< "Sizes of site id, charge and referencePositions are not the same. "
<< abort(FatalError);
}
checkSiteListSizes();
}
inline Foam::molecule::constantProperties::constantProperties
(
const dictionary& dict
)
:
siteReferencePositions_(dict.lookup("siteReferencePositions")),
siteCharges_(dict.lookup("siteCharges")),
siteIds_(List<word>(dict.lookup("siteIds")).size(), -1),
momentOfInertia_(dict.lookup("momentOfInertia")),
mass_(readScalar(dict.lookup("mass")))
{
checkSiteListSizes();
Info<< nl << "siteIds not set yet." << endl;
}
@ -101,6 +107,25 @@ inline Foam::molecule::molecule
}
// * * * constantProperties Private Member Functions * * * * * * * * * * * * //
inline void Foam::molecule::constantProperties::checkSiteListSizes() const
{
if
(
siteIds_.size() != siteReferencePositions_.size()
|| siteIds_.size() != siteCharges_.size()
)
{
FatalErrorIn("molecule::constantProperties::checkSiteListSizes")
<< "Sizes of site id, charge and "
<< "referencePositions are not the same. "
<< nl << abort(FatalError);
}
}
// * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
inline const Foam::List<Foam::vector>&
@ -124,6 +149,13 @@ inline const Foam::List<Foam::label>&
}
inline Foam::List<Foam::label>&
Foam::molecule::constantProperties::siteIds()
{
return siteIds_;
}
inline const Foam::diagTensor&
Foam::molecule::constantProperties::momentOfInertia() const
{

View File

@ -44,6 +44,57 @@ Foam::scalar Foam::moleculeCloud::vacuumPermittivity = 8.854187817e-12;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::moleculeCloud::buildConstProps()
{
Info<< nl << "Reading moleculeProperties dictionary." << endl;
const List<word>& idList(pot_.idList());
constPropList_.setSize(idList.size());
const List<word>& allSiteIdNames(pot_.allSiteIdNames());
forAll(idList, i)
{
const word& id(idList[i]);
const dictionary& molDict(moleculePropertiesDict.subDict(id));
List<word> siteIdNames = molDict.lookup("siteIds");
List<label> siteIds(siteIdNames.size());
forAll(siteIdNames, sI)
{
const word& siteId = siteIdNames[sI];
siteIds[sI] = findIndex(allSiteIdNames, siteId);
if (siteIds[sI] == -1)
{
FatalErrorIn("moleculeCloud.C") << nl
<< siteId << " site not found."
<< nl << abort(FatalError);
}
}
molecule::constantProperties& constProp = constPropList_[i];
constProp = molecule::constantProperties(molDict);
constProp.siteIds() = siteIds;
Info<< constPropList_[i].siteReferencePositions()
<< nl << constPropList_[i].siteCharges()
<< nl << constPropList_[i].siteIds()
<< nl << constPropList_[i].momentOfInertia()
<< nl << constPropList_[i].mass()
<< nl << constPropList_[i].nSites()
<< endl;
}
}
void Foam::moleculeCloud::buildCellOccupancy()
{
forAll(cellOccupancy_, cO)
@ -422,6 +473,8 @@ Foam::moleculeCloud::moleculeCloud
{
molecule::readFields(*this);
buildConstProps();
buildCellOccupancy();
removeHighEnergyOverlaps();

View File

@ -73,6 +73,8 @@ private:
// Private Member Functions
void buildConstProps();
//- Determine which molecules are in which cells
void buildCellOccupancy();
@ -141,7 +143,8 @@ public:
const potential& pot
);
//- Construct given polyMesh and fields of components. Intended for use
//- Construct given polyMesh and fields of components.
//- Intended for use
//- by the molConfig utility
// moleculeCloud
// (
@ -184,7 +187,8 @@ public:
inline const List<molecule::constantProperties> constProps() const;
inline const molecule::constantProperties& constProps(label id) const;
inline const molecule::constantProperties&
constProps(label id) const;
// Member Operators

View File

@ -101,7 +101,7 @@ inline void Foam::moleculeCloud::evaluatePair
// Acceleration increment for molReal
molReal->a() += fRealRef/massReal;
// Adding a contribution 19f 1/2 of the potential energy
// Adding a contribution of 1/2 of the potential energy
// from this interaction
molReal->potentialEnergy() +=
0.5*pairPot.energy(idReal, idRef, rRealRefMag);

View File

@ -31,7 +31,22 @@ License
Foam::potential::potential(const polyMesh& mesh)
:
mesh_(mesh)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::potential::~potential()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::potential::potential::readPotentialDict
(
const List<word>& allSiteIdNames
)
{
Info<< nl << "Reading potential dictionary:" << endl;
@ -49,6 +64,52 @@ Foam::potential::potential(const polyMesh& mesh)
idList_ = List<word>(idListDict.lookup("idList"));
IOdictionary moleculePropertiesDict
(
IOobject
(
"moleculeProperties",
mesh_.time().constant(),
mesh_,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
);
DynamicList<word> allSiteIdNames;
forAll(idList, i)
{
const word& id(idList[i]);
if (!moleculePropertiesDict.found(id))
{
FatalErrorIn("potential.C") << nl
<< id << " molecule subDict not found"
<< nl << abort(FatalError);
}
const dictionary& molDict(moleculePropertiesDict.subDict(id));
List<word> siteIdNames = molDict.lookup("siteIds");
forAll(siteIdNames, sI)
{
const word& siteId = siteIdNames[sI];
if(findIndex(allSiteIdNames, siteId) == -1)
{
allSiteIdNames.append(siteId);
}
}
}
// CREATE allSiteIdNames_ DATA MEMBER
allSiteIdNames_.transfer(allSiteIdNames.shrink());
Info<< nl << "Unique site ids found: " << allSiteIdNames << endl;
List<word> tetherIdList(0);
if (idListDict.found("tetherIdList"))
@ -85,7 +146,7 @@ Foam::potential::potential(const polyMesh& mesh)
}
}
// ****************************************************************************
// *************************************************************************
// Pair potentials
if (!potentialDict.found("pair"))
@ -99,7 +160,7 @@ Foam::potential::potential(const polyMesh& mesh)
pairPotentials_.buildPotentials(idList_, pairDict, mesh_);
// ****************************************************************************
// *************************************************************************
// Tether potentials
if (tetherIdList.size())
@ -116,7 +177,7 @@ Foam::potential::potential(const polyMesh& mesh)
tetherPotentials_.buildPotentials(idList_, tetherDict, tetherIdList);
}
// ****************************************************************************
// *************************************************************************
// External Forces
gravity_ = vector::zero;
@ -128,7 +189,7 @@ Foam::potential::potential(const polyMesh& mesh)
const dictionary& externalDict = potentialDict.subDict("external");
// ************************************************************************
// *********************************************************************
// gravity
if (externalDict.found("gravity"))
@ -140,16 +201,6 @@ Foam::potential::potential(const polyMesh& mesh)
Info << nl << tab << "gravity = " << gravity_ << endl;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::potential::~potential()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// ************************************************************************* //