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