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-2016 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  WarningIn("Foam::isoSurfaceCell::generateTriPoints(..)")
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>
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  triPoints
562  );
563 }
564 
565 
566 // ************************************************************************* //
#define WarningIn(functionName)
Report a warning using Foam::Warning.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
const Field< point > & points() const
Return reference to global points.
const cellList & cells() const
label nCells() const
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))
void Swap(T &a, T &b)
Definition: Swap.H:43
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
Pre-declare SubField and related Field type.
Definition: Field.H:57
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
face triFace(3)
List< label > labelList
A List of labels.
Definition: labelList.H:56
const labelList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:818
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
label size() const
Return the number of elements in the UList.
tmp< Field< Type > > interpolate(const Field< Type > &cCoords, const Field< Type > &pCoords) const
Interpolates cCoords,pCoords.
label nPoints() const
dimensioned< scalar > mag(const dimensioned< Type > &)
A class for managing temporary objects.
Definition: PtrList.H:54
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1017
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004