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, coordinates, celli, tetFacei, tetPti),
245 specialPosition_(specialPosition),
246 potentialEnergy_(0.0),
251 sitePositions_(constProps.
nSites())
267 const vector& specialPosition,
280 specialPosition_(specialPosition),
281 potentialEnergy_(0.0),
286 sitePositions_(constProps.
nSites())
294 inline void Foam::molecule::constantProperties::checkSiteListSizes()
const 298 siteIds_.
size() != siteReferencePositions_.
size()
299 || siteIds_.
size() != siteCharges_.
size()
303 <<
"Sizes of site id, charge and " 304 <<
"referencePositions are not the same. " 310 inline void Foam::molecule::constantProperties::setInteracionSiteBools
322 const word&
id(siteIds[i]);
324 pairPotentialSites_[i] = (
findIndex(pairPotSiteIds,
id) > -1);
326 electrostaticSites_[i] = (
mag(siteCharges_[i]) > vSmall);
331 inline bool Foam::molecule::constantProperties::linearMoleculeTest()
const 333 if (siteIds_.
size() == 2)
338 vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0];
340 refDir /=
mag(refDir);
345 i < siteReferencePositions_.size();
349 vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1];
353 if (
mag(refDir & dir) < 1 - small)
368 return siteReferencePositions_;
403 return pairPotentialSites_;
417 << sId <<
" site not found." 421 return pairPotentialSites_[
s];
428 return electrostaticSites_;
442 << sId <<
" site not found." 446 return electrostaticSites_[
s];
453 return momentOfInertia_;
459 return ((momentOfInertia_.
xx() < 0) && (momentOfInertia_.
yy() > 0));
465 return (momentOfInertia_.
zz() < 0);
494 return siteIds_.
size();
575 return sitePositions_;
581 return sitePositions_;
587 return specialPosition_;
593 return specialPosition_;
599 return potentialEnergy_;
605 return potentialEnergy_;
scalar potentialEnergy() const
const List< bool > & electrostaticSites() const
#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.
const List< vector > & sitePositions() const
DiagTensor< scalar > diagTensor
A scalar version of the templated DiagTensor.
const tensor & rf() const
const diagTensor & momentOfInertia() const
Class to hold molecule constant properties.
void setSitePositions(const constantProperties &constProps)
A list of keyword definitions, which are a keyword followed by any number of values (e...
const List< vector > & siteForces() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
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)
const List< label > & siteIds() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
const vector & pi() const
Vector< scalar > vector
A scalar version of the templated Vector.
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
const Field< vector > & siteReferencePositions() const
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))
const fileName & name() const
Return the dictionary name.
molecule(const polyMesh &mesh, const barycentric &coordinates, 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.
Pre-declare SubField and related Field type.
A class for handling words, derived from string.
particle(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from components.
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.
const List< bool > & pairPotentialSites() const
bool pointMolecule() const
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
dimensionedTensor eigenVectors(const dimensionedTensor &dt)
void setSize(const label)
Reset size of List.
bool electrostaticSite(label sId) const
label degreesOfFreedom() const
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Mesh consisting of general polyhedral cells.
bool linearMolecule() const
bool pairPotentialSite(label sId) const
const List< scalar > & siteMasses() const
const vector & tau() const
Tensor< scalar > tensor
Tensor of scalars.
const List< scalar > & siteCharges() const
const vector & specialPosition() const
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.