tetOverlapVolume.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) 2012-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 "tetOverlapVolume.H"
27 #include "polyMesh.H"
28 #include "OFstream.H"
29 #include "treeBoundBox.H"
30 #include "indexedOctree.H"
31 #include "treeDataCell.H"
32 #include "cutTriTet.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
39 }
40 
41 
42 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
43 
45 {}
46 
47 
48 // * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
49 
50 Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol
51 (
52  const tetPointRef& tetA,
53  const tetPointRef& tetB
54 ) const
55 {
56  typedef cutTriTet::noOp noOp;
57 
58  // A maximum of three cuts are made (the tets that result from the final cut
59  // are not stored), and each cut can create at most three tets. The
60  // temporary storage must therefore extend to 3^3 = 27 tets.
61  typedef cutTetList<27> tetListType;
62  static tetListType cutTetList1, cutTetList2;
63 
64  // face 0
65  const plane pl0(tetB.b(), tetB.d(), tetB.c());
66  if (!pl0.valid())
67  {
68  return 0;
69  }
70  const FixedList<point, 4> t({tetA.a(), tetA.b(), tetA.c(), tetA.d()});
71  cutTetList1.clear();
72  tetCut(t, pl0, cutTriTet::appendOp<tetListType>(cutTetList1), noOp());
73  if (cutTetList1.size() == 0)
74  {
75  return 0;
76  }
77 
78  // face 1
79  const plane pl1(tetB.a(), tetB.c(), tetB.d());
80  if (!pl1.valid())
81  {
82  return 0;
83  }
84  cutTetList2.clear();
85  for (label i = 0; i < cutTetList1.size(); i++)
86  {
87  const FixedList<point, 4>& t = cutTetList1[i];
88  tetCut(t, pl1, cutTriTet::appendOp<tetListType>(cutTetList2), noOp());
89  }
90  if (cutTetList2.size() == 0)
91  {
92  return 0;
93  }
94 
95  // face 2
96  const plane pl2(tetB.a(), tetB.d(), tetB.b());
97  if (!pl2.valid())
98  {
99  return 0;
100  }
101  cutTetList1.clear();
102  for (label i = 0; i < cutTetList2.size(); i++)
103  {
104  const FixedList<point, 4>& t = cutTetList2[i];
105  tetCut(t, pl2, cutTriTet::appendOp<tetListType>(cutTetList1), noOp());
106  }
107  if (cutTetList1.size() == 0)
108  {
109  return 0;
110  }
111 
112  // face 3
113  const plane pl3(tetB.a(), tetB.b(), tetB.c());
114  if (!pl3.valid())
115  {
116  return 0;
117  }
118  scalar v = 0;
119  for (label i = 0; i < cutTetList1.size(); i++)
120  {
121  const FixedList<point, 4>& t = cutTetList1[i];
122  v += tetCut(t, pl3, cutTriTet::volumeOp(), noOp());
123  }
124 
125  return v;
126 }
127 
128 
129 Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
130 (
131  const pointField& points,
132  const face& f,
133  const point& fc
134 ) const
135 {
136  treeBoundBox bb(fc, fc);
137  forAll(f, fp)
138  {
139  const point& pt = points[f[fp]];
140  bb.min() = min(bb.min(), pt);
141  bb.max() = max(bb.max(), pt);
142  }
143  return bb;
144 }
145 
146 
147 // * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * * //
148 
150 (
151  const polyMesh& fromMesh,
152  const polyMesh& toMesh,
153  const label iTo
154 ) const
155 {
156  const indexedOctree<treeDataCell>& treeA = fromMesh.cellTree();
157 
158  treeBoundBox bbB(toMesh.points(), toMesh.cellPoints()[iTo]);
159 
160  return treeA.findBox(bbB);
161 }
162 
163 
165 (
166  const primitiveMesh& mesh,
167  const label cellI
168 ) const
169 {
170  const cell& cFaces = mesh.cells()[cellI];
171  const point& cc = mesh.cellCentres()[cellI];
172 
173  scalar vol = 0.0;
174 
175  forAll(cFaces, cF)
176  {
177  label faceI = cFaces[cF];
178 
179  const face& f = mesh.faces()[faceI];
180 
181  bool own = (mesh.faceOwner()[faceI] == cellI);
182 
183  label tetBasePtI = 0;
184 
185  const point& tetBasePt = mesh.points()[f[tetBasePtI]];
186 
187  for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
188  {
189  label facePtI = (tetPtI + tetBasePtI) % f.size();
190  label otherFacePtI = f.fcIndex(facePtI);
191 
192  label pt0I = -1;
193  label pt1I = -1;
194 
195  if (own)
196  {
197  pt0I = f[facePtI];
198  pt1I = f[otherFacePtI];
199  }
200  else
201  {
202  pt0I = f[otherFacePtI];
203  pt1I = f[facePtI];
204  }
205 
206  const tetPointRef tet
207  (
208  cc,
209  tetBasePt,
210  mesh.points()[pt0I],
211  mesh.points()[pt1I]
212  );
213 
214  vol += tet.mag();
215  }
216  }
217 
218  return vol;
219 }
220 
221 
223 (
224  const primitiveMesh& meshA,
225  const label cellAI,
226  const primitiveMesh& meshB,
227  const label cellBI,
228  const treeBoundBox& cellBbB,
229  const scalar threshold
230 ) const
231 {
232  const cell& cFacesA = meshA.cells()[cellAI];
233  const point& ccA = meshA.cellCentres()[cellAI];
234 
235  const cell& cFacesB = meshB.cells()[cellBI];
236  const point& ccB = meshB.cellCentres()[cellBI];
237 
238  scalar vol = 0.0;
239 
240  forAll(cFacesA, cFA)
241  {
242  label faceAI = cFacesA[cFA];
243 
244  const face& fA = meshA.faces()[faceAI];
245  const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
246  if (!pyrA.overlaps(cellBbB))
247  {
248  continue;
249  }
250 
251  bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
252 
253  label tetBasePtAI = 0;
254 
255  const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
256 
257  for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
258  {
259  label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
260  label otherFacePtAI = fA.fcIndex(facePtAI);
261 
262  label pt0I = -1;
263  label pt1I = -1;
264 
265  if (ownA)
266  {
267  pt0I = fA[facePtAI];
268  pt1I = fA[otherFacePtAI];
269  }
270  else
271  {
272  pt0I = fA[otherFacePtAI];
273  pt1I = fA[facePtAI];
274  }
275 
276  const tetPointRef tetA
277  (
278  ccA,
279  tetBasePtA,
280  meshA.points()[pt0I],
281  meshA.points()[pt1I]
282  );
283  const treeBoundBox tetABb(tetA.bounds());
284 
285 
286  // Loop over tets of cellB
287  forAll(cFacesB, cFB)
288  {
289  label faceBI = cFacesB[cFB];
290 
291  const face& fB = meshB.faces()[faceBI];
292  const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
293  if (!pyrB.overlaps(pyrA))
294  {
295  continue;
296  }
297 
298  bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
299 
300  label tetBasePtBI = 0;
301 
302  const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
303 
304  for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
305  {
306  label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
307  label otherFacePtBI = fB.fcIndex(facePtBI);
308 
309  label pt0I = -1;
310  label pt1I = -1;
311 
312  if (ownB)
313  {
314  pt0I = fB[facePtBI];
315  pt1I = fB[otherFacePtBI];
316  }
317  else
318  {
319  pt0I = fB[otherFacePtBI];
320  pt1I = fB[facePtBI];
321  }
322 
323  const tetPointRef tetB
324  (
325  ccB,
326  tetBasePtB,
327  meshB.points()[pt0I],
328  meshB.points()[pt1I]
329  );
330 
331  if (!tetB.bounds().overlaps(tetABb))
332  {
333  continue;
334  }
335 
336  vol += tetTetOverlapVol(tetA, tetB);
337 
338  if (vol > threshold)
339  {
340  return true;
341  }
342  }
343  }
344  }
345  }
346 
347  return false;
348 }
349 
350 
352 (
353  const primitiveMesh& meshA,
354  const label cellAI,
355  const primitiveMesh& meshB,
356  const label cellBI,
357  const treeBoundBox& cellBbB
358 ) const
359 {
360  const cell& cFacesA = meshA.cells()[cellAI];
361  const point& ccA = meshA.cellCentres()[cellAI];
362 
363  const cell& cFacesB = meshB.cells()[cellBI];
364  const point& ccB = meshB.cellCentres()[cellBI];
365 
366  scalar vol = 0.0;
367 
368  forAll(cFacesA, cFA)
369  {
370  label faceAI = cFacesA[cFA];
371 
372  const face& fA = meshA.faces()[faceAI];
373  const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
374  if (!pyrA.overlaps(cellBbB))
375  {
376  continue;
377  }
378 
379  bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
380 
381  label tetBasePtAI = 0;
382 
383  const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
384 
385  for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
386  {
387  label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
388  label otherFacePtAI = fA.fcIndex(facePtAI);
389 
390  label pt0I = -1;
391  label pt1I = -1;
392 
393  if (ownA)
394  {
395  pt0I = fA[facePtAI];
396  pt1I = fA[otherFacePtAI];
397  }
398  else
399  {
400  pt0I = fA[otherFacePtAI];
401  pt1I = fA[facePtAI];
402  }
403 
404  const tetPointRef tetA
405  (
406  ccA,
407  tetBasePtA,
408  meshA.points()[pt0I],
409  meshA.points()[pt1I]
410  );
411  const treeBoundBox tetABb(tetA.bounds());
412 
413 
414  // Loop over tets of cellB
415  forAll(cFacesB, cFB)
416  {
417  label faceBI = cFacesB[cFB];
418 
419  const face& fB = meshB.faces()[faceBI];
420  const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
421  if (!pyrB.overlaps(pyrA))
422  {
423  continue;
424  }
425 
426  bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
427 
428  label tetBasePtBI = 0;
429 
430  const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
431 
432  for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
433  {
434  label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
435  label otherFacePtBI = fB.fcIndex(facePtBI);
436 
437  label pt0I = -1;
438  label pt1I = -1;
439 
440  if (ownB)
441  {
442  pt0I = fB[facePtBI];
443  pt1I = fB[otherFacePtBI];
444  }
445  else
446  {
447  pt0I = fB[otherFacePtBI];
448  pt1I = fB[facePtBI];
449  }
450 
451  const tetPointRef tetB
452  (
453  ccB,
454  tetBasePtB,
455  meshB.points()[pt0I],
456  meshB.points()[pt1I]
457  );
458  if (!tetB.bounds().overlaps(tetABb))
459  {
460  continue;
461  }
462 
463  vol += tetTetOverlapVol(tetA, tetB);
464  }
465  }
466  }
467  }
468 
469  return vol;
470 }
471 
472 
473 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
bool overlaps(const boundBox &) const
Overlaps/touches boundingBox?
Definition: boundBoxI.H:126
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:60
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:76
Non-pointer based hierarchical recursive searching.
Definition: indexedOctree.H:72
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:61
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1313
const indexedOctree< treeDataCell > & cellTree() const
Return the cell search tree.
Definition: polyMesh.C:1057
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:75
virtual const faceList & faces() const =0
Return faces.
const vectorField & cellCentres() const
virtual const labelList & faceOwner() const =0
Face face-owner addressing.
const labelListList & cellPoints() const
virtual const pointField & points() const =0
Return mesh points.
const cellList & cells() const
Calculates the overlap volume of two cells using tetrahedral decomposition.
scalar cellCellOverlapVolumeMinDecomp(const primitiveMesh &meshA, const label cellAI, const primitiveMesh &meshB, const label cellBI, const treeBoundBox &cellBbB) const
Calculates the overlap volume.
scalar cellVolumeMinDecomp(const primitiveMesh &mesh, const label celli) const
Calculates the cell volume.
labelList overlappingCells(const polyMesh &meshA, const polyMesh &meshB, const label cellBI) const
Return a list of cells in meshA which overlaps with cellBI in.
bool cellCellOverlapMinDecomp(const primitiveMesh &meshA, const label cellAI, const primitiveMesh &meshB, const label cellBI, const treeBoundBox &cellBbB, const scalar threshold=0.0) const
Return true if olverlap volume is greater than threshold.
tetOverlapVolume()
Null constructor.
A tetrahedron primitive.
Definition: tetrahedron.H:82
boundBox bounds() const
Calculate the bounding box.
Definition: tetrahedronI.H:463
const Point & a() const
Return vertices.
Definition: tetrahedronI.H:72
const Point & c() const
Definition: tetrahedronI.H:86
scalar mag() const
Return volume.
Definition: tetrahedronI.H:170
const Point & b() const
Definition: tetrahedronI.H:79
const Point & d() const
Definition: tetrahedronI.H:93
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:90
bool overlaps(const boundBox &) const
Overlaps other bounding box?
Definition: boundBoxI.H:126
const pointField & points
Namespace for OpenFOAM.
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
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
vector point
Point is a vector.
Definition: point.H:41
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
defineTypeNameAndDebug(combustionModel, 0)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
cutTriTet::opAddResult< AboveOp, BelowOp >::type tetCut(const FixedList< Point, 4 > &tet, const FixedList< scalar, 4 > &level, const AboveOp &aboveOp, const BelowOp &belowOp)
As triCut, but for a tetrahedron.
labelList f(nPoints)