tetOverlapVolumeTemplates.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2015 OpenFOAM Foundation
9 Copyright (C) 2016 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "tetOverlapVolume.H"
30#include "primitiveMesh.H"
31
32// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
33
34template<class tetPointsOp>
35void Foam::tetOverlapVolume::tetTetOverlap
36(
37 const tetPoints& tetA,
38 const tetPoints& tetB,
39 tetPointsOp& insideOp
40)
41{
42 static tetPointRef::tetIntersectionList insideTets;
43 label nInside = 0;
44 static tetPointRef::tetIntersectionList cutInsideTets;
45 label nCutInside = 0;
46
47 tetPointRef::storeOp inside(insideTets, nInside);
48 tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
49 tetPointRef::dummyOp outside;
50
51 // Precompute the tet face areas and exit early if any face area is
52 // too small
53 static FixedList<vector, 4> tetAFaceAreas;
54 static FixedList<scalar, 4> tetAMag2FaceAreas;
55 tetPointRef tetATet = tetA.tet();
56 for (label facei = 0; facei < 4; ++facei)
57 {
58 tetAFaceAreas[facei] = -tetATet.tri(facei).areaNormal();
59 tetAMag2FaceAreas[facei] = magSqr(tetAFaceAreas[facei]);
60 if (tetAMag2FaceAreas[facei] < ROOTVSMALL)
61 {
62 return;
63 }
64 }
65
66 static FixedList<vector, 4> tetBFaceAreas;
67 static FixedList<scalar, 4> tetBMag2FaceAreas;
68 tetPointRef tetBTet = tetB.tet();
69 for (label facei = 0; facei < 4; ++facei)
70 {
71 tetBFaceAreas[facei] = -tetBTet.tri(facei).areaNormal();
72 tetBMag2FaceAreas[facei] = magSqr(tetBFaceAreas[facei]);
73 if (tetBMag2FaceAreas[facei] < ROOTVSMALL)
74 {
75 return;
76 }
77 }
78
79
80 // Face 0
81 {
82 vector n = tetBFaceAreas[0]/Foam::sqrt(tetBMag2FaceAreas[0]);
83 plane pl0(tetBTet.tri(0).a(), n, false);
84
85 tetA.tet().sliceWithPlane(pl0, cutInside, outside);
86 if (nCutInside == 0)
87 {
88 return;
89 }
90 }
91
92 // Face 1
93 {
94 vector n = tetBFaceAreas[1]/Foam::sqrt(tetBMag2FaceAreas[1]);
95 plane pl1(tetBTet.tri(1).a(), n, false);
96
97 nInside = 0;
98 for (label i = 0; i < nCutInside; i++)
99 {
100 const tetPointRef t = cutInsideTets[i].tet();
101 t.sliceWithPlane(pl1, inside, outside);
102 }
103 if (nInside == 0)
104 {
105 return;
106 }
107 }
108
109 // Face 2
110 {
111 vector n = tetBFaceAreas[2]/Foam::sqrt(tetBMag2FaceAreas[2]);
112 plane pl2(tetBTet.tri(2).a(), n, false);
113
114 nCutInside = 0;
115 for (label i = 0; i < nInside; i++)
116 {
117 const tetPointRef t = insideTets[i].tet();
118 t.sliceWithPlane(pl2, cutInside, outside);
119 }
120 if (nCutInside == 0)
121 {
122 return;
123 }
124 }
125
126 // Face 3
127 {
128 vector n = tetBFaceAreas[3]/Foam::sqrt(tetBMag2FaceAreas[3]);
129 plane pl3(tetBTet.tri(3).a(), n, false);
130 for (label i = 0; i < nCutInside; i++)
131 {
132 const tetPointRef t = cutInsideTets[i].tet();
133 t.sliceWithPlane(pl3, insideOp, outside);
134 }
135 }
136}
137
138
139template<class tetsOp>
140void Foam::tetOverlapVolume::cellCellOverlapMinDecomp
141(
142 const primitiveMesh& meshA,
143 const label cellAI,
144
145 const primitiveMesh& meshB,
146 const label cellBI,
147 const treeBoundBox& cellBbB,
148 tetsOp& combineTetsOp
149)
150{
151 const cell& cFacesA = meshA.cells()[cellAI];
152 const point& ccA = meshA.cellCentres()[cellAI];
153
154 const cell& cFacesB = meshB.cells()[cellBI];
155 const point& ccB = meshB.cellCentres()[cellBI];
156
157 forAll(cFacesA, cFA)
158 {
159 label faceAI = cFacesA[cFA];
160
161 const face& fA = meshA.faces()[faceAI];
162 const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
163 if (!pyrA.overlaps(cellBbB))
164 {
165 continue;
166 }
167
168 bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
169
170 label tetBasePtAI = 0;
171
172 const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
173
174 for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
175 {
176 label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
177 label otherFacePtAI = fA.fcIndex(facePtAI);
178
179 label pt0I = -1;
180 label pt1I = -1;
181
182 if (ownA)
183 {
184 pt0I = fA[facePtAI];
185 pt1I = fA[otherFacePtAI];
186 }
187 else
188 {
189 pt0I = fA[otherFacePtAI];
190 pt1I = fA[facePtAI];
191 }
192
193 const tetPoints tetA
194 (
195 ccA,
196 tetBasePtA,
197 meshA.points()[pt0I],
198 meshA.points()[pt1I]
199 );
200 const treeBoundBox tetABb(tetA.bounds());
201
202 // Loop over tets of cellB
203 forAll(cFacesB, cFB)
204 {
205 label faceBI = cFacesB[cFB];
206
207 const face& fB = meshB.faces()[faceBI];
208 const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
209 if (!pyrB.overlaps(pyrA))
210 {
211 continue;
212 }
213
214 bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
215
216 label tetBasePtBI = 0;
217
218 const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
219
220 for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
221 {
222 label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
223 label otherFacePtBI = fB.fcIndex(facePtBI);
224
225 label pt0I = -1;
226 label pt1I = -1;
227
228 if (ownB)
229 {
230 pt0I = fB[facePtBI];
231 pt1I = fB[otherFacePtBI];
232 }
233 else
234 {
235 pt0I = fB[otherFacePtBI];
236 pt1I = fB[facePtBI];
237 }
238
239 const tetPoints tetB
240 (
241 ccB,
242 tetBasePtB,
243 meshB.points()[pt0I],
244 meshB.points()[pt1I]
245 );
246 if (!tetB.bounds().overlaps(tetABb))
247 {
248 continue;
249 }
250
251 if (combineTetsOp(tetA, tetB))
252 {
253 return;
254 }
255 }
256 }
257 }
258 }
259}
260
261
262// ************************************************************************* //
label n
static constexpr direction size() noexcept
The number of elements in the VectorSpace = Ncmpts.
Definition: VectorSpace.H:176
FixedList< tetPoints, 200 > tetIntersectionList
Storage type for tets originating from intersecting tets.
Definition: tetrahedron.H:94
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
tetrahedron< point, const point & > tetPointRef
A tetrahedron using referred points.
Definition: tetPointRef.H:47
vector point
Point is a vector.
Definition: point.H:43
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333