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 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
32 template<class Type>
33 Type Foam::isoSurfaceCell::generatePoint
34 (
35  const DynamicList<Type>& snappedPoints,
36 
37  const scalar s0,
38  const Type& p0,
39  const label p0Index,
40 
41  const scalar s1,
42  const Type& p1,
43  const label p1Index
44 ) const
45 {
46  scalar d = s1-s0;
47 
48  if (mag(d) > VSMALL)
49  {
50  scalar s = (iso_-s0)/d;
51 
52  if (s >= 0.5 && s <= 1 && p1Index != -1)
53  {
54  return snappedPoints[p1Index];
55  }
56  else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
57  {
58  return snappedPoints[p0Index];
59  }
60  else
61  {
62  return s*p1 + (1.0-s)*p0;
63  }
64  }
65  else
66  {
67  scalar s = 0.4999;
68 
69  return s*p1 + (1.0-s)*p0;
70  }
71 }
72 
73 
74 template<class Type>
75 void Foam::isoSurfaceCell::generateTriPoints
76 (
77  const DynamicList<Type>& snapped,
78 
79  const scalar isoVal0,
80  const scalar s0,
81  const Type& p0,
82  const label p0Index,
83 
84  const scalar isoVal1,
85  const scalar s1,
86  const Type& p1,
87  const label p1Index,
88 
89  const scalar isoVal2,
90  const scalar s2,
91  const Type& p2,
92  const label p2Index,
93 
94  const scalar isoVal3,
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 0x0E:
128  case 0x01:
129  {
130  // 0 is common point. Orient such that normal points in positive
131  // gradient direction
132  if (isoVal0 >= isoVal1)
133  {
134  pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
135  pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
136  pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
137  }
138  else
139  {
140  pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
141  pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
142  pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
143  }
144  }
145  break;
146 
147  case 0x0D:
148  case 0x02:
149  {
150  // 1 is common point
151  if (isoVal1 >= isoVal0)
152  {
153  pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
154  pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
155  pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
156  }
157  else
158  {
159  pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
160  pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
161  pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
162  }
163  }
164  break;
165 
166  case 0x0C:
167  case 0x03:
168  {
169  Type s02 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
170  Type s13 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
171 
172  if (isoVal0 >= isoVal3)
173  {
174  pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
175  pts.append(s02);
176  pts.append(s13);
177  pts.append(s13);
178  pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
179  pts.append(s02);
180  }
181  else
182  {
183  pts.append(s02);
184  pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
185  pts.append(s13);
186  pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
187  pts.append(s13);
188  pts.append(s02);
189  }
190  }
191  break;
192 
193  case 0x0B:
194  case 0x04:
195  {
196  // 2 is common point
197  if (isoVal2 >= isoVal0)
198  {
199  pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
200  pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
201  pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
202  }
203  else
204  {
205  pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
206  pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
207  pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
208  }
209  }
210  break;
211 
212  case 0x0A:
213  case 0x05:
214  {
215  Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
216  Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
217 
218  if (isoVal3 >= isoVal0)
219  {
220  pts.append(s01);
221  pts.append(s23);
222  pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
223  pts.append(s01);
224  pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
225  pts.append(s23);
226  }
227  else
228  {
229  pts.append(s23);
230  pts.append(s01);
231  pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
232  pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
233  pts.append(s01);
234  pts.append(s23);
235  }
236  }
237  break;
238 
239  case 0x09:
240  case 0x06:
241  {
242  Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
243  Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
244 
245  if (isoVal3 >= isoVal1)
246  {
247  pts.append(s01);
248  pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
249  pts.append(s23);
250  pts.append(s01);
251  pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
252  pts.append(s23);
253  }
254  else
255  {
256  pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
257  pts.append(s01);
258  pts.append(s23);
259  pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
260  pts.append(s01);
261  pts.append(s23);
262  }
263  }
264  break;
265 
266  case 0x07:
267  case 0x08:
268  {
269  // 3 is common point
270  if (isoVal3 >= isoVal0)
271  {
272  pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
273  pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
274  pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
275  }
276  else
277  {
278  pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
279  pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
280  pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
281  }
282  }
283  break;
284  }
285 }
286 
287 
288 template<class Type>
289 void Foam::isoSurfaceCell::generateTriPoints
290 (
291  const scalarField& cVals,
292  const scalarField& pVals,
293 
294  const Field<Type>& cCoords,
295  const Field<Type>& pCoords,
296 
297  const DynamicList<Type>& snappedPoints,
298  const labelList& snappedCc,
299  const labelList& snappedPoint,
300 
301  DynamicList<Type>& triPoints,
302  DynamicList<label>& triMeshCells
303 ) const
304 {
305  tetMatcher tet;
306  label countNotFoundTets = 0;
307 
308  forAll(mesh_.cells(), cellI)
309  {
310  if (cellCutType_[cellI] != NOTCUT)
311  {
312  label oldNPoints = triPoints.size();
313 
314  const cell& cFaces = mesh_.cells()[cellI];
315 
316  if (tet.isA(mesh_, cellI))
317  {
318  // For tets don't do cell-centre decomposition, just use the
319  // tet points and values
320 
321  const face& f0 = mesh_.faces()[cFaces[0]];
322 
323  // Get the other point
324  const face& f1 = mesh_.faces()[cFaces[1]];
325  label oppositeI = -1;
326  forAll(f1, fp)
327  {
328  oppositeI = f1[fp];
329 
330  if (findIndex(f0, oppositeI) == -1)
331  {
332  break;
333  }
334  }
335 
336  // Start off from positive volume tet to make sure we
337  // generate outwards pointing tets
338  if (mesh_.faceOwner()[cFaces[0]] == cellI)
339  {
340  generateTriPoints
341  (
342  snappedPoints,
343 
344  pVals_[f0[1]],
345  pVals[f0[1]],
346  pCoords[f0[1]],
347  snappedPoint[f0[1]],
348 
349  pVals_[f0[0]],
350  pVals[f0[0]],
351  pCoords[f0[0]],
352  snappedPoint[f0[0]],
353 
354  pVals_[f0[2]],
355  pVals[f0[2]],
356  pCoords[f0[2]],
357  snappedPoint[f0[2]],
358 
359  pVals_[oppositeI],
360  pVals[oppositeI],
361  pCoords[oppositeI],
362  snappedPoint[oppositeI],
363 
364  triPoints
365  );
366  }
367  else
368  {
369  generateTriPoints
370  (
371  snappedPoints,
372 
373  pVals_[f0[0]],
374  pVals[f0[0]],
375  pCoords[f0[0]],
376  snappedPoint[f0[0]],
377 
378  pVals_[f0[1]],
379  pVals[f0[1]],
380  pCoords[f0[1]],
381  snappedPoint[f0[1]],
382 
383  pVals_[f0[2]],
384  pVals[f0[2]],
385  pCoords[f0[2]],
386  snappedPoint[f0[2]],
387 
388  pVals_[oppositeI],
389  pVals[oppositeI],
390  pCoords[oppositeI],
391  snappedPoint[oppositeI],
392 
393  triPoints
394  );
395  }
396  }
397  else
398  {
399  forAll(cFaces, cFaceI)
400  {
401  label faceI = cFaces[cFaceI];
402  const face& f = mesh_.faces()[faceI];
403 
404  label fp0 = mesh_.tetBasePtIs()[faceI];
405 
406  // Skip undefined tets
407  if (fp0 < 0)
408  {
409  fp0 = 0;
410  countNotFoundTets++;
411  }
412 
413  label fp = f.fcIndex(fp0);
414 
415  for (label i = 2; i < f.size(); i++)
416  {
417  label nextFp = f.fcIndex(fp);
418  triFace tri(f[fp0], f[fp], f[nextFp]);
419 
420  // Start off from positive volume tet to make sure we
421  // generate outwards pointing tets
422  if (mesh_.faceOwner()[faceI] == cellI)
423  {
424  generateTriPoints
425  (
426  snappedPoints,
427 
428  pVals_[tri[1]],
429  pVals[tri[1]],
430  pCoords[tri[1]],
431  snappedPoint[tri[1]],
432 
433  pVals_[tri[0]],
434  pVals[tri[0]],
435  pCoords[tri[0]],
436  snappedPoint[tri[0]],
437 
438  pVals_[tri[2]],
439  pVals[tri[2]],
440  pCoords[tri[2]],
441  snappedPoint[tri[2]],
442 
443  cVals_[cellI],
444  cVals[cellI],
445  cCoords[cellI],
446  snappedCc[cellI],
447 
448  triPoints
449  );
450  }
451  else
452  {
453  generateTriPoints
454  (
455  snappedPoints,
456 
457  pVals_[tri[0]],
458  pVals[tri[0]],
459  pCoords[tri[0]],
460  snappedPoint[tri[0]],
461 
462  pVals_[tri[1]],
463  pVals[tri[1]],
464  pCoords[tri[1]],
465  snappedPoint[tri[1]],
466 
467  pVals_[tri[2]],
468  pVals[tri[2]],
469  pCoords[tri[2]],
470  snappedPoint[tri[2]],
471 
472  cVals_[cellI],
473  cVals[cellI],
474  cCoords[cellI],
475  snappedCc[cellI],
476 
477  triPoints
478  );
479  }
480 
481  fp = nextFp;
482  }
483  }
484  }
485 
486 
487  // Every three triPoints is a cell
488  label nCells = (triPoints.size()-oldNPoints)/3;
489  for (label i = 0; i < nCells; i++)
490  {
491  triMeshCells.append(cellI);
492  }
493  }
494  }
495 
496  if (countNotFoundTets > 0)
497  {
498  WarningIn("Foam::isoSurfaceCell::generateTriPoints")
499  << "Could not find " << countNotFoundTets
500  << " tet base points, which may lead to inverted triangles."
501  << endl;
502  }
503 
504  triPoints.shrink();
505  triMeshCells.shrink();
506 }
507 
508 
509 template<class Type>
512 (
513  const scalarField& cVals,
514  const scalarField& pVals,
515  const Field<Type>& cCoords,
516  const Field<Type>& pCoords
517 ) const
518 {
519  DynamicList<Type> triPoints(nCutCells_);
520  DynamicList<label> triMeshCells(nCutCells_);
521 
522  // Dummy snap data
523  DynamicList<Type> snappedPoints;
524  labelList snappedCc(mesh_.nCells(), -1);
525  labelList snappedPoint(mesh_.nPoints(), -1);
526 
527 
528  generateTriPoints
529  (
530  cVals,
531  pVals,
532 
533  cCoords,
534  pCoords,
535 
536  snappedPoints,
537  snappedCc,
538  snappedPoint,
539 
540  triPoints,
541  triMeshCells
542  );
543 
544 
545  // One value per point
546  tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
547  Field<Type>& values = tvalues();
548 
549  forAll(triPoints, i)
550  {
551  label mergedPointI = triPointMergeMap_[i];
552 
553  if (mergedPointI >= 0)
554  {
555  values[mergedPointI] = triPoints[i];
556  }
557  }
558 
559  return tvalues;
560 }
561 
562 
563 template<class Type>
566 (
567  const Field<Type>& cCoords,
568  const Field<Type>& pCoords
569 ) const
570 {
571  DynamicList<Type> triPoints(nCutCells_);
572  DynamicList<label> triMeshCells(nCutCells_);
573 
574  // Dummy snap data
575  DynamicList<Type> snappedPoints;
576  labelList snappedCc(mesh_.nCells(), -1);
577  labelList snappedPoint(mesh_.nPoints(), -1);
578 
579 
580  generateTriPoints
581  (
582  cVals_,
583  pVals_,
584 
585  cCoords,
586  pCoords,
587 
588  snappedPoints,
589  snappedCc,
590  snappedPoint,
591 
592  triPoints,
593  triMeshCells
594  );
595 
596 
597  // One value per point
598  tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
599  Field<Type>& values = tvalues();
600 
601  forAll(triPoints, i)
602  {
603  label mergedPointI = triPointMergeMap_[i];
604 
605  if (mergedPointI >= 0)
606  {
607  values[mergedPointI] = triPoints[i];
608  }
609  }
610 
611  return tvalues;
612 }
613 
614 
615 // ************************************************************************* //
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
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 ))
dimensioned< scalar > mag(const dimensioned< Type > &)
const labelList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:864
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
label nCells() const
const cellList & cells() const
face triFace(3)
tmp< Field< Type > > interpolate(const scalarField &cVals, const scalarField &pVals, const Field< Type > &cCoords, const Field< Type > &pCoords) const
Interpolates cCoords,pCoords. Takes the original fields.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define WarningIn(functionName)
Report a warning using Foam::Warning.
#define forAll(list, i)
Definition: UList.H:421
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1073
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
List< label > labelList
A List of labels.
Definition: labelList.H:56
label size() const
Return the number of elements in the UList.
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1060
A class for managing temporary objects.
Definition: PtrList.H:118
const Field< point > & points() const
Return reference to global points.
label nPoints() const