Go to the documentation of this file.
33 siteReferencePositions_(),
37 pairPotentialSites_(),
38 electrostaticSites_(),
39 momentOfInertia_(
Zero),
49 siteReferencePositions_(
dict.lookup(
"siteReferencePositions")),
50 siteMasses_(
dict.lookup(
"siteMasses")),
51 siteCharges_(
dict.lookup(
"siteCharges")),
53 pairPotentialSites_(),
54 electrostaticSites_(),
60 setInteracionSiteBools
66 mass_ =
sum(siteMasses_);
73 forAll(siteReferencePositions_, i)
75 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
78 centreOfMass /= mass_;
80 siteReferencePositions_ -= centreOfMass;
82 if (siteIds_.size() == 1)
86 siteReferencePositions_[0] =
Zero;
90 else if (linearMoleculeTest())
99 siteReferencePositions_[1] - siteReferencePositions_[0]
104 siteReferencePositions_ = (
Q & siteReferencePositions_);
111 forAll(siteReferencePositions_, i)
113 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
116 centreOfMass /= mass_;
118 siteReferencePositions_ -= centreOfMass;
122 forAll(siteReferencePositions_, i)
124 const vector&
p(siteReferencePositions_[i]);
145 forAll(siteReferencePositions_, i)
147 const vector&
p(siteReferencePositions_[i]);
151 p.y()*
p.y() +
p.z()*
p.z(), -
p.x()*
p.y(), -
p.x()*
p.z(),
152 p.x()*
p.x() +
p.z()*
p.z(), -
p.y()*
p.z(),
153 p.x()*
p.x() +
p.y()*
p.y()
160 <<
"An eigenvalue of the inertia tensor is zero, the molecule "
161 <<
dict.name() <<
" is not a valid 6DOF shape."
180 siteReferencePositions_ = (
Q & siteReferencePositions_);
189 forAll(siteReferencePositions_, i)
191 centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
194 centreOfMass /= mass_;
196 siteReferencePositions_ -= centreOfMass;
203 forAll(siteReferencePositions_, i)
205 const vector&
p(siteReferencePositions_[i]);
209 p.y()*
p.y() +
p.z()*
p.z(), -
p.x()*
p.y(), -
p.x()*
p.z(),
210 p.x()*
p.x() +
p.z()*
p.z(), -
p.y()*
p.z(),
211 p.x()*
p.x() +
p.y()*
p.y()
230 const label tetFacei,
251 potentialEnergy_(0.0),
256 sitePositions_(constProps.
nSites())
286 potentialEnergy_(0.0),
291 sitePositions_(constProps.
nSites())
299 inline void Foam::molecule::constantProperties::checkSiteListSizes()
const
303 siteIds_.size() != siteReferencePositions_.size()
304 || siteIds_.size() != siteCharges_.size()
308 <<
"Sizes of site id, charge and "
309 <<
"referencePositions are not the same. "
315 inline void Foam::molecule::constantProperties::setInteracionSiteBools
317 const List<word>& siteIds,
318 const List<word>& pairPotSiteIds
321 pairPotentialSites_.setSize(siteIds_.size());
323 electrostaticSites_.setSize(siteIds_.size());
327 const word&
id = siteIds[i];
329 pairPotentialSites_[i] = pairPotSiteIds.found(
id);
330 electrostaticSites_[i] = (
mag(siteCharges_[i]) > VSMALL);
335 inline bool Foam::molecule::constantProperties::linearMoleculeTest()
const
337 if (siteIds_.size() == 2)
345 siteReferencePositions_[1] - siteReferencePositions_[0]
351 i < siteReferencePositions_.size();
358 siteReferencePositions_[i] - siteReferencePositions_[i-1]
361 if (
mag(refDir & dir) < 1 - SMALL)
376 return siteReferencePositions_;
411 return pairPotentialSites_;
420 label
s = siteIds_.find(sId);
425 << sId <<
" site not found."
429 return pairPotentialSites_[
s];
436 return electrostaticSites_;
445 label
s = siteIds_.find(sId);
450 << sId <<
" site not found."
454 return electrostaticSites_[
s];
461 return momentOfInertia_;
467 return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
473 return (momentOfInertia_.zz() < 0);
479 if (linearMolecule())
483 else if (pointMolecule())
502 return siteIds_.size();
583 return sitePositions_;
589 return sitePositions_;
595 return specialPosition_;
601 return specialPosition_;
607 return potentialEnergy_;
613 return potentialEnergy_;
const List< label > & siteIds() const
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.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
static constexpr const zero Zero
Global zero (0)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
bool pairPotentialSite(label sId) const
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
const vector & specialPosition() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
const vector & tau() const
Mesh consisting of general polyhedral cells.
#define forAll(list, i)
Loop across all elements in list.
void setSitePositions(const constantProperties &constProps)
const List< scalar > & siteCharges() const
const barycentric & coordinates() const
Return current particle coordinates.
bool pointMolecule() const
DiagTensor< scalar > diagTensor
DiagTensor of scalars, i.e. DiagTensor<scalar>.
Generic templated field type.
const List< vector > & sitePositions() const
messageStream Info
Information stream (uses stdout - output is on the master only)
const diagTensor & momentOfInertia() const
const Field< vector > & siteReferencePositions() const
const tensor & rf() const
const List< bool > & pairPotentialSites() const
bool linearMolecule() const
const List< vector > & siteForces() const
const List< bool > & electrostaticSites() const
vector position() const
Return current particle position.
SymmTensor< scalar > symmTensor
SymmTensor of scalars, i.e. SymmTensor<scalar>.
const polyMesh & mesh() const
Return the mesh database.
label degreesOfFreedom() const
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
errorManip< error > abort(error &err)
const List< scalar > & siteMasses() const
Vector< scalar > vector
A scalar version of the templated Vector.
scalar potentialEnergy() const
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
const dimensionedScalar e
Elementary charge.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Class to hold molecule constant properties.
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
particle(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from components.
bool electrostaticSite(label sId) const
const vector & pi() const