binaryTree.C
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) 2016-2022 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 #include "binaryTree.H"
27 #include "SortableList.H"
28 
29 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
30 
31 bool Foam::binaryTree::inSubTree
32 (
33  const scalarField& phiq,
34  binaryNode* y,
35  chemPointISAT* x
36 )
37 {
38  if ((n2ndSearch_ < max2ndSearch_) && (y!=nullptr))
39  {
40  scalar vPhi = 0;
41  const scalarField& v = y->v();
42  const scalar a = y->a();
43  // compute v*phi
44  for (label i=0; i<phiq.size(); i++)
45  {
46  vPhi += phiq[i]*v[i];
47  }
48  if (vPhi <= a)// on the left side of the node
49  {
50  if (y->nodeLeft() == nullptr)// left is a chemPoint
51  {
52  n2ndSearch_++;
53  if (y->leafLeft()->inEOA(phiq))
54  {
55  x = y->leafLeft();
56  return true;
57  }
58  }
59  else // the left side is a node
60  {
61  if (inSubTree(phiq, y->nodeLeft(),x))
62  {
63  return true;
64  }
65  }
66 
67  // not on the left side, try the right side
68  if ((n2ndSearch_ < max2ndSearch_) && y->nodeRight() == nullptr)
69  {
70  n2ndSearch_++;
71  // we reach the end of the subTree we can return the result
72  if (y->leafRight()->inEOA(phiq))
73  {
74  x = y->leafRight();
75  return true;
76  }
77  else
78  {
79  x = nullptr;
80  return false;
81  }
82  }
83  else // test for n2ndSearch is done in the call of inSubTree
84  {
85  return inSubTree(phiq, y->nodeRight(),x);
86  }
87  }
88  else // on right side (symmetric of above)
89  {
90  if (y->nodeRight() == nullptr)
91  {
92  n2ndSearch_++;
93  if (y->leafRight()->inEOA(phiq))
94  {
95  return true;
96  }
97  }
98  else // the right side is a node
99  {
100  if (inSubTree(phiq, y->nodeRight(),x))
101  {
102  x = y->leafRight();
103  return true;
104  }
105  }
106  // if we reach this point, the retrieve has
107  // failed on the right side, explore the left side
108  if ((n2ndSearch_ < max2ndSearch_) && y->nodeLeft() == nullptr)
109  {
110  n2ndSearch_++;
111  if (y->leafLeft()->inEOA(phiq))
112  {
113  x = y->leafLeft();
114  return true;
115  }
116  else
117  {
118  x = nullptr;
119  return false;
120  }
121  }
122  else
123  {
124  return inSubTree(phiq, y->nodeLeft(),x);
125  }
126  }
127  }
128  else
129  {
130  return false;
131  }
132 }
133 
134 
135 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
136 
138 (
140  dictionary coeffsDict
141 )
142 :
143  table_(table),
144  root_(nullptr),
145  maxNLeafs_(coeffsDict.lookup<label>("maxNLeafs")),
146  size_(0),
147  n2ndSearch_(0),
148  max2ndSearch_(coeffsDict.lookupOrDefault("max2ndSearch",0)),
149  coeffsDict_(coeffsDict)
150 {}
151 
152 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
153 
155 (
156  const scalarField& phiq,
157  const scalarField& Rphiq,
158  const scalarSquareMatrix& A,
159  const scalarField& scaleFactor,
160  const scalar& epsTol,
161  const label nCols,
162  const label nActive,
164 )
165 {
166  if (size_ == 0) // no points are stored
167  {
168  // create an empty binary node and point root_ to it
169  root_ = new binaryNode();
170  // create the new chemPoint which holds the composition point
171  // phiq and the data to initialise the EOA
172  chemPointISAT* newChemPoint =
173  new chemPointISAT
174  (
175  table_,
176  phiq,
177  Rphiq,
178  A,
179  scaleFactor,
180  epsTol,
181  nCols,
182  nActive,
183  coeffsDict_,
184  root_
185  );
186  root_->leafLeft() = newChemPoint;
187  }
188  else // at least one point stored
189  {
190  // no reference chemPoint, a BT search is required
191  if (phi0 == nullptr)
192  {
193  binaryTreeSearch(phiq, root_,phi0);
194  }
195  // access to the parent node of the chemPoint
196  binaryNode* parentNode = phi0->node();
197 
198  // create the new chemPoint which holds the composition point
199  // phiq and the data to initialise the EOA
200  chemPointISAT* newChemPoint =
201  new chemPointISAT
202  (
203  table_,
204  phiq,
205  Rphiq,
206  A,
207  scaleFactor,
208  epsTol,
209  nCols,
210  nActive,
211  coeffsDict_
212  );
213  // insert new node on the parent node in the position of the
214  // previously stored leaf (phi0)
215  // the new node contains phi0 on the left and phiq on the right
216  // the hyper plane is computed in the binaryNode constructor
217  binaryNode* newNode;
218  if (size_>1)
219  {
220  newNode = new binaryNode(phi0, newChemPoint, parentNode);
221  // make the parent of phi0 point to the newly created node
222  insertNode(phi0, newNode);
223  }
224  else // size_ == 1 (because not equal to 0)
225  {
226  // when size is 1, the binaryNode is without hyperplane
227  deleteDemandDrivenData(root_);
228  newNode = new binaryNode(phi0, newChemPoint, nullptr);
229  root_ = newNode;
230  }
231 
232  phi0->node() = newNode;
233  newChemPoint->node()=newNode;
234  }
235  size_++;
236 }
237 
238 
240 (
241  const scalarField& phiq,
242  chemPointISAT*& x
243 )
244 {
245  // initialise n2ndSearch_
246  n2ndSearch_ = 0;
247  if ((n2ndSearch_ < max2ndSearch_) && (size_ > 1))
248  {
249  chemPointISAT* xS = chemPSibling(x);
250  if (xS != nullptr)
251  {
252  n2ndSearch_++;
253  if (xS->inEOA(phiq))
254  {
255  x = xS;
256  return true;
257  }
258  }
259  else if (inSubTree(phiq, nodeSibling(x),x))
260  {
261  return true;
262  }
263  // if we reach this point, no leafs were found at this depth or lower
264  // we move upward in the tree
265  binaryNode* y = x->node();
266  while((y->parent()!= nullptr) && (n2ndSearch_ < max2ndSearch_))
267  {
268  xS = chemPSibling(y);
269  if (xS != nullptr)
270  {
271  n2ndSearch_++;
272  if (xS->inEOA(phiq))
273  {
274  x=xS;
275  return true;
276  }
277  }
278  else if (inSubTree(phiq, nodeSibling(y),x))
279  {
280  return true;
281  }
282  y = y->parent();
283  }
284  // if we reach this point it is either because
285  // we did not find another covering EOA in the entire tree or
286  // we reach the maximum number of secondary search
287  return false;
288  }
289  else
290  {
291  return false;
292  }
293 }
294 
295 
297 {
298  if (size_ == 1) // only one point is stored
299  {
301  deleteDemandDrivenData(root_);
302  }
303  else if (size_ > 1)
304  {
305  binaryNode* z = phi0->node();
306  binaryNode* x;
307  chemPointISAT* siblingPhi0 = chemPSibling(phi0);
308 
309  if (siblingPhi0 != nullptr)// the sibling of phi0 is a chemPoint
310  {
311  // z was root (only two chemPoints in the tree)
312  if (z->parent() == nullptr)
313  {
314  root_ = new binaryNode();
315  root_->leafLeft()=siblingPhi0;
316  siblingPhi0->node()=root_;
317  }
318  else if (z == z->parent()->nodeLeft())
319  {
320  z->parent()->leafLeft() = siblingPhi0;
321  z->parent()->nodeLeft() = nullptr;
322  siblingPhi0->node() = z->parent();
323  }
324  else if (z == z->parent()->nodeRight())
325  {
326  z->parent()->leafRight() = siblingPhi0;
327  z->parent()->nodeRight() = nullptr;
328  siblingPhi0->node() = z->parent();
329  }
330  else
331  {
333  << "wrong addressing of the initial leaf"
334  << exit(FatalError);
335  }
336  }
337  else
338  {
339  x = nodeSibling(phi0);
340  if (x !=nullptr)
341  {
342  transplant(z, x);
343  }
344  else
345  {
347  << "inconsistent structure of the tree, no leaf and no node"
348  << exit(FatalError);
349  }
350  }
353  }
354  size_--;
355 }
356 
357 
359 {
360  //1) walk through the entire tree by starting with the tree's most left
361  // chemPoint
362  chemPointISAT* x = treeMin();
363  List<chemPointISAT*> chemPoints(size_);
364  label chemPointISATi=0;
365  while (x != nullptr)
366  {
367  chemPoints[chemPointISATi++] = x;
368  x = treeSuccessor(x);
369  }
370 
371  //2) compute the mean composition
372  scalarField mean(treeMin()->phi().size(), Zero);
373  forAll(chemPoints, j)
374  {
375  const scalarField& phij = chemPoints[j]->phi();
376  mean += phij;
377  }
378  mean /= size_;
379 
380  //3) compute the variance for each space direction
381  List<scalar> variance(treeMin()->phi().size(), Zero);
382  forAll(chemPoints, j)
383  {
384  const scalarField& phij = chemPoints[j]->phi();
385  forAll(variance, vi)
386  {
387  variance[vi] += sqr(phij[vi]-mean[vi]);
388  }
389  }
390 
391  //4) analyze what is the direction of the maximal variance
392  scalar maxVariance(-1.0);
393  label maxDir(-1);
394  forAll(variance, vi)
395  {
396  if (maxVariance < variance[vi])
397  {
398  maxVariance = variance[vi];
399  maxDir = vi;
400  }
401  }
402 
403  // maxDir indicates the direction of maximum variance
404  // we create the new root node by taking the two extreme points
405  // in this direction if these extreme points were not deleted in the
406  // cleaning that come before the balance function they are still important
407  // and the tree should therefore take them into account
408  SortableList<scalar> phiMaxDir(chemPoints.size(),0.0);
409  forAll(chemPoints, j)
410  {
411  phiMaxDir[j] = chemPoints[j]->phi()[maxDir];
412  }
413 
414  phiMaxDir.sort();
415  // delete reference to all node since the tree is reshaped
416  deleteAllNode();
417  root_ = nullptr;
418 
419  // add the node for the two extremum
420  binaryNode* newNode = new binaryNode
421  (
422  chemPoints[phiMaxDir.indices()[0]],
423  chemPoints[phiMaxDir.indices()[phiMaxDir.size()-1]],
424  nullptr
425  );
426  root_ = newNode;
427 
428  chemPoints[phiMaxDir.indices()[0]]->node() = newNode;
429  chemPoints[phiMaxDir.indices()[phiMaxDir.size()-1]]->node() = newNode;
430 
431  for (label cpi=1; cpi<chemPoints.size()-1; cpi++)
432  {
434  binaryTreeSearch
435  (
436  chemPoints[phiMaxDir.indices()[cpi]]->phi(),
437  root_,
438  phi0
439  );
440  // add the chemPoint
441  binaryNode* nodeToAdd = new binaryNode
442  (
443  phi0,
444  chemPoints[phiMaxDir.indices()[cpi]],
445  phi0->node()
446  );
447 
448  // make the parent of phi0 point to the newly created node
449  insertNode(phi0, nodeToAdd);
450  phi0->node() = nodeToAdd;
451  chemPoints[phiMaxDir.indices()[cpi]]->node() = nodeToAdd;
452  }
453 }
454 
455 
457 {
458  if (size_>1)
459  {
460  if (x == x->node()->leafLeft())
461  {
462  if (x->node()->nodeRight() == nullptr)
463  {
464  return x->node()->leafRight();
465  }
466  else
467  {
468  return treeMin(x->node()->nodeRight());
469  }
470  }
471  else if (x == x->node()->leafRight())
472  {
473  binaryNode* y = x->node();
474  while((y->parent() !=nullptr))
475  {
476  if (y == y->parent()->nodeLeft())
477  {
478  if (y->parent()->nodeRight() == nullptr)
479  {
480  return y->parent()->leafRight();
481  }
482  else
483  {
484  return treeMin(y->parent()->nodeRight());
485  }
486  }
487  y=y->parent();
488  }
489  // when we reach this point, y points to the root and
490  // never entered in the if loop (coming from the right)
491  // so we are at the tree maximum and there is no successor
492  return nullptr;
493  }
494  else
495  {
497  << "inconsistent structure of the tree, no leaf and no node"
498  << exit(FatalError);
499  return nullptr;
500  }
501  }
502 
503  return nullptr;
504 }
505 
506 
507 // ************************************************************************* //
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
scalar y
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:91
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: SortableList.H:55
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
Definition: SortableList.H:96
void sort()
(stable) sort the list (if changed after construction time)
Definition: SortableList.C:112
Node of the binary tree.
Definition: binaryNode.H:49
binaryNode *& parent()
Definition: binaryNode.H:158
chemPointISAT *& leafLeft()
Access.
Definition: binaryNode.H:138
binaryNode *& nodeRight()
Definition: binaryNode.H:153
chemPointISAT *& leafRight()
Definition: binaryNode.H:143
binaryNode *& nodeLeft()
Definition: binaryNode.H:148
void deleteLeaf(chemPointISAT *&phi0)
Delete a leaf from the binary tree and reshape the binary tree for.
Definition: binaryTree.C:296
bool secondaryBTSearch(const scalarField &phiq, chemPointISAT *&x)
Definition: binaryTree.C:240
chemPointISAT * treeSuccessor(chemPointISAT *x)
Definition: binaryTree.C:456
void insertNewLeaf(const scalarField &phiq, const scalarField &Rphiq, const scalarSquareMatrix &A, const scalarField &scaleFactor, const scalar &epsTol, const label nCols, const label nActive, chemPointISAT *&phi0)
Definition: binaryTree.C:155
binaryTree(chemistryTabulationMethods::ISAT &table, dictionary coeffsDict)
Constructors.
Definition: binaryTree.C:138
void balance()
Cheap balance function.
Definition: binaryTree.C:358
Leaf of the binary tree. The chemPoint stores the composition 'phi', the mapping of this composition ...
bool inEOA(const scalarField &phiq)
To RETRIEVE the mapping from the stored chemPoint phi, the query.
binaryNode *& node()
Implementation of the ISAT (In-situ adaptive tabulation), for chemistry calculation.
Definition: ISAT.H:63
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
volScalarField scalarField(fieldObject, mesh)
const dimensionedScalar phi0
Magnetic flux quantum: default SI units: [Wb].
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static const zero Zero
Definition: zero.H:97
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
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void deleteDemandDrivenData(DataType *&dataPtr)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
error FatalError