isoSurfaceCellTemplates.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) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2020-2021 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 "isoSurfaceCell.H"
30#include "isoSurfacePoint.H"
31#include "polyMesh.H"
32#include "tetMatcher.H"
33
34// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
35
36template<class Type>
37Type Foam::isoSurfaceCell::generatePoint
38(
39 const DynamicList<Type>& snappedPoints,
40
41 const scalar s0,
42 const Type& p0,
43 const label p0Index,
44
45 const scalar s1,
46 const Type& p1,
47 const label p1Index
48) const
49{
50 const scalar d = s1-s0;
51
52 if (mag(d) > VSMALL)
53 {
54 const scalar s = (iso_-s0)/d;
55
56 if (s >= 0.5 && s <= 1 && p1Index != -1)
57 {
58 return snappedPoints[p1Index];
59 }
60 else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
61 {
62 return snappedPoints[p0Index];
63 }
64 else
65 {
66 return s*p1 + (1.0-s)*p0;
67 }
68 }
69 else
70 {
71 constexpr scalar s = 0.4999;
72
73 return s*p1 + (1.0-s)*p0;
74 }
75}
76
77
78template<class Type>
79void Foam::isoSurfaceCell::generateTriPoints
80(
81 const DynamicList<Type>& snapped,
82
83 const scalar s0,
84 const Type& p0,
85 const label p0Index,
86
87 const scalar s1,
88 const Type& p1,
89 const label p1Index,
90
91 const scalar s2,
92 const Type& p2,
93 const label p2Index,
94
95 const scalar s3,
96 const Type& p3,
97 const label p3Index,
98
99 DynamicList<Type>& pts
100) const
101{
102 int triIndex = 0;
103 if (s0 < iso_)
104 {
105 triIndex |= 1;
106 }
107 if (s1 < iso_)
108 {
109 triIndex |= 2;
110 }
111 if (s2 < iso_)
112 {
113 triIndex |= 4;
114 }
115 if (s3 < iso_)
116 {
117 triIndex |= 8;
118 }
119
120 // Form the vertices of the triangles for each case
121 switch (triIndex)
122 {
123 case 0x00:
124 case 0x0F:
125 break;
126
127 case 0x01:
128 case 0x0E:
129 {
130 pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
131 pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
132 pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
133 if (triIndex == 0x0E)
134 {
135 // Flip normals
136 const label sz = pts.size();
137 std::swap(pts[sz-2], pts[sz-1]);
138 }
139 }
140 break;
141
142 case 0x02:
143 case 0x0D:
144 {
145 pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
146 pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
147 pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
148
149 if (triIndex == 0x0D)
150 {
151 // Flip normals
152 const label sz = pts.size();
153 std::swap(pts[sz-2], pts[sz-1]);
154 }
155 }
156 break;
157
158 case 0x03:
159 case 0x0C:
160 {
161 Type p0p2 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
162 Type p1p3 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
163
164 pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
165 pts.append(p1p3);
166 pts.append(p0p2);
167
168 pts.append(p1p3);
169 pts.append
170 (
171 generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
172 );
173 pts.append(p0p2);
174
175 if (triIndex == 0x0C)
176 {
177 // Flip normals
178 const label sz = pts.size();
179 std::swap(pts[sz-5], pts[sz-4]);
180 std::swap(pts[sz-2], pts[sz-1]);
181 }
182 }
183 break;
184
185 case 0x04:
186 case 0x0B:
187 {
188 pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
189 pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
190 pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
191
192 if (triIndex == 0x0B)
193 {
194 // Flip normals
195 const label sz = pts.size();
196 std::swap(pts[sz-2], pts[sz-1]);
197 }
198 }
199 break;
200
201 case 0x05:
202 case 0x0A:
203 {
204 Type p0p1 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
205 Type p2p3 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
206
207 pts.append(p0p1);
208 pts.append(p2p3);
209 pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
210
211 pts.append(p0p1);
212 pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
213 pts.append(p2p3);
214
215 if (triIndex == 0x0A)
216 {
217 // Flip normals
218 const label sz = pts.size();
219 std::swap(pts[sz-5], pts[sz-4]);
220 std::swap(pts[sz-2], pts[sz-1]);
221 }
222 }
223 break;
224
225 case 0x06:
226 case 0x09:
227 {
228 Type p0p1 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
229 Type p2p3 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
230
231 pts.append(p0p1);
232 pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
233 pts.append(p2p3);
234
235 pts.append(p0p1);
236 pts.append(p2p3);
237 pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
238
239 if (triIndex == 0x09)
240 {
241 // Flip normals
242 const label sz = pts.size();
243 std::swap(pts[sz-5], pts[sz-4]);
244 std::swap(pts[sz-2], pts[sz-1]);
245 }
246 }
247 break;
248
249 case 0x08:
250 case 0x07:
251 {
252 pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
253 pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
254 pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
255 if (triIndex == 0x07)
256 {
257 // Flip normals
258 const label sz = pts.size();
259 std::swap(pts[sz-2], pts[sz-1]);
260 }
261 }
262 break;
263 }
264}
265
266
267template<class Type>
268void Foam::isoSurfaceCell::generateTriPoints
269(
270 const scalarField& cVals,
271 const scalarField& pVals,
272
273 const Field<Type>& cCoords,
274 const Field<Type>& pCoords,
275
276 const DynamicList<Type>& snappedPoints,
277 const labelList& snappedCc,
278 const labelList& snappedPoint,
279
280 DynamicList<Type>& triPoints,
281 DynamicList<label>& triMeshCells
282) const
283{
284 label countNotFoundTets = 0;
285
286 forAll(mesh_.cells(), celli)
287 {
288 if ((cellCutType_[celli] & cutType::ANYCUT) != 0)
289 {
290 label oldNPoints = triPoints.size();
291
292 const cell& cFaces = mesh_.cells()[celli];
293
294 if (tetMatcher::test(mesh_, celli))
295 {
296 // For tets don't do cell-centre decomposition, just use the
297 // tet points and values
298
299 const face& f0 = mesh_.faces()[cFaces[0]];
300
301 // Get the other point
302 const face& f1 = mesh_.faces()[cFaces[1]];
303 label oppositeI = -1;
304 forAll(f1, fp)
305 {
306 oppositeI = f1[fp];
307
308 if (!f0.found(oppositeI))
309 {
310 break;
311 }
312 }
313
314 // Start off from positive volume tet to make sure we
315 // generate outwards pointing tets
316 if (mesh_.faceOwner()[cFaces[0]] == celli)
317 {
318 generateTriPoints
319 (
320 snappedPoints,
321
322 pVals[f0[1]],
323 pCoords[f0[1]],
324 snappedPoint[f0[1]],
325
326 pVals[f0[0]],
327 pCoords[f0[0]],
328 snappedPoint[f0[0]],
329
330 pVals[f0[2]],
331 pCoords[f0[2]],
332 snappedPoint[f0[2]],
333
334 pVals[oppositeI],
335 pCoords[oppositeI],
336 snappedPoint[oppositeI],
337
338 triPoints
339 );
340 }
341 else
342 {
343 generateTriPoints
344 (
345 snappedPoints,
346
347 pVals[f0[0]],
348 pCoords[f0[0]],
349 snappedPoint[f0[0]],
350
351 pVals[f0[1]],
352 pCoords[f0[1]],
353 snappedPoint[f0[1]],
354
355 pVals[f0[2]],
356 pCoords[f0[2]],
357 snappedPoint[f0[2]],
358
359 pVals[oppositeI],
360 pCoords[oppositeI],
361 snappedPoint[oppositeI],
362
363 triPoints
364 );
365 }
366 }
367 else
368 {
369 forAll(cFaces, cFacei)
370 {
371 label facei = cFaces[cFacei];
372 const face& f = mesh_.faces()[facei];
373
374 label fp0 = mesh_.tetBasePtIs()[facei];
375
376 // Skip undefined tets
377 if (fp0 < 0)
378 {
379 fp0 = 0;
380 countNotFoundTets++;
381 }
382
383 label fp = f.fcIndex(fp0);
384 for (label i = 2; i < f.size(); i++)
385 {
386 label nextFp = f.fcIndex(fp);
387 triFace tri(f[fp0], f[fp], f[nextFp]);
388
389 // Start off from positive volume tet to make sure we
390 // generate outwards pointing tets
391 if (mesh_.faceOwner()[facei] == celli)
392 {
393 generateTriPoints
394 (
395 snappedPoints,
396
397 pVals[tri[1]],
398 pCoords[tri[1]],
399 snappedPoint[tri[1]],
400
401 pVals[tri[0]],
402 pCoords[tri[0]],
403 snappedPoint[tri[0]],
404
405 pVals[tri[2]],
406 pCoords[tri[2]],
407 snappedPoint[tri[2]],
408
409 cVals[celli],
410 cCoords[celli],
411 snappedCc[celli],
412
413 triPoints
414 );
415 }
416 else
417 {
418 generateTriPoints
419 (
420 snappedPoints,
421
422 pVals[tri[0]],
423 pCoords[tri[0]],
424 snappedPoint[tri[0]],
425
426 pVals[tri[1]],
427 pCoords[tri[1]],
428 snappedPoint[tri[1]],
429
430 pVals[tri[2]],
431 pCoords[tri[2]],
432 snappedPoint[tri[2]],
433
434 cVals[celli],
435 cCoords[celli],
436 snappedCc[celli],
437
438 triPoints
439 );
440 }
441
442 fp = nextFp;
443 }
444 }
445 }
446
447
448 // Every three triPoints is a cell
449 label nCells = (triPoints.size()-oldNPoints)/3;
450 for (label i = 0; i < nCells; i++)
451 {
452 triMeshCells.append(celli);
453 }
454 }
455 }
456
457 if (countNotFoundTets > 0)
458 {
459 WarningInFunction
460 << "Could not find " << countNotFoundTets
461 << " tet base points, which may lead to inverted triangles."
462 << endl;
463 }
464
465 triPoints.shrink();
466 triMeshCells.shrink();
467}
468
469
470template<class Type>
471Foam::tmp<Foam::Field<Type>>
472Foam::isoSurfaceCell::interpolateTemplate
473(
474 const Field<Type>& cCoords,
475 const Field<Type>& pCoords
476) const
477{
478 DynamicList<Type> triPoints(3*nCutCells_);
479 DynamicList<label> triMeshCells(nCutCells_);
480
481 // Dummy snap data
482 DynamicList<Type> snappedPoints;
483 labelList snappedCc(mesh_.nCells(), -1);
484 labelList snappedPoint(mesh_.nPoints(), -1);
485
486 generateTriPoints
487 (
488 cVals_,
489 pVals_,
490
491 cCoords,
492 pCoords,
493
494 snappedPoints,
495 snappedCc,
496 snappedPoint,
497
498 triPoints,
499 triMeshCells
500 );
501
502 return isoSurfacePoint::interpolate
503 (
504 this->points().size(),
505 triPointMergeMap_,
506 interpolatedPoints_,
507 interpolatedOldPoints_,
508 interpolationWeights_,
509 triPoints
510 );
511}
512
513
514// ************************************************************************* //
Foam::DynamicList
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
Foam::Field
Generic templated field type.
Definition: Field.H:82
Foam::UList::size
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Foam::UList::fcIndex
label fcIndex(const label i) const noexcept
Definition: UListI.H:60
Foam::isoSurfaceBase::iso_
const scalar iso_
Iso value.
Definition: isoSurfaceBase.H:111
Foam::sampledSurface::interpolate
bool interpolate() const noexcept
Same as isPointData()
Definition: sampledSurface.H:549
Foam::tetMatcher::test
static bool test(const UList< face > &faces)
Definition: tetMatcher.C:89
Foam::tmp
A class for managing temporary objects.
Definition: tmp.H:65
Foam::triPoints
Triangle storage. Null constructable (unfortunately triangle<point, point> is not)
Definition: triPoints.H:55
p0
const volScalarField & p0
Definition: EEqn.H:36
points
const pointField & points
Definition: gmvOutputHeader.H:1
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:66
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
triFace
face triFace(3)
f
labelList f(nPoints)
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333