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-2018 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 "cut.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(tetOverlapVolume, 0);
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  // A maximum of three cuts are made (the tets that result from the final cut
57  // are not stored), and each cut can create at most three tets. The
58  // temporary storage must therefore extend to 3^3 = 27 tets.
59  typedef cutTetList<27> tetListType;
60  static tetListType cutTetList1, cutTetList2;
61 
62  // face 0
63  const plane pl0(tetB.b(), tetB.d(), tetB.c());
64  const FixedList<point, 4> t({tetA.a(), tetA.b(), tetA.c(), tetA.d()});
65  cutTetList1.clear();
66  tetCut(t, pl0, cut::appendOp<tetListType>(cutTetList1), cut::noOp());
67  if (cutTetList1.size() == 0)
68  {
69  return 0;
70  }
71 
72  // face 1
73  const plane pl1(tetB.a(), tetB.c(), tetB.d());
74  cutTetList2.clear();
75  for (label i = 0; i < cutTetList1.size(); i++)
76  {
77  const FixedList<point, 4>& t = cutTetList1[i];
78  tetCut(t, pl1, cut::appendOp<tetListType>(cutTetList2), cut::noOp());
79  }
80  if (cutTetList2.size() == 0)
81  {
82  return 0;
83  }
84 
85  // face 2
86  const plane pl2(tetB.a(), tetB.d(), tetB.b());
87  cutTetList1.clear();
88  for (label i = 0; i < cutTetList2.size(); i++)
89  {
90  const FixedList<point, 4>& t = cutTetList2[i];
91  tetCut(t, pl2, cut::appendOp<tetListType>(cutTetList1), cut::noOp());
92  }
93  if (cutTetList1.size() == 0)
94  {
95  return 0;
96  }
97 
98  // face 3
99  const plane pl3(tetB.a(), tetB.b(), tetB.c());
100  scalar v = 0;
101  for (label i = 0; i < cutTetList1.size(); i++)
102  {
103  const FixedList<point, 4>& t = cutTetList1[i];
104  v += tetCut(t, pl3, cut::volumeOp(), cut::noOp());
105  }
106 
107  return v;
108 }
109 
110 
111 Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
112 (
113  const pointField& points,
114  const face& f,
115  const point& fc
116 ) const
117 {
118  treeBoundBox bb(fc, fc);
119  forAll(f, fp)
120  {
121  const point& pt = points[f[fp]];
122  bb.min() = min(bb.min(), pt);
123  bb.max() = max(bb.max(), pt);
124  }
125  return bb;
126 }
127 
128 
129 // * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * * //
130 
132 (
133  const primitiveMesh& meshA,
134  const label cellAI,
135  const primitiveMesh& meshB,
136  const label cellBI,
137  const treeBoundBox& cellBbB,
138  const scalar threshold
139 ) const
140 {
141  const cell& cFacesA = meshA.cells()[cellAI];
142  const point& ccA = meshA.cellCentres()[cellAI];
143 
144  const cell& cFacesB = meshB.cells()[cellBI];
145  const point& ccB = meshB.cellCentres()[cellBI];
146 
147  scalar vol = 0.0;
148 
149  forAll(cFacesA, cFA)
150  {
151  label faceAI = cFacesA[cFA];
152 
153  const face& fA = meshA.faces()[faceAI];
154  const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
155  if (!pyrA.overlaps(cellBbB))
156  {
157  continue;
158  }
159 
160  bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
161 
162  label tetBasePtAI = 0;
163 
164  const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
165 
166  for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
167  {
168  label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
169  label otherFacePtAI = fA.fcIndex(facePtAI);
170 
171  label pt0I = -1;
172  label pt1I = -1;
173 
174  if (ownA)
175  {
176  pt0I = fA[facePtAI];
177  pt1I = fA[otherFacePtAI];
178  }
179  else
180  {
181  pt0I = fA[otherFacePtAI];
182  pt1I = fA[facePtAI];
183  }
184 
185  const tetPointRef tetA
186  (
187  ccA,
188  tetBasePtA,
189  meshA.points()[pt0I],
190  meshA.points()[pt1I]
191  );
192  const treeBoundBox tetABb(tetA.bounds());
193 
194 
195  // Loop over tets of cellB
196  forAll(cFacesB, cFB)
197  {
198  label faceBI = cFacesB[cFB];
199 
200  const face& fB = meshB.faces()[faceBI];
201  const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
202  if (!pyrB.overlaps(pyrA))
203  {
204  continue;
205  }
206 
207  bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
208 
209  label tetBasePtBI = 0;
210 
211  const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
212 
213  for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
214  {
215  label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
216  label otherFacePtBI = fB.fcIndex(facePtBI);
217 
218  label pt0I = -1;
219  label pt1I = -1;
220 
221  if (ownB)
222  {
223  pt0I = fB[facePtBI];
224  pt1I = fB[otherFacePtBI];
225  }
226  else
227  {
228  pt0I = fB[otherFacePtBI];
229  pt1I = fB[facePtBI];
230  }
231 
232  const tetPointRef tetB
233  (
234  ccB,
235  tetBasePtB,
236  meshB.points()[pt0I],
237  meshB.points()[pt1I]
238  );
239 
240  if (!tetB.bounds().overlaps(tetABb))
241  {
242  continue;
243  }
244 
245  vol += tetTetOverlapVol(tetA, tetB);
246 
247  if (vol > threshold)
248  {
249  return true;
250  }
251  }
252  }
253  }
254  }
255 
256  return false;
257 }
258 
259 
261 (
262  const primitiveMesh& meshA,
263  const label cellAI,
264  const primitiveMesh& meshB,
265  const label cellBI,
266  const treeBoundBox& cellBbB
267 ) const
268 {
269  const cell& cFacesA = meshA.cells()[cellAI];
270  const point& ccA = meshA.cellCentres()[cellAI];
271 
272  const cell& cFacesB = meshB.cells()[cellBI];
273  const point& ccB = meshB.cellCentres()[cellBI];
274 
275  scalar vol = 0.0;
276 
277  forAll(cFacesA, cFA)
278  {
279  label faceAI = cFacesA[cFA];
280 
281  const face& fA = meshA.faces()[faceAI];
282  const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
283  if (!pyrA.overlaps(cellBbB))
284  {
285  continue;
286  }
287 
288  bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
289 
290  label tetBasePtAI = 0;
291 
292  const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
293 
294  for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
295  {
296  label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
297  label otherFacePtAI = fA.fcIndex(facePtAI);
298 
299  label pt0I = -1;
300  label pt1I = -1;
301 
302  if (ownA)
303  {
304  pt0I = fA[facePtAI];
305  pt1I = fA[otherFacePtAI];
306  }
307  else
308  {
309  pt0I = fA[otherFacePtAI];
310  pt1I = fA[facePtAI];
311  }
312 
313  const tetPointRef tetA
314  (
315  ccA,
316  tetBasePtA,
317  meshA.points()[pt0I],
318  meshA.points()[pt1I]
319  );
320  const treeBoundBox tetABb(tetA.bounds());
321 
322 
323  // Loop over tets of cellB
324  forAll(cFacesB, cFB)
325  {
326  label faceBI = cFacesB[cFB];
327 
328  const face& fB = meshB.faces()[faceBI];
329  const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
330  if (!pyrB.overlaps(pyrA))
331  {
332  continue;
333  }
334 
335  bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
336 
337  label tetBasePtBI = 0;
338 
339  const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
340 
341  for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
342  {
343  label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
344  label otherFacePtBI = fB.fcIndex(facePtBI);
345 
346  label pt0I = -1;
347  label pt1I = -1;
348 
349  if (ownB)
350  {
351  pt0I = fB[facePtBI];
352  pt1I = fB[otherFacePtBI];
353  }
354  else
355  {
356  pt0I = fB[otherFacePtBI];
357  pt1I = fB[facePtBI];
358  }
359 
360  const tetPointRef tetB
361  (
362  ccB,
363  tetBasePtB,
364  meshB.points()[pt0I],
365  meshB.points()[pt1I]
366  );
367  if (!tetB.bounds().overlaps(tetABb))
368  {
369  continue;
370  }
371 
372  vol += tetTetOverlapVol(tetA, tetB);
373  }
374  }
375  }
376  }
377 
378  return vol;
379 }
380 
381 
383 (
384  const polyMesh& fromMesh,
385  const polyMesh& toMesh,
386  const label iTo
387 ) const
388 {
389  const indexedOctree<treeDataCell>& treeA = fromMesh.cellTree();
390 
391  treeBoundBox bbB(toMesh.points(), toMesh.cellPoints()[iTo]);
392 
393  return treeA.findBox(bbB);
394 }
395 
396 
397 // ************************************************************************* //
bool overlaps(const boundBox &) const
Overlaps/touches boundingBox?
Definition: boundBoxI.H:120
A tetrahedron primitive.
Definition: tetrahedron.H:61
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:54
const Point & c() const
Definition: tetrahedronI.H:86
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
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
const Point & a() const
Return vertices.
Definition: tetrahedronI.H:72
virtual const pointField & points() const =0
Return mesh points.
const cellList & cells() const
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1211
const Point & d() const
Definition: tetrahedronI.H:93
const Point & b() const
Definition: tetrahedronI.H:79
tetOverlapVolume()
Null constructor.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
const vectorField & cellCentres() const
defineTypeNameAndDebug(combustionModel, 0)
const point & max() const
Maximum point defining the bounding box.
Definition: boundBoxI.H:60
const indexedOctree< treeDataCell > & cellTree() const
Return the cell search tree.
Definition: polyMesh.C:962
boundBox bounds() const
Calculate the bounding box.
Definition: tetrahedron.C:340
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
virtual const faceList & faces() const =0
Return faces.
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
virtual const labelList & faceOwner() const =0
Face face-owner addressing.
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.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
const point & min() const
Minimum point defining the bounding box.
Definition: boundBoxI.H:54
labelList overlappingCells(const polyMesh &meshA, const polyMesh &meshB, const label cellBI) const
Return a list of cells in meshA which overlaps with cellBI in.
scalar cellCellOverlapVolumeMinDecomp(const primitiveMesh &meshA, const label cellAI, const primitiveMesh &meshB, const label cellBI, const treeBoundBox &cellBbB) const
Calculates the overlap volume.
const labelListList & cellPoints() const
Namespace for OpenFOAM.
cut::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.