FreeFOAM The Cross-Platform CFD Toolkit
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) 2008-2010 OpenCFD Ltd.
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(vector::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] = vector::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 = vector::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 = diagTensor::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(tensor::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  {
154  FatalErrorIn("molecule::constantProperties::constantProperties")
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 = vector::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 = tensor::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 Cloud<molecule>& c,
223  const vector& position,
224  const label celli,
225  const tensor& Q,
226  const vector& v,
227  const vector& a,
228  const vector& pi,
229  const vector& tau,
230  const vector& specialPosition,
231  const constantProperties& constProps,
232  const label special,
233  const label id
234 
235 )
236 :
237  Particle<molecule>(c, position, celli),
238  Q_(Q),
239  v_(v),
240  a_(a),
241  pi_(pi),
242  tau_(tau),
243  specialPosition_(specialPosition),
244  potentialEnergy_(0.0),
245  rf_(tensor::zero),
246  special_(special),
247  id_(id),
248  siteForces_(constProps.nSites(), vector::zero),
249  sitePositions_(constProps.nSites())
250 {
251  setSitePositions(constProps);
252 }
253 
254 
255 // * * * constantProperties Private Member Functions * * * * * * * * * * * * //
256 
257 
258 inline void Foam::molecule::constantProperties::checkSiteListSizes() const
259 {
260  if
261  (
262  siteIds_.size() != siteReferencePositions_.size()
263  || siteIds_.size() != siteCharges_.size()
264  )
265  {
266  FatalErrorIn("molecule::constantProperties::checkSiteListSizes")
267  << "Sizes of site id, charge and "
268  << "referencePositions are not the same. "
269  << nl << abort(FatalError);
270  }
271 }
272 
273 
274 inline void Foam::molecule::constantProperties::setInteracionSiteBools
275 (
276  const List<word>& siteIds,
277  const List<word>& pairPotSiteIds
278 )
279 {
280  pairPotentialSites_.setSize(siteIds_.size());
281 
282  electrostaticSites_.setSize(siteIds_.size());
283 
284  forAll(siteIds_, i)
285  {
286  const word& id(siteIds[i]);
287 
288  pairPotentialSites_[i] = (findIndex(pairPotSiteIds, id) > -1);
289 
290  electrostaticSites_[i] = (mag(siteCharges_[i]) > VSMALL);
291  }
292 }
293 
294 
295 inline bool Foam::molecule::constantProperties::linearMoleculeTest() const
296 {
297  if (siteIds_.size() == 2)
298  {
299  return true;
300  }
301 
302  vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0];
303 
304  refDir /= mag(refDir);
305 
306  for
307  (
308  label i = 2;
309  i < siteReferencePositions_.size();
310  i++
311  )
312  {
313  vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1];
314 
315  dir /= mag(dir);
316 
317  if (mag(refDir & dir) < 1 - SMALL)
318  {
319  return false;
320  }
321  }
322 
323  return true;
324 }
325 
326 
327 // * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
328 
329 inline const Foam::Field<Foam::vector>&
331 {
332  return siteReferencePositions_;
333 }
334 
335 
336 inline const Foam::List<Foam::scalar>&
338 {
339  return siteMasses_;
340 }
341 
342 
343 inline const Foam::List<Foam::scalar>&
345 {
346  return siteCharges_;
347 }
348 
349 
350 inline const Foam::List<Foam::label>&
352 {
353  return siteIds_;
354 }
355 
356 
359 {
360  return siteIds_;
361 }
362 
363 
364 inline const Foam::List<bool>&
366 {
367  return pairPotentialSites_;
368 }
369 
370 
372 (
373  label sId
374 ) const
375 {
376  label s = findIndex(siteIds_, sId);
377 
378  if (s == -1)
379  {
380  FatalErrorIn("moleculeI.H") << nl
381  << sId << " site not found."
382  << nl << abort(FatalError);
383  }
384 
385  return pairPotentialSites_[s];
386 }
387 
388 
389 inline const Foam::List<bool>&
391 {
392  return electrostaticSites_;
393 }
394 
396 (
397  label sId
398 ) const
399 {
400  label s = findIndex(siteIds_, sId);
401 
402  if (s == -1)
403  {
405  (
406  "molecule::constantProperties::electrostaticSite(label)"
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 // * * * * * * * * * * * trackData Member Functions * * * * * * * * * * * * //
464 
466 {
467  return molCloud_;
468 }
469 
470 
471 inline Foam::label Foam::molecule::trackData::part() const
472 {
473  return part_;
474 }
475 
476 
477 // * * * * * * * * * * * * molecule Member Functions * * * * * * * * * * * * //
478 
479 inline const Foam::tensor& Foam::molecule::Q() const
480 {
481  return Q_;
482 }
483 
484 
486 {
487  return Q_;
488 }
489 
490 
491 inline const Foam::vector& Foam::molecule::v() const
492 {
493  return v_;
494 }
495 
496 
498 {
499  return v_;
500 }
501 
502 
503 inline const Foam::vector& Foam::molecule::a() const
504 {
505  return a_;
506 }
507 
508 
510 {
511  return a_;
512 }
513 
514 
515 inline const Foam::vector& Foam::molecule::pi() const
516 {
517  return pi_;
518 }
519 
520 
522 {
523  return pi_;
524 }
525 
526 
527 inline const Foam::vector& Foam::molecule::tau() const
528 {
529  return tau_;
530 }
531 
532 
534 {
535  return tau_;
536 }
537 
538 
540 {
541  return siteForces_;
542 }
543 
544 
546 {
547  return siteForces_;
548 }
549 
550 
552 {
553  return sitePositions_;
554 }
555 
556 
558 {
559  return sitePositions_;
560 }
561 
562 
564 {
565  return specialPosition_;
566 }
567 
568 
570 {
571  return specialPosition_;
572 }
573 
574 
575 inline Foam::scalar Foam::molecule::potentialEnergy() const
576 {
577  return potentialEnergy_;
578 }
579 
580 
581 inline Foam::scalar& Foam::molecule::potentialEnergy()
582 {
583  return potentialEnergy_;
584 }
585 
586 
587 inline const Foam::tensor& Foam::molecule::rf() const
588 {
589  return rf_;
590 }
591 
592 
594 {
595  return rf_;
596 }
597 
598 
599 inline Foam::label Foam::molecule::special() const
600 {
601  return special_;
602 }
603 
604 
605 inline bool Foam::molecule::tethered() const
606 {
607  return special_ == SPECIAL_TETHERED;
608 }
609 
610 
611 inline Foam::label Foam::molecule::id() const
612 {
613  return id_;
614 }
615 
616 
617 // ************************ vim: set sw=4 sts=4 et: ************************ //