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 -------------------------------------------------------------------------------
11 License
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 
36 template<class Type>
37 Type 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 
78 template<class Type>
79 void 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 
267 template<class Type>
268 void 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 
470 template<class Type>
471 Foam::tmp<Foam::Field<Type>>
472 Foam::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::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
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))
Definition: gmvOutputSpray.H:25
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::DynamicList< Type >
isoSurfaceCell.H
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
polyMesh.H
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::triPoints
Triangle storage. Null constructable (unfortunately triangle<point, point> is not)
Definition: triPoints.H:52
Foam::Field
Generic templated field type.
Definition: Field.H:63
isoSurfacePoint.H
Foam::isoSurfaceBase::iso_
const scalar iso_
Iso value.
Definition: isoSurfaceBase.H:111
f
labelList f(nPoints)
Foam::List< label >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
tetMatcher.H
Foam::tetMatcher::test
static bool test(const UList< face > &faces)
Definition: tetMatcher.C:89
p0
const volScalarField & p0
Definition: EEqn.H:36
triFace
face triFace(3)
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328