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 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "isoSurfaceCell.H"
29 #include "isoSurface.H"
30 #include "polyMesh.H"
31 #include "tetMatcher.H"
32 
33 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
34 
35 template<class Type>
36 Type Foam::isoSurfaceCell::generatePoint
37 (
38  const DynamicList<Type>& snappedPoints,
39 
40  const scalar s0,
41  const Type& p0,
42  const label p0Index,
43 
44  const scalar s1,
45  const Type& p1,
46  const label p1Index
47 ) const
48 {
49  const scalar d = s1-s0;
50 
51  if (mag(d) > VSMALL)
52  {
53  const scalar s = (iso_-s0)/d;
54 
55  if (s >= 0.5 && s <= 1 && p1Index != -1)
56  {
57  return snappedPoints[p1Index];
58  }
59  else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
60  {
61  return snappedPoints[p0Index];
62  }
63  else
64  {
65  return s*p1 + (1.0-s)*p0;
66  }
67  }
68  else
69  {
70  scalar s = 0.4999;
71 
72  return s*p1 + (1.0-s)*p0;
73  }
74 }
75 
76 
77 template<class Type>
78 void Foam::isoSurfaceCell::generateTriPoints
79 (
80  const DynamicList<Type>& snapped,
81 
82  const scalar s0,
83  const Type& p0,
84  const label p0Index,
85 
86  const scalar s1,
87  const Type& p1,
88  const label p1Index,
89 
90  const scalar s2,
91  const Type& p2,
92  const label p2Index,
93 
94  const scalar s3,
95  const Type& p3,
96  const label p3Index,
97 
98  DynamicList<Type>& pts
99 ) const
100 {
101  int triIndex = 0;
102  if (s0 < iso_)
103  {
104  triIndex |= 1;
105  }
106  if (s1 < iso_)
107  {
108  triIndex |= 2;
109  }
110  if (s2 < iso_)
111  {
112  triIndex |= 4;
113  }
114  if (s3 < iso_)
115  {
116  triIndex |= 8;
117  }
118 
119  // Form the vertices of the triangles for each case
120  switch (triIndex)
121  {
122  case 0x00:
123  case 0x0F:
124  break;
125 
126  case 0x01:
127  case 0x0E:
128  {
129  pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
130  pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
131  pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
132  if (triIndex == 0x0E)
133  {
134  // Flip normals
135  const label sz = pts.size();
136  Swap(pts[sz-2], pts[sz-1]);
137  }
138  }
139  break;
140 
141  case 0x02:
142  case 0x0D:
143  {
144  pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
145  pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
146  pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
147 
148  if (triIndex == 0x0D)
149  {
150  // Flip normals
151  const label sz = pts.size();
152  Swap(pts[sz-2], pts[sz-1]);
153  }
154  }
155  break;
156 
157  case 0x03:
158  case 0x0C:
159  {
160  Type p0p2 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
161  Type p1p3 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
162 
163  pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
164  pts.append(p1p3);
165  pts.append(p0p2);
166 
167  pts.append(p1p3);
168  pts.append
169  (
170  generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
171  );
172  pts.append(p0p2);
173 
174  if (triIndex == 0x0C)
175  {
176  // Flip normals
177  label sz = pts.size();
178  Swap(pts[sz-5], pts[sz-4]);
179  Swap(pts[sz-2], pts[sz-1]);
180  }
181  }
182  break;
183 
184  case 0x04:
185  case 0x0B:
186  {
187  pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
188  pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
189  pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
190 
191  if (triIndex == 0x0B)
192  {
193  // Flip normals
194  const label sz = pts.size();
195  Swap(pts[sz-2], pts[sz-1]);
196  }
197  }
198  break;
199 
200  case 0x05:
201  case 0x0A:
202  {
203  Type p0p1 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
204  Type p2p3 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
205 
206  pts.append(p0p1);
207  pts.append(p2p3);
208  pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
209 
210  pts.append(p0p1);
211  pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
212  pts.append(p2p3);
213 
214  if (triIndex == 0x0A)
215  {
216  // Flip normals
217  const label sz = pts.size();
218  Swap(pts[sz-5], pts[sz-4]);
219  Swap(pts[sz-2], pts[sz-1]);
220  }
221  }
222  break;
223 
224  case 0x06:
225  case 0x09:
226  {
227  Type p0p1 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
228  Type p2p3 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
229 
230  pts.append(p0p1);
231  pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
232  pts.append(p2p3);
233 
234  pts.append(p0p1);
235  pts.append(p2p3);
236  pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
237 
238  if (triIndex == 0x09)
239  {
240  // Flip normals
241  const label sz = pts.size();
242  Swap(pts[sz-5], pts[sz-4]);
243  Swap(pts[sz-2], pts[sz-1]);
244  }
245  }
246  break;
247 
248  case 0x08:
249  case 0x07:
250  {
251  pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
252  pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
253  pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
254  if (triIndex == 0x07)
255  {
256  // Flip normals
257  const label sz = pts.size();
258  Swap(pts[sz-2], pts[sz-1]);
259  }
260  }
261  break;
262  }
263 }
264 
265 
266 template<class Type>
267 void Foam::isoSurfaceCell::generateTriPoints
268 (
269  const scalarField& cVals,
270  const scalarField& pVals,
271 
272  const Field<Type>& cCoords,
273  const Field<Type>& pCoords,
274 
275  const DynamicList<Type>& snappedPoints,
276  const labelList& snappedCc,
277  const labelList& snappedPoint,
278 
279  DynamicList<Type>& triPoints,
280  DynamicList<label>& triMeshCells
281 ) const
282 {
283  tetMatcher tet;
284  label countNotFoundTets = 0;
285 
286  forAll(mesh_.cells(), celli)
287  {
288  if (cellCutType_[celli] != NOTCUT)
289  {
290  label oldNPoints = triPoints.size();
291 
292  const cell& cFaces = mesh_.cells()[celli];
293 
294  if (tet.isA(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::interpolate
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 isoSurface::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:71
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:59
Foam::DynamicList
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:55
isoSurfaceCell.H
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:350
polyMesh.H
Foam::Swap
void Swap(DynamicList< T, SizeMin1 > &a, DynamicList< T, SizeMin2 > &b)
Definition: DynamicListI.H:909
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
Foam::isoSurfaceBase::iso_
const scalar iso_
Iso value.
Definition: isoSurfaceBase.H:66
isoSurface.H
f
labelList f(nPoints)
Foam::List< label >
Foam::isoSurfaceCell::interpolate
tmp< Field< Type > > interpolate(const Field< Type > &cCoords, const Field< Type > &pCoords) const
Interpolates cCoords, pCoords.
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
tetMatcher.H
p0
const volScalarField & p0
Definition: EEqn.H:36
triFace
face triFace(3)
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:298