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 | Copyright (C) 2011-2015 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 "isoSurfaceCell.H"
27 #include "polyMesh.H"
28 #include "tetMatcher.H"
29 #include "isoSurface.H"
30 
31 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32 
33 template<class Type>
34 Type Foam::isoSurfaceCell::generatePoint
35 (
36  const DynamicList<Type>& snappedPoints,
37 
38  const scalar s0,
39  const Type& p0,
40  const label p0Index,
41 
42  const scalar s1,
43  const Type& p1,
44  const label p1Index
45 ) const
46 {
47  scalar d = s1-s0;
48 
49  if (mag(d) > VSMALL)
50  {
51  scalar s = (iso_-s0)/d;
52 
53  if (s >= 0.5 && s <= 1 && p1Index != -1)
54  {
55  return snappedPoints[p1Index];
56  }
57  else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
58  {
59  return snappedPoints[p0Index];
60  }
61  else
62  {
63  return s*p1 + (1.0-s)*p0;
64  }
65  }
66  else
67  {
68  scalar s = 0.4999;
69 
70  return s*p1 + (1.0-s)*p0;
71  }
72 }
73 
74 
75 template<class Type>
76 void Foam::isoSurfaceCell::generateTriPoints
77 (
78  const DynamicList<Type>& snapped,
79 
80  const scalar s0,
81  const Type& p0,
82  const label p0Index,
83 
84  const scalar s1,
85  const Type& p1,
86  const label p1Index,
87 
88  const scalar s2,
89  const Type& p2,
90  const label p2Index,
91 
92  const scalar s3,
93  const Type& p3,
94  const label p3Index,
95 
96  DynamicList<Type>& pts
97 ) const
98 {
99  int triIndex = 0;
100  if (s0 < iso_)
101  {
102  triIndex |= 1;
103  }
104  if (s1 < iso_)
105  {
106  triIndex |= 2;
107  }
108  if (s2 < iso_)
109  {
110  triIndex |= 4;
111  }
112  if (s3 < iso_)
113  {
114  triIndex |= 8;
115  }
116 
117  // Form the vertices of the triangles for each case
118  switch (triIndex)
119  {
120  case 0x00:
121  case 0x0F:
122  break;
123 
124  case 0x01:
125  case 0x0E:
126  {
127  pts.append
128  (
129  generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index)
130  );
131  pts.append
132  (
133  generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
134  );
135  pts.append
136  (
137  generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
138  );
139  if (triIndex == 0x0E)
140  {
141  // Flip normals
142  label sz = pts.size();
143  Swap(pts[sz-2], pts[sz-1]);
144  }
145  }
146  break;
147 
148  case 0x02:
149  case 0x0D:
150  {
151  pts.append
152  (
153  generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index)
154  );
155  pts.append
156  (
157  generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
158  );
159  pts.append
160  (
161  generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
162  );
163 
164  if (triIndex == 0x0D)
165  {
166  // Flip normals
167  label sz = pts.size();
168  Swap(pts[sz-2], pts[sz-1]);
169  }
170  }
171  break;
172 
173  case 0x03:
174  case 0x0C:
175  {
176  Type p0p2 =
177  generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
178  Type p1p3 =
179  generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
180 
181  pts.append
182  (
183  generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
184  );
185  pts.append(p1p3);
186  pts.append(p0p2);
187 
188  pts.append(p1p3);
189  pts.append
190  (
191  generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
192  );
193  pts.append(p0p2);
194 
195  if (triIndex == 0x0C)
196  {
197  // Flip normals
198  label sz = pts.size();
199  Swap(pts[sz-5], pts[sz-4]);
200  Swap(pts[sz-2], pts[sz-1]);
201  }
202  }
203  break;
204 
205  case 0x04:
206  case 0x0B:
207  {
208  pts.append
209  (
210  generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index)
211  );
212  pts.append
213  (
214  generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index)
215  );
216  pts.append
217  (
218  generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index)
219  );
220 
221  if (triIndex == 0x0B)
222  {
223  // Flip normals
224  label sz = pts.size();
225  Swap(pts[sz-2], pts[sz-1]);
226  }
227  }
228  break;
229 
230  case 0x05:
231  case 0x0A:
232  {
233  Type p0p1 =
234  generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
235  Type p2p3 =
236  generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
237 
238  pts.append(p0p1);
239  pts.append(p2p3);
240  pts.append
241  (
242  generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
243  );
244 
245  pts.append(p0p1);
246  pts.append
247  (
248  generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
249  );
250  pts.append(p2p3);
251 
252  if (triIndex == 0x0A)
253  {
254  // Flip normals
255  label sz = pts.size();
256  Swap(pts[sz-5], pts[sz-4]);
257  Swap(pts[sz-2], pts[sz-1]);
258  }
259  }
260  break;
261 
262  case 0x06:
263  case 0x09:
264  {
265  Type p0p1 =
266  generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
267  Type p2p3 =
268  generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
269 
270  pts.append(p0p1);
271  pts.append
272  (
273  generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
274  );
275  pts.append(p2p3);
276 
277  pts.append(p0p1);
278  pts.append(p2p3);
279  pts.append
280  (
281  generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
282  );
283 
284  if (triIndex == 0x09)
285  {
286  // Flip normals
287  label sz = pts.size();
288  Swap(pts[sz-5], pts[sz-4]);
289  Swap(pts[sz-2], pts[sz-1]);
290  }
291  }
292  break;
293 
294  case 0x08:
295  case 0x07:
296  {
297  pts.append
298  (
299  generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index)
300  );
301  pts.append
302  (
303  generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index)
304  );
305  pts.append
306  (
307  generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index)
308  );
309  if (triIndex == 0x07)
310  {
311  // Flip normals
312  label sz = pts.size();
313  Swap(pts[sz-2], pts[sz-1]);
314  }
315  }
316  break;
317  }
318 }
319 
320 
321 template<class Type>
322 void Foam::isoSurfaceCell::generateTriPoints
323 (
324  const scalarField& cVals,
325  const scalarField& pVals,
326 
327  const Field<Type>& cCoords,
328  const Field<Type>& pCoords,
329 
330  const DynamicList<Type>& snappedPoints,
331  const labelList& snappedCc,
332  const labelList& snappedPoint,
333 
334  DynamicList<Type>& triPoints,
335  DynamicList<label>& triMeshCells
336 ) const
337 {
338  tetMatcher tet;
339  label countNotFoundTets = 0;
340 
341  forAll(mesh_.cells(), cellI)
342  {
343  if (cellCutType_[cellI] != NOTCUT)
344  {
345  label oldNPoints = triPoints.size();
346 
347  const cell& cFaces = mesh_.cells()[cellI];
348 
349  if (tet.isA(mesh_, cellI))
350  {
351  // For tets don't do cell-centre decomposition, just use the
352  // tet points and values
353 
354  const face& f0 = mesh_.faces()[cFaces[0]];
355 
356  // Get the other point
357  const face& f1 = mesh_.faces()[cFaces[1]];
358  label oppositeI = -1;
359  forAll(f1, fp)
360  {
361  oppositeI = f1[fp];
362 
363  if (findIndex(f0, oppositeI) == -1)
364  {
365  break;
366  }
367  }
368 
369  // Start off from positive volume tet to make sure we
370  // generate outwards pointing tets
371  if (mesh_.faceOwner()[cFaces[0]] == cellI)
372  {
373  generateTriPoints
374  (
375  snappedPoints,
376 
377  pVals[f0[1]],
378  pCoords[f0[1]],
379  snappedPoint[f0[1]],
380 
381  pVals[f0[0]],
382  pCoords[f0[0]],
383  snappedPoint[f0[0]],
384 
385  pVals[f0[2]],
386  pCoords[f0[2]],
387  snappedPoint[f0[2]],
388 
389  pVals[oppositeI],
390  pCoords[oppositeI],
391  snappedPoint[oppositeI],
392 
393  triPoints
394  );
395  }
396  else
397  {
398  generateTriPoints
399  (
400  snappedPoints,
401 
402  pVals[f0[0]],
403  pCoords[f0[0]],
404  snappedPoint[f0[0]],
405 
406  pVals[f0[1]],
407  pCoords[f0[1]],
408  snappedPoint[f0[1]],
409 
410  pVals[f0[2]],
411  pCoords[f0[2]],
412  snappedPoint[f0[2]],
413 
414  pVals[oppositeI],
415  pCoords[oppositeI],
416  snappedPoint[oppositeI],
417 
418  triPoints
419  );
420  }
421  }
422  else
423  {
424  forAll(cFaces, cFaceI)
425  {
426  label faceI = cFaces[cFaceI];
427  const face& f = mesh_.faces()[faceI];
428 
429  label fp0 = mesh_.tetBasePtIs()[faceI];
430 
431  // Skip undefined tets
432  if (fp0 < 0)
433  {
434  fp0 = 0;
435  countNotFoundTets++;
436  }
437 
438  label fp = f.fcIndex(fp0);
439  for (label i = 2; i < f.size(); i++)
440  {
441  label nextFp = f.fcIndex(fp);
442  triFace tri(f[fp0], f[fp], f[nextFp]);
443 
444  // Start off from positive volume tet to make sure we
445  // generate outwards pointing tets
446  if (mesh_.faceOwner()[faceI] == cellI)
447  {
448  generateTriPoints
449  (
450  snappedPoints,
451 
452  pVals[tri[1]],
453  pCoords[tri[1]],
454  snappedPoint[tri[1]],
455 
456  pVals[tri[0]],
457  pCoords[tri[0]],
458  snappedPoint[tri[0]],
459 
460  pVals[tri[2]],
461  pCoords[tri[2]],
462  snappedPoint[tri[2]],
463 
464  cVals[cellI],
465  cCoords[cellI],
466  snappedCc[cellI],
467 
468  triPoints
469  );
470  }
471  else
472  {
473  generateTriPoints
474  (
475  snappedPoints,
476 
477  pVals[tri[0]],
478  pCoords[tri[0]],
479  snappedPoint[tri[0]],
480 
481  pVals[tri[1]],
482  pCoords[tri[1]],
483  snappedPoint[tri[1]],
484 
485  pVals[tri[2]],
486  pCoords[tri[2]],
487  snappedPoint[tri[2]],
488 
489  cVals[cellI],
490  cCoords[cellI],
491  snappedCc[cellI],
492 
493  triPoints
494  );
495  }
496 
497  fp = nextFp;
498  }
499  }
500  }
501 
502 
503  // Every three triPoints is a cell
504  label nCells = (triPoints.size()-oldNPoints)/3;
505  for (label i = 0; i < nCells; i++)
506  {
507  triMeshCells.append(cellI);
508  }
509  }
510  }
511 
512  if (countNotFoundTets > 0)
513  {
514  WarningInFunction
515  << "Could not find " << countNotFoundTets
516  << " tet base points, which may lead to inverted triangles."
517  << endl;
518  }
519 
520  triPoints.shrink();
521  triMeshCells.shrink();
522 }
523 
524 
525 template<class Type>
526 Foam::tmp<Foam::Field<Type> >
527 Foam::isoSurfaceCell::interpolate
528 (
529  const Field<Type>& cCoords,
530  const Field<Type>& pCoords
531 ) const
532 {
533  DynamicList<Type> triPoints(3*nCutCells_);
534  DynamicList<label> triMeshCells(nCutCells_);
535 
536  // Dummy snap data
537  DynamicList<Type> snappedPoints;
538  labelList snappedCc(mesh_.nCells(), -1);
539  labelList snappedPoint(mesh_.nPoints(), -1);
540 
541  generateTriPoints
542  (
543  cVals_,
544  pVals_,
545 
546  cCoords,
547  pCoords,
548 
549  snappedPoints,
550  snappedCc,
551  snappedPoint,
552 
553  triPoints,
554  triMeshCells
555  );
556 
557  return isoSurface::interpolate
558  (
559  points().size(),
560  triPointMergeMap_,
561  interpolatedPoints_,
562  interpolatedOldPoints_,
563  interpolationWeights_,
564  triPoints
565  );
566 }
567 
568 
569 // ************************************************************************* //
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:406
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:118
Foam::findIndex
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
Foam::DynamicList
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:56
isoSurfaceCell.H
Foam::tetMatcher
A cellMatcher for tet cells.
Definition: tetMatcher.H:51
Foam::MULES::interpolate
tmp< surfaceScalarField > interpolate(const RhoType &rho)
Definition: IMULESTemplates.C:40
Foam::isoSurfaceCell::generateTriPoints
void generateTriPoints(const DynamicList< Type > &snapped, const scalar s0, const Type &p0, const label p0Index, const scalar s1, const Type &p1, const label p1Index, const scalar s2, const Type &p2, const label p2Index, const scalar s3, const Type &p3, const label p3Index, DynamicList< Type > &points) const
Definition: isoSurfaceCellTemplates.C:77
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Foam::mag
dimensioned< scalar > mag(const dimensioned< Type > &)
polyMesh.H
Foam::triPoints
Triangle storage. Null constructable (unfortunately triangle<point, point> is not)
Definition: triPoints.H:49
Foam::isoSurfaceCell::generatePoint
Type generatePoint(const DynamicList< Type > &snappedPoints, const scalar s0, const Type &p0, const label p0Index, const scalar s1, const Type &p1, const label p1Index) const
Generate single point by interpolation or snapping.
Definition: isoSurfaceCellTemplates.C:35
Foam::label
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Foam::Field
Pre-declare SubField and related Field type.
Definition: Field.H:57
f1
scalar f1
Definition: createFields.H:28
Foam::DynamicList::shrink
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:258
s
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().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))
Foam::DynamicList::append
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
isoSurface.H
Foam::triFace
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
f
labelList f(nPoints)
Foam::List
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
Foam::tetMatcher::isA
virtual bool isA(const primitiveMesh &mesh, const label cellI)
Exact match. Uses faceSizeMatch.
Definition: tetMatcher.C:218
Foam::isoSurfaceCell::interpolate
tmp< Field< Type > > interpolate(const Field< Type > &cCoords, const Field< Type > &pCoords) const
Interpolates cCoords,pCoords.
Foam::FixedList::size
label size() const
Return the number of elements in the FixedList.
Definition: FixedListI.H:408
points
const pointField & points
Definition: gmvOutputHeader.H:1
tetMatcher.H
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Foam::Swap
void Swap(T &a, T &b)
Definition: Swap.H:43
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:259
Foam::cell
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.17