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 | Copyright (C) 2011-2016 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  setInteracionSiteBools
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 vector& position,
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, position, 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 
257 // * * * * * * * constantProperties Private Member Functions * * * * * * * * //
258 
259 inline void Foam::molecule::constantProperties::checkSiteListSizes() const
260 {
261  if
262  (
263  siteIds_.size() != siteReferencePositions_.size()
264  || siteIds_.size() != siteCharges_.size()
265  )
266  {
268  << "Sizes of site id, charge and "
269  << "referencePositions are not the same. "
270  << nl << abort(FatalError);
271  }
272 }
273 
274 
275 inline void Foam::molecule::constantProperties::setInteracionSiteBools
276 (
277  const List<word>& siteIds,
278  const List<word>& pairPotSiteIds
279 )
280 {
281  pairPotentialSites_.setSize(siteIds_.size());
282 
283  electrostaticSites_.setSize(siteIds_.size());
284 
285  forAll(siteIds_, i)
286  {
287  const word& id(siteIds[i]);
288 
289  pairPotentialSites_[i] = (findIndex(pairPotSiteIds, id) > -1);
290 
291  electrostaticSites_[i] = (mag(siteCharges_[i]) > VSMALL);
292  }
293 }
294 
295 
296 inline bool Foam::molecule::constantProperties::linearMoleculeTest() const
297 {
298  if (siteIds_.size() == 2)
299  {
300  return true;
301  }
302 
303  vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0];
304 
305  refDir /= mag(refDir);
306 
307  for
308  (
309  label i = 2;
310  i < siteReferencePositions_.size();
311  i++
312  )
313  {
314  vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1];
315 
316  dir /= mag(dir);
317 
318  if (mag(refDir & dir) < 1 - SMALL)
319  {
320  return false;
321  }
322  }
323 
324  return true;
325 }
326 
327 
328 // * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
329 
330 inline const Foam::Field<Foam::vector>&
332 {
333  return siteReferencePositions_;
334 }
335 
336 
337 inline const Foam::List<Foam::scalar>&
339 {
340  return siteMasses_;
341 }
342 
343 
344 inline const Foam::List<Foam::scalar>&
346 {
347  return siteCharges_;
348 }
349 
350 
351 inline const Foam::List<Foam::label>&
353 {
354  return siteIds_;
355 }
356 
357 
360 {
361  return siteIds_;
362 }
363 
364 
365 inline const Foam::List<bool>&
367 {
368  return pairPotentialSites_;
369 }
370 
371 
373 (
374  label sId
375 ) const
376 {
377  label s = findIndex(siteIds_, sId);
378 
379  if (s == -1)
380  {
382  << sId << " site not found."
383  << nl << abort(FatalError);
384  }
385 
386  return pairPotentialSites_[s];
387 }
388 
389 
390 inline const Foam::List<bool>&
392 {
393  return electrostaticSites_;
394 }
395 
396 
398 (
399  label sId
400 ) const
401 {
402  label s = findIndex(siteIds_, sId);
403 
404  if (s == -1)
405  {
407  << sId << " site not found."
408  << nl << abort(FatalError);
409  }
410 
411  return electrostaticSites_[s];
412 }
413 
414 
415 inline const Foam::diagTensor&
417 {
418  return momentOfInertia_;
419 }
420 
421 
423 {
424  return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
425 }
426 
427 
429 {
430  return (momentOfInertia_.zz() < 0);
431 }
432 
433 
435 {
436  if (linearMolecule())
437  {
438  return 5;
439  }
440  else if (pointMolecule())
441  {
442  return 3;
443  }
444  else
445  {
446  return 6;
447  }
448 }
449 
450 
451 inline Foam::scalar Foam::molecule::constantProperties::mass() const
452 {
453  return mass_;
454 }
455 
456 
458 {
459  return siteIds_.size();
460 
461 }
462 
463 
464 // * * * * * * * * * * * * molecule Member Functions * * * * * * * * * * * * //
465 
466 inline const Foam::tensor& Foam::molecule::Q() const
467 {
468  return Q_;
469 }
470 
471 
473 {
474  return Q_;
475 }
476 
477 
478 inline const Foam::vector& Foam::molecule::v() const
479 {
480  return v_;
481 }
482 
483 
485 {
486  return v_;
487 }
488 
489 
490 inline const Foam::vector& Foam::molecule::a() const
491 {
492  return a_;
493 }
494 
495 
497 {
498  return a_;
499 }
500 
501 
502 inline const Foam::vector& Foam::molecule::pi() const
503 {
504  return pi_;
505 }
506 
507 
509 {
510  return pi_;
511 }
512 
513 
514 inline const Foam::vector& Foam::molecule::tau() const
515 {
516  return tau_;
517 }
518 
519 
521 {
522  return tau_;
523 }
524 
525 
527 {
528  return siteForces_;
529 }
530 
531 
533 {
534  return siteForces_;
535 }
536 
537 
539 {
540  return sitePositions_;
541 }
542 
543 
545 {
546  return sitePositions_;
547 }
548 
549 
551 {
552  return specialPosition_;
553 }
554 
555 
557 {
558  return specialPosition_;
559 }
560 
561 
562 inline Foam::scalar Foam::molecule::potentialEnergy() const
563 {
564  return potentialEnergy_;
565 }
566 
567 
568 inline Foam::scalar& Foam::molecule::potentialEnergy()
569 {
570  return potentialEnergy_;
571 }
572 
573 
574 inline const Foam::tensor& Foam::molecule::rf() const
575 {
576  return rf_;
577 }
578 
579 
581 {
582  return rf_;
583 }
584 
585 
587 {
588  return special_;
589 }
590 
591 
592 inline bool Foam::molecule::tethered() const
593 {
594  return special_ == SPECIAL_TETHERED;
595 }
596 
597 
599 {
600  return id_;
601 }
602 
603 
604 // ************************************************************************* //
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.
Definition: moleculeI.H:221
const tensor & Q() const
Definition: moleculeI.H:466
const Cmpt & z() const
Definition: VectorI.H:87
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
const Cmpt & yy() const
Definition: TensorI.H:181
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
DiagTensor< scalar > diagTensor
A scalar version of the templated DiagTensor.
Definition: diagTensor.H:48
Class to hold molecule constant properties.
Definition: molecule.H:89
void setSitePositions(const constantProperties &constProps)
Definition: molecule.C:227
const Cmpt & x() const
Definition: VectorI.H:75
const double e
Elementary charge.
Definition: doubleFloat.H:78
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
const List< bool > & pairPotentialSites() const
Definition: moleculeI.H:366
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const List< vector > & siteForces() const
Definition: moleculeI.H:526
const Cmpt & zz() const
Definition: DiagTensorI.H:90
const tensor & rf() const
Definition: moleculeI.H:574
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
const Cmpt & xx() const
Definition: DiagTensorI.H:78
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
dimensionedVector eigenValues(const dimensionedTensor &dt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
particle(const polyMesh &mesh, const vector &position, const label celli, const label tetFacei, const label tetPtI)
Construct from components.
Definition: particle.C:48
const fileName & name() const
Return the dictionary name.
Definition: dictionary.H:103
Vector< Cmpt > x() const
Definition: TensorI.H:279
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
Definition: transform.H:47
const List< scalar > & siteCharges() const
Definition: moleculeI.H:345
Vector< Cmpt > y() const
Definition: TensorI.H:286
const List< label > & siteIds() const
Definition: moleculeI.H:352
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.
Definition: Field.H:57
const diagTensor & momentOfInertia() const
Definition: moleculeI.H:416
A class for handling words, derived from string.
Definition: word.H:59
const Field< vector > & siteReferencePositions() const
Definition: moleculeI.H:331
const vector & v() const
Definition: moleculeI.H:478
const vector & tau() const
Definition: moleculeI.H:514
const Cmpt & xx() const
Definition: TensorI.H:153
const Cmpt & y() const
Definition: VectorI.H:81
static const zero Zero
Definition: zero.H:91
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Vector< Cmpt > z() const
Definition: TensorI.H:293
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
static const char nl
Definition: Ostream.H:262
bool tethered() const
Definition: moleculeI.H:592
const Cmpt & zz() const
Definition: TensorI.H:209
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.
Definition: List.C:295
scalar potentialEnergy() const
Definition: moleculeI.H:562
const vector & specialPosition() const
Definition: moleculeI.H:550
const vector & pi() const
Definition: moleculeI.H:502
label id() const
Definition: moleculeI.H:598
bool pairPotentialSite(label sId) const
Definition: moleculeI.H:373
const List< bool > & electrostaticSites() const
Definition: moleculeI.H:391
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
label special() const
Definition: moleculeI.H:586
volScalarField & p
const List< vector > & sitePositions() const
Definition: moleculeI.H:538
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
bool electrostaticSite(label sId) const
Definition: moleculeI.H:398
const List< scalar > & siteMasses() const
Definition: moleculeI.H:338
const vector & a() const
Definition: moleculeI.H:490
const Cmpt & yy() const
Definition: DiagTensorI.H:84
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451