31 siteMasses_(
List<scalar>(0)),
32 siteCharges_(
List<scalar>(0)),
34 pairPotentialSites_(
List<bool>(false)),
35 electrostaticSites_(
List<bool>(false)),
46 siteReferencePositions_(dict.
lookup(
"siteReferencePositions")),
47 siteMasses_(dict.
lookup(
"siteMasses")),
48 siteCharges_(dict.
lookup(
"siteCharges")),
50 pairPotentialSites_(),
51 electrostaticSites_(),
57 setInteracionSiteBools
63 mass_ =
sum(siteMasses_);
70 forAll(siteReferencePositions_, i)
72 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
75 centreOfMass /= mass_;
77 siteReferencePositions_ -= centreOfMass;
79 if (siteIds_.
size() == 1)
83 siteReferencePositions_[0] =
Zero;
87 else if (linearMoleculeTest())
93 vector dir = siteReferencePositions_[1] - siteReferencePositions_[0];
99 siteReferencePositions_ = (Q & siteReferencePositions_);
106 forAll(siteReferencePositions_, i)
108 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
111 centreOfMass /= mass_;
113 siteReferencePositions_ -= centreOfMass;
117 forAll(siteReferencePositions_, i)
119 const vector&
p(siteReferencePositions_[i]);
140 forAll(siteReferencePositions_, i)
142 const vector&
p(siteReferencePositions_[i]);
144 momOfInertia += siteMasses_[i]*
tensor 146 p.
y()*p.
y() + p.
z()*p.
z(), -p.
x()*p.
y(), -p.
x()*p.
z(),
147 -p.
y()*p.
x(), p.
x()*p.
x() + p.
z()*p.
z(), -p.
y()*p.
z(),
148 -p.
z()*p.
x(), -p.
z()*p.
y(), p.
x()*p.
x() + p.
y()*p.
y()
155 <<
"An eigenvalue of the inertia tensor is zero, the molecule " 156 << dict.
name() <<
" is not a valid 6DOF shape." 175 siteReferencePositions_ = (Q & siteReferencePositions_);
184 forAll(siteReferencePositions_, i)
186 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
189 centreOfMass /= mass_;
191 siteReferencePositions_ -= centreOfMass;
198 forAll(siteReferencePositions_, i)
200 const vector&
p(siteReferencePositions_[i]);
202 momOfInertia += siteMasses_[i]*
tensor 204 p.
y()*p.
y() + p.
z()*p.
z(), -p.
x()*p.
y(), -p.
x()*p.
z(),
205 -p.
y()*p.
x(), p.
x()*p.
x() + p.
z()*p.
z(), -p.
y()*p.
z(),
206 -p.
z()*p.
x(), -p.
z()*p.
y(), p.
x()*p.
x() + p.
y()*p.
y()
225 const label tetFacei,
232 const vector& specialPosition,
239 particle(mesh, position, celli, tetFacei, tetPtI),
245 specialPosition_(specialPosition),
246 potentialEnergy_(0.0),
251 sitePositions_(constProps.
nSites())
259 inline void Foam::molecule::constantProperties::checkSiteListSizes()
const 263 siteIds_.
size() != siteReferencePositions_.
size()
264 || siteIds_.
size() != siteCharges_.
size()
268 <<
"Sizes of site id, charge and " 269 <<
"referencePositions are not the same. " 275 inline void Foam::molecule::constantProperties::setInteracionSiteBools
287 const word&
id(siteIds[i]);
289 pairPotentialSites_[i] = (
findIndex(pairPotSiteIds,
id) > -1);
291 electrostaticSites_[i] = (
mag(siteCharges_[i]) > VSMALL);
296 inline bool Foam::molecule::constantProperties::linearMoleculeTest()
const 298 if (siteIds_.
size() == 2)
303 vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0];
305 refDir /=
mag(refDir);
310 i < siteReferencePositions_.size();
314 vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1];
318 if (
mag(refDir & dir) < 1 - SMALL)
333 return siteReferencePositions_;
368 return pairPotentialSites_;
382 << sId <<
" site not found." 386 return pairPotentialSites_[
s];
393 return electrostaticSites_;
407 << sId <<
" site not found." 411 return electrostaticSites_[
s];
418 return momentOfInertia_;
424 return ((momentOfInertia_.
xx() < 0) && (momentOfInertia_.
yy() > 0));
430 return (momentOfInertia_.
zz() < 0);
459 return siteIds_.
size();
540 return sitePositions_;
546 return sitePositions_;
552 return specialPosition_;
558 return specialPosition_;
564 return potentialEnergy_;
570 return potentialEnergy_;
molecule(const polyMesh &mesh, const vector &position, const label celli, const label tetFacei, const label tetPtI, const tensor &Q, const vector &v, const vector &a, const vector &pi, const vector &tau, const vector &specialPosition, const constantProperties &constProps, const label special, const label id)
Construct from components.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
DiagTensor< scalar > diagTensor
A scalar version of the templated DiagTensor.
label degreesOfFreedom() const
Class to hold molecule constant properties.
void setSitePositions(const constantProperties &constProps)
const double e
Elementary charge.
bool pointMolecule() const
A list of keyword definitions, which are a keyword followed by any number of values (e...
const List< bool > & pairPotentialSites() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const List< vector > & siteForces() const
const tensor & rf() const
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void size(const label)
Override size to be inconsistent with allocated storage.
dimensionedVector eigenValues(const dimensionedTensor &dt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Vector< scalar > vector
A scalar version of the templated Vector.
particle(const polyMesh &mesh, const vector &position, const label celli, const label tetFacei, const label tetPtI)
Construct from components.
const fileName & name() const
Return the dictionary name.
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
const List< scalar > & siteCharges() const
const List< label > & siteIds() const
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Pre-declare SubField and related Field type.
const diagTensor & momentOfInertia() const
A class for handling words, derived from string.
const Field< vector > & siteReferencePositions() const
const vector & tau() const
errorManip< error > abort(error &err)
Templated 3D Vector derived from VectorSpace adding construction from 3 components, element access using x(), y() and z() member functions and the inner-product (dot-product) and cross product operators.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
dimensionedTensor eigenVectors(const dimensionedTensor &dt)
void setSize(const label)
Reset size of List.
scalar potentialEnergy() const
const vector & specialPosition() const
const vector & pi() const
bool pairPotentialSite(label sId) const
const List< bool > & electrostaticSites() const
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
bool linearMolecule() const
const List< vector > & sitePositions() const
Tensor< scalar > tensor
Tensor of scalars.
bool electrostaticSite(label sId) const
const List< scalar > & siteMasses() const
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.