moleculeI.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
30 
32 :
33  siteReferencePositions_(),
34  siteMasses_(),
35  siteCharges_(),
36  siteIds_(),
37  pairPotentialSites_(),
38  electrostaticSites_(),
39  momentOfInertia_(Zero),
40  mass_(0)
41 {}
42 
43 
45 (
46  const dictionary& dict
47 )
48 :
49  siteReferencePositions_(dict.lookup("siteReferencePositions")),
50  siteMasses_(dict.lookup("siteMasses")),
51  siteCharges_(dict.lookup("siteCharges")),
52  siteIds_(List<word>(dict.lookup("siteIds")).size(), -1),
53  pairPotentialSites_(),
54  electrostaticSites_(),
55  momentOfInertia_(),
56  mass_()
57 {
58  checkSiteListSizes();
59 
60  setInteracionSiteBools
61  (
62  List<word>(dict.lookup("siteIds")),
63  List<word>(dict.lookup("pairPotentialSiteIds"))
64  );
65 
66  mass_ = sum(siteMasses_);
67 
68  vector centreOfMass(Zero);
69 
70  // Calculate the centre of mass of the body and subtract it from each
71  // position
72 
73  forAll(siteReferencePositions_, i)
74  {
75  centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
76  }
77 
78  centreOfMass /= mass_;
79 
80  siteReferencePositions_ -= centreOfMass;
81 
82  if (siteIds_.size() == 1)
83  {
84  // Single site molecule - no rotational motion.
85 
86  siteReferencePositions_[0] = Zero;
87 
88  momentOfInertia_ = diagTensor(-1, -1, -1);
89  }
90  else if (linearMoleculeTest())
91  {
92  // Linear molecule.
93 
94  Info<< nl << "Linear molecule." << endl;
95 
96  const vector dir =
98  (
99  siteReferencePositions_[1] - siteReferencePositions_[0]
100  );
101 
102  tensor Q = rotationTensor(dir, vector(1,0,0));
103 
104  siteReferencePositions_ = (Q & siteReferencePositions_);
105 
106  // The rotation was around the centre of mass but remove any
107  // components that have crept in due to floating point errors
108 
109  centreOfMass = Zero;
110 
111  forAll(siteReferencePositions_, i)
112  {
113  centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
114  }
115 
116  centreOfMass /= mass_;
117 
118  siteReferencePositions_ -= centreOfMass;
119 
120  diagTensor momOfInertia = Zero;
121 
122  forAll(siteReferencePositions_, i)
123  {
124  const vector& p(siteReferencePositions_[i]);
125 
126  momOfInertia +=
127  siteMasses_[i]*diagTensor(0, p.x()*p.x(), p.x()*p.x());
128  }
129 
130  momentOfInertia_ = diagTensor
131  (
132  -1,
133  momOfInertia.yy(),
134  momOfInertia.zz()
135  );
136  }
137  else
138  {
139  // Fully 6DOF molecule
140 
141  // Calculate the inertia tensor in the current orientation
142 
143  symmTensor momOfInertia(Zero);
144 
145  forAll(siteReferencePositions_, i)
146  {
147  const vector& p(siteReferencePositions_[i]);
148 
149  momOfInertia += siteMasses_[i]*symmTensor
150  (
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()
154  );
155  }
156 
157  if (eigenValues(momOfInertia).x() < VSMALL)
158  {
160  << "An eigenvalue of the inertia tensor is zero, the molecule "
161  << dict.name() << " is not a valid 6DOF shape."
162  << nl << abort(FatalError);
163  }
164 
165  // Normalise the inertia tensor magnitude to avoid SMALL numbers in the
166  // components causing problems
167 
168  momOfInertia /= eigenValues(momOfInertia).x();
169 
170  tensor e(eigenVectors(momOfInertia));
171 
172  // Calculate the transformation between the principle axes to the
173  // global axes
174 
175  tensor Q =
176  vector(1,0,0)*e.x() + vector(0,1,0)*e.y() + vector(0,0,1)*e.z();
177 
178  // Transform the site positions
179 
180  siteReferencePositions_ = (Q & siteReferencePositions_);
181 
182  // Recalculating the moment of inertia with the new site positions
183 
184  // The rotation was around the centre of mass but remove any
185  // components that have crept in due to floating point errors
186 
187  centreOfMass = Zero;
188 
189  forAll(siteReferencePositions_, i)
190  {
191  centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
192  }
193 
194  centreOfMass /= mass_;
195 
196  siteReferencePositions_ -= centreOfMass;
197 
198  // Calculate the moment of inertia in the principle component
199  // reference frame
200 
201  momOfInertia = Zero;
202 
203  forAll(siteReferencePositions_, i)
204  {
205  const vector& p(siteReferencePositions_[i]);
206 
207  momOfInertia += siteMasses_[i]*symmTensor
208  (
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()
212  );
213  }
214 
215  momentOfInertia_ = diagTensor
216  (
217  momOfInertia.xx(),
218  momOfInertia.yy(),
219  momOfInertia.zz()
220  );
221  }
222 }
223 
224 
226 (
227  const polyMesh& mesh,
228  const barycentric& coordinates,
229  const label celli,
230  const label tetFacei,
231  const label tetPti,
232  const tensor& Q,
233  const vector& v,
234  const vector& a,
235  const vector& pi,
236  const vector& tau,
237  const vector& specialPosition,
238  const constantProperties& constProps,
239  const label special,
240  const label id
241 
242 )
243 :
244  particle(mesh, coordinates, celli, tetFacei, tetPti),
245  Q_(Q),
246  v_(v),
247  a_(a),
248  pi_(pi),
249  tau_(tau),
250  specialPosition_(specialPosition),
251  potentialEnergy_(0.0),
252  rf_(Zero),
253  special_(special),
254  id_(id),
255  siteForces_(constProps.nSites(), Zero),
256  sitePositions_(constProps.nSites())
257 {
258  setSitePositions(constProps);
259 }
260 
261 
263 (
264  const polyMesh& mesh,
265  const vector& position,
266  const label celli,
267  const tensor& Q,
268  const vector& v,
269  const vector& a,
270  const vector& pi,
271  const vector& tau,
272  const vector& specialPosition,
273  const constantProperties& constProps,
274  const label special,
275  const label id
276 
277 )
278 :
279  particle(mesh, position, celli),
280  Q_(Q),
281  v_(v),
282  a_(a),
283  pi_(pi),
284  tau_(tau),
285  specialPosition_(specialPosition),
286  potentialEnergy_(0.0),
287  rf_(Zero),
288  special_(special),
289  id_(id),
290  siteForces_(constProps.nSites(), Zero),
291  sitePositions_(constProps.nSites())
292 {
293  setSitePositions(constProps);
294 }
295 
296 
297 // * * * * * * * constantProperties Private Member Functions * * * * * * * * //
298 
299 inline void Foam::molecule::constantProperties::checkSiteListSizes() const
300 {
301  if
302  (
303  siteIds_.size() != siteReferencePositions_.size()
304  || siteIds_.size() != siteCharges_.size()
305  )
306  {
308  << "Sizes of site id, charge and "
309  << "referencePositions are not the same. "
310  << nl << abort(FatalError);
311  }
312 }
313 
314 
315 inline void Foam::molecule::constantProperties::setInteracionSiteBools
316 (
317  const List<word>& siteIds,
318  const List<word>& pairPotSiteIds
319 )
320 {
321  pairPotentialSites_.setSize(siteIds_.size());
322 
323  electrostaticSites_.setSize(siteIds_.size());
324 
325  forAll(siteIds_, i)
326  {
327  const word& id = siteIds[i];
328 
329  pairPotentialSites_[i] = pairPotSiteIds.found(id);
330  electrostaticSites_[i] = (mag(siteCharges_[i]) > VSMALL);
331  }
332 }
333 
334 
335 inline bool Foam::molecule::constantProperties::linearMoleculeTest() const
336 {
337  if (siteIds_.size() == 2)
338  {
339  return true;
340  }
341 
342  const vector refDir =
343  normalised
344  (
345  siteReferencePositions_[1] - siteReferencePositions_[0]
346  );
347 
348  for
349  (
350  label i = 2;
351  i < siteReferencePositions_.size();
352  i++
353  )
354  {
355  const vector dir =
356  normalised
357  (
358  siteReferencePositions_[i] - siteReferencePositions_[i-1]
359  );
360 
361  if (mag(refDir & dir) < 1 - SMALL)
362  {
363  return false;
364  }
365  }
366 
367  return true;
368 }
369 
370 
371 // * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
372 
373 inline const Foam::Field<Foam::vector>&
375 {
376  return siteReferencePositions_;
377 }
378 
379 
380 inline const Foam::List<Foam::scalar>&
382 {
383  return siteMasses_;
384 }
385 
386 
387 inline const Foam::List<Foam::scalar>&
389 {
390  return siteCharges_;
391 }
392 
393 
394 inline const Foam::List<Foam::label>&
396 {
397  return siteIds_;
398 }
399 
400 
403 {
404  return siteIds_;
405 }
406 
407 
408 inline const Foam::List<bool>&
410 {
411  return pairPotentialSites_;
412 }
413 
414 
416 (
417  label sId
418 ) const
419 {
420  label s = siteIds_.find(sId);
421 
422  if (s == -1)
423  {
425  << sId << " site not found."
426  << nl << abort(FatalError);
427  }
428 
429  return pairPotentialSites_[s];
430 }
431 
432 
433 inline const Foam::List<bool>&
435 {
436  return electrostaticSites_;
437 }
438 
439 
441 (
442  label sId
443 ) const
444 {
445  label s = siteIds_.find(sId);
446 
447  if (s == -1)
448  {
450  << sId << " site not found."
451  << nl << abort(FatalError);
452  }
453 
454  return electrostaticSites_[s];
455 }
456 
457 
458 inline const Foam::diagTensor&
460 {
461  return momentOfInertia_;
462 }
463 
464 
466 {
467  return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
468 }
469 
470 
472 {
473  return (momentOfInertia_.zz() < 0);
474 }
475 
476 
478 {
479  if (linearMolecule())
480  {
481  return 5;
482  }
483  else if (pointMolecule())
484  {
485  return 3;
486  }
487  else
488  {
489  return 6;
490  }
491 }
492 
493 
494 inline Foam::scalar Foam::molecule::constantProperties::mass() const
495 {
496  return mass_;
497 }
498 
499 
501 {
502  return siteIds_.size();
503 
504 }
505 
506 
507 // * * * * * * * * * * * * molecule Member Functions * * * * * * * * * * * * //
508 
509 inline const Foam::tensor& Foam::molecule::Q() const
510 {
511  return Q_;
512 }
513 
514 
516 {
517  return Q_;
518 }
519 
520 
521 inline const Foam::vector& Foam::molecule::v() const
522 {
523  return v_;
524 }
525 
526 
528 {
529  return v_;
530 }
531 
532 
533 inline const Foam::vector& Foam::molecule::a() const
534 {
535  return a_;
536 }
537 
538 
540 {
541  return a_;
542 }
543 
544 
545 inline const Foam::vector& Foam::molecule::pi() const
546 {
547  return pi_;
548 }
549 
550 
552 {
553  return pi_;
554 }
555 
556 
557 inline const Foam::vector& Foam::molecule::tau() const
558 {
559  return tau_;
560 }
561 
562 
564 {
565  return tau_;
566 }
567 
568 
570 {
571  return siteForces_;
572 }
573 
574 
576 {
577  return siteForces_;
578 }
579 
580 
582 {
583  return sitePositions_;
584 }
585 
586 
588 {
589  return sitePositions_;
590 }
591 
592 
594 {
595  return specialPosition_;
596 }
597 
598 
600 {
601  return specialPosition_;
602 }
603 
604 
605 inline Foam::scalar Foam::molecule::potentialEnergy() const
606 {
607  return potentialEnergy_;
608 }
609 
610 
611 inline Foam::scalar& Foam::molecule::potentialEnergy()
612 {
613  return potentialEnergy_;
614 }
615 
616 
617 inline const Foam::tensor& Foam::molecule::rf() const
618 {
619  return rf_;
620 }
621 
622 
624 {
625  return rf_;
626 }
627 
628 
629 inline Foam::label Foam::molecule::special() const
630 {
631  return special_;
632 }
633 
634 
635 inline bool Foam::molecule::tethered() const
636 {
637  return special_ == SPECIAL_TETHERED;
638 }
639 
640 
641 inline Foam::label Foam::molecule::id() const
642 {
643  return id_;
644 }
645 
646 
647 // ************************************************************************* //
Foam::molecule::constantProperties::siteIds
const List< label > & siteIds() const
Definition: moleculeI.H:395
Foam::Tensor< scalar >
Foam::SymmTensor< scalar >
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::DiagTensor< scalar >
Foam::molecule::molecule
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.
Definition: moleculeI.H:226
s
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))
Definition: gmvOutputSpray.H:25
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::eigenValues
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
Definition: dimensionedTensor.C:149
Foam::molecule::constantProperties::pairPotentialSite
bool pairPotentialSite(label sId) const
Definition: moleculeI.H:416
Foam::eigenVectors
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
Definition: dimensionedTensor.C:160
Foam::molecule::specialPosition
const vector & specialPosition() const
Definition: moleculeI.H:593
Foam::molecule::id
label id() const
Definition: moleculeI.H:641
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::molecule::tau
const vector & tau() const
Definition: moleculeI.H:557
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::molecule::setSitePositions
void setSitePositions(const constantProperties &constProps)
Definition: molecule.C:226
Foam::molecule::constantProperties::siteCharges
const List< scalar > & siteCharges() const
Definition: moleculeI.H:388
Foam::particle::coordinates
const barycentric & coordinates() const
Return current particle coordinates.
Definition: particleI.H:143
Foam::SymmTensor::xx
const Cmpt & xx() const
Definition: SymmTensorI.H:84
Foam::molecule::constantProperties::pointMolecule
bool pointMolecule() const
Definition: moleculeI.H:471
Foam::molecule::special
label special() const
Definition: moleculeI.H:629
Foam::diagTensor
DiagTensor< scalar > diagTensor
DiagTensor of scalars, i.e. DiagTensor<scalar>.
Definition: diagTensor.H:53
Foam::Field
Generic templated field type.
Definition: Field.H:63
Foam::molecule::sitePositions
const List< vector > & sitePositions() const
Definition: moleculeI.H:581
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::molecule::constantProperties::momentOfInertia
const diagTensor & momentOfInertia() const
Definition: moleculeI.H:459
Foam::molecule::constantProperties::siteReferencePositions
const Field< vector > & siteReferencePositions() const
Definition: moleculeI.H:374
Foam::molecule::rf
const tensor & rf() const
Definition: moleculeI.H:617
Foam::molecule::constantProperties::pairPotentialSites
const List< bool > & pairPotentialSites() const
Definition: moleculeI.H:409
Foam::molecule::constantProperties::linearMolecule
bool linearMolecule() const
Definition: moleculeI.H:465
Foam::molecule::siteForces
const List< vector > & siteForces() const
Definition: moleculeI.H:569
Foam::molecule::constantProperties::electrostaticSites
const List< bool > & electrostaticSites() const
Definition: moleculeI.H:434
Foam::particle::position
vector position() const
Return current particle position.
Definition: particleI.H:314
Foam::Barycentric< scalar >
Foam::symmTensor
SymmTensor< scalar > symmTensor
SymmTensor of scalars, i.e. SymmTensor<scalar>.
Definition: symmTensor.H:59
Foam::particle::mesh
const polyMesh & mesh() const
Return the mesh database.
Definition: particleI.H:137
Foam::molecule::constantProperties::degreesOfFreedom
label degreesOfFreedom() const
Definition: moleculeI.H:477
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::molecule::constantProperties::nSites
label nSites() const
Definition: moleculeI.H:500
Foam::molecule::constantProperties::siteMasses
const List< scalar > & siteMasses() const
Definition: moleculeI.H:381
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::molecule::Q
const tensor & Q() const
Definition: moleculeI.H:509
Foam::molecule::potentialEnergy
scalar potentialEnergy() const
Definition: moleculeI.H:605
Foam::normalised
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:487
Foam::molecule::a
const vector & a() const
Definition: moleculeI.H:533
Foam::DiagTensor::zz
const Cmpt & zz() const
Definition: DiagTensorI.H:88
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::molecule::v
const vector & v() const
Definition: moleculeI.H:521
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::SymmTensor::yy
const Cmpt & yy() const
Definition: SymmTensorI.H:108
Foam::molecule::constantProperties::mass
scalar mass() const
Definition: moleculeI.H:494
Foam::Vector
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
Definition: Vector.H:62
Foam::List< word >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
Foam::molecule::constantProperties
Class to hold molecule constant properties.
Definition: molecule.H:91
Foam::SymmTensor::zz
const Cmpt & zz() const
Definition: SymmTensorI.H:132
Foam::DiagTensor::yy
const Cmpt & yy() const
Definition: DiagTensorI.H:82
Foam::molecule::SPECIAL_TETHERED
Definition: molecule.H:83
Foam::rotationTensor
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
Definition: transform.H:51
Foam::molecule::tethered
bool tethered() const
Definition: moleculeI.H:635
Foam::particle::particle
particle(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from components.
Definition: particle.C:516
Foam::molecule::constantProperties::constantProperties
constantProperties()
Definition: moleculeI.H:31
Foam::molecule::constantProperties::electrostaticSite
bool electrostaticSite(label sId) const
Definition: moleculeI.H:441
Foam::molecule::pi
const vector & pi() const
Definition: moleculeI.H:545