FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/polyMesh.H>
28 #include <OpenFOAM/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 s0,
80  const Type& p0,
81  const label p0Index,
82 
83  const scalar s1,
84  const Type& p1,
85  const label p1Index,
86 
87  const scalar s2,
88  const Type& p2,
89  const label p2Index,
90 
91  const scalar s3,
92  const Type& p3,
93  const label p3Index,
94 
95  DynamicList<Type>& points
96 ) const
97 {
98  int triIndex = 0;
99  if (s0 < iso_)
100  {
101  triIndex |= 1;
102  }
103  if (s1 < iso_)
104  {
105  triIndex |= 2;
106  }
107  if (s2 < iso_)
108  {
109  triIndex |= 4;
110  }
111  if (s3 < iso_)
112  {
113  triIndex |= 8;
114  }
115 
116  /* Form the vertices of the triangles for each case */
117  switch (triIndex)
118  {
119  case 0x00:
120  case 0x0F:
121  break;
122 
123  case 0x0E:
124  case 0x01:
125  points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
126  points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
127  points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
128  break;
129 
130  case 0x0D:
131  case 0x02:
132  points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
133  points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
134  points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
135  break;
136 
137  case 0x0C:
138  case 0x03:
139  {
140  Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
141  Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
142 
143  points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
144  points.append(tp1);
145  points.append(tp2);
146  points.append(tp2);
147  points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
148  points.append(tp1);
149  }
150  break;
151 
152  case 0x0B:
153  case 0x04:
154  {
155  points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
156  points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
157  points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
158  }
159  break;
160 
161  case 0x0A:
162  case 0x05:
163  {
164  Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
165  Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
166 
167  points.append(tp0);
168  points.append(tp1);
169  points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
170  points.append(tp0);
171  points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
172  points.append(tp1);
173  }
174  break;
175 
176  case 0x09:
177  case 0x06:
178  {
179  Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
180  Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
181 
182  points.append(tp0);
183  points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
184  points.append(tp1);
185  points.append(tp0);
186  points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
187  points.append(tp1);
188  }
189  break;
190 
191  case 0x07:
192  case 0x08:
193  points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
194  points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
195  points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
196  break;
197  }
198 }
199 
200 
201 template<class Type>
202 void Foam::isoSurfaceCell::generateTriPoints
203 (
204  const scalarField& cVals,
205  const scalarField& pVals,
206 
207  const Field<Type>& cCoords,
208  const Field<Type>& pCoords,
209 
210  const DynamicList<Type>& snappedPoints,
211  const labelList& snappedCc,
212  const labelList& snappedPoint,
213 
214  DynamicList<Type>& triPoints,
215  DynamicList<label>& triMeshCells
216 ) const
217 {
218  tetMatcher tet;
219 
220  forAll(mesh_.cells(), cellI)
221  {
222  if (cellCutType_[cellI] != NOTCUT)
223  {
224  label oldNPoints = triPoints.size();
225 
226  const cell& cFaces = mesh_.cells()[cellI];
227 
228  if (tet.isA(mesh_, cellI))
229  {
230  // For tets don't do cell-centre decomposition, just use the
231  // tet points and values
232 
233  const face& f0 = mesh_.faces()[cFaces[0]];
234 
235  // Get the other point
236  const face& f1 = mesh_.faces()[cFaces[1]];
237  label oppositeI = -1;
238  forAll(f1, fp)
239  {
240  oppositeI = f1[fp];
241 
242  if (findIndex(f0, oppositeI) == -1)
243  {
244  break;
245  }
246  }
247 
248  generateTriPoints
249  (
250  snappedPoints,
251  pVals[f0[0]],
252  pCoords[f0[0]],
253  snappedPoint[f0[0]],
254 
255  pVals[f0[1]],
256  pCoords[f0[1]],
257  snappedPoint[f0[1]],
258 
259  pVals[f0[2]],
260  pCoords[f0[2]],
261  snappedPoint[f0[2]],
262 
263  pVals[oppositeI],
264  pCoords[oppositeI],
265  snappedPoint[oppositeI],
266 
267  triPoints
268  );
269  }
270  else
271  {
272  const cell& cFaces = mesh_.cells()[cellI];
273 
274  forAll(cFaces, cFaceI)
275  {
276  label faceI = cFaces[cFaceI];
277  const face& f = mesh_.faces()[faceI];
278 
279  for (label fp = 1; fp < f.size() - 1; fp++)
280  {
281  triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
282  //List<triFace> tris(triangulate(f));
283  //forAll(tris, i)
284  //{
285  // const triFace& tri = tris[i];
286 
287 
288  generateTriPoints
289  (
290  snappedPoints,
291 
292  pVals[tri[0]],
293  pCoords[tri[0]],
294  snappedPoint[tri[0]],
295 
296  pVals[tri[1]],
297  pCoords[tri[1]],
298  snappedPoint[tri[1]],
299 
300  pVals[tri[2]],
301  pCoords[tri[2]],
302  snappedPoint[tri[2]],
303 
304  cVals[cellI],
305  cCoords[cellI],
306  snappedCc[cellI],
307 
308  triPoints
309  );
310  }
311  }
312  }
313 
314 
315  // Every three triPoints is a cell
316  label nCells = (triPoints.size()-oldNPoints)/3;
317  for (label i = 0; i < nCells; i++)
318  {
319  triMeshCells.append(cellI);
320  }
321  }
322  }
323 
324  triPoints.shrink();
325  triMeshCells.shrink();
326 }
327 
328 
329 template <class Type>
332 (
333  const scalarField& cVals,
334  const scalarField& pVals,
335  const Field<Type>& cCoords,
336  const Field<Type>& pCoords
337 ) const
338 {
339  DynamicList<Type> triPoints(nCutCells_);
340  DynamicList<label> triMeshCells(nCutCells_);
341 
342  // Dummy snap data
343  DynamicList<Type> snappedPoints;
344  labelList snappedCc(mesh_.nCells(), -1);
345  labelList snappedPoint(mesh_.nPoints(), -1);
346 
347 
348  generateTriPoints
349  (
350  cVals,
351  pVals,
352 
353  cCoords,
354  pCoords,
355 
356  snappedPoints,
357  snappedCc,
358  snappedPoint,
359 
360  triPoints,
361  triMeshCells
362  );
363 
364 
365  // One value per point
366  tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
367  Field<Type>& values = tvalues();
368 
369  forAll(triPoints, i)
370  {
371  label mergedPointI = triPointMergeMap_[i];
372 
373  if (mergedPointI >= 0)
374  {
375  values[mergedPointI] = triPoints[i];
376  }
377  }
378 
379  return tvalues;
380 }
381 
382 
383 // ************************ vim: set sw=4 sts=4 et: ************************ //