37 Type Foam::isoSurfaceCell::generatePoint
39 const DynamicList<Type>& snappedPoints,
50 const scalar d = s1-s0;
54 const scalar
s = (
iso_-s0)/d;
56 if (
s >= 0.5 &&
s <= 1 && p1Index != -1)
58 return snappedPoints[p1Index];
60 else if (
s >= 0.0 &&
s <= 0.5 && p0Index != -1)
62 return snappedPoints[p0Index];
66 return s*p1 + (1.0-
s)*
p0;
71 constexpr scalar
s = 0.4999;
73 return s*p1 + (1.0-
s)*
p0;
79 void Foam::isoSurfaceCell::generateTriPoints
81 const DynamicList<Type>& snapped,
99 DynamicList<Type>&
pts 130 pts.
append(generatePoint(snapped,s0,
p0,p0Index,s1,p1,p1Index));
131 pts.
append(generatePoint(snapped,s0,
p0,p0Index,s2,p2,p2Index));
132 pts.
append(generatePoint(snapped,s0,
p0,p0Index,s3,p3,p3Index));
133 if (triIndex == 0x0E)
137 std::swap(
pts[sz-2],
pts[sz-1]);
145 pts.
append(generatePoint(snapped,s1,p1,p1Index,s0,
p0,p0Index));
146 pts.
append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
147 pts.
append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
149 if (triIndex == 0x0D)
153 std::swap(
pts[sz-2],
pts[sz-1]);
161 Type p0p2 = generatePoint(snapped,s0,
p0,p0Index,s2,p2,p2Index);
162 Type p1p3 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
164 pts.
append(generatePoint(snapped,s0,
p0,p0Index,s3,p3,p3Index));
171 generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
175 if (triIndex == 0x0C)
179 std::swap(
pts[sz-5],
pts[sz-4]);
180 std::swap(
pts[sz-2],
pts[sz-1]);
188 pts.
append(generatePoint(snapped,s2,p2,p2Index,s0,
p0,p0Index));
189 pts.
append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
190 pts.
append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
192 if (triIndex == 0x0B)
196 std::swap(
pts[sz-2],
pts[sz-1]);
204 Type p0p1 = generatePoint(snapped,s0,
p0,p0Index,s1,p1,p1Index);
205 Type p2p3 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
209 pts.
append(generatePoint(snapped,s0,
p0,p0Index,s3,p3,p3Index));
212 pts.
append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
215 if (triIndex == 0x0A)
219 std::swap(
pts[sz-5],
pts[sz-4]);
220 std::swap(
pts[sz-2],
pts[sz-1]);
228 Type p0p1 = generatePoint(snapped,s0,
p0,p0Index,s1,p1,p1Index);
229 Type p2p3 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
232 pts.
append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
237 pts.
append(generatePoint(snapped,s0,
p0,p0Index,s2,p2,p2Index));
239 if (triIndex == 0x09)
243 std::swap(
pts[sz-5],
pts[sz-4]);
244 std::swap(
pts[sz-2],
pts[sz-1]);
252 pts.
append(generatePoint(snapped,s3,p3,p3Index,s0,
p0,p0Index));
253 pts.
append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
254 pts.
append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
255 if (triIndex == 0x07)
259 std::swap(
pts[sz-2],
pts[sz-1]);
268 void Foam::isoSurfaceCell::generateTriPoints
273 const Field<Type>& cCoords,
274 const Field<Type>& pCoords,
276 const DynamicList<Type>& snappedPoints,
280 DynamicList<Type>& triPoints,
281 DynamicList<label>& triMeshCells
284 label countNotFoundTets = 0;
286 forAll(mesh_.cells(), celli)
288 if ((cellCutType_[celli] & cutType::ANYCUT) != 0)
290 label oldNPoints = triPoints.size();
292 const cell& cFaces = mesh_.cells()[celli];
299 const face& f0 = mesh_.faces()[cFaces[0]];
302 const face& f1 = mesh_.faces()[cFaces[1]];
303 label oppositeI = -1;
308 if (!f0.found(oppositeI))
316 if (mesh_.faceOwner()[cFaces[0]] == celli)
336 snappedPoint[oppositeI],
361 snappedPoint[oppositeI],
371 label facei = cFaces[cFacei];
372 const face&
f = mesh_.faces()[facei];
374 label fp0 = mesh_.tetBasePtIs()[facei];
384 for (label i = 2; i <
f.
size(); i++)
391 if (mesh_.faceOwner()[facei] == celli)
399 snappedPoint[tri[1]],
403 snappedPoint[tri[0]],
407 snappedPoint[tri[2]],
424 snappedPoint[tri[0]],
428 snappedPoint[tri[1]],
432 snappedPoint[tri[2]],
449 label nCells = (triPoints.size()-oldNPoints)/3;
450 for (label i = 0; i < nCells; i++)
452 triMeshCells.append(celli);
457 if (countNotFoundTets > 0)
460 <<
"Could not find " << countNotFoundTets
461 <<
" tet base points, which may lead to inverted triangles." 466 triMeshCells.shrink();
472 Foam::isoSurfaceCell::interpolateTemplate
484 labelList snappedPoint(mesh_.nPoints(), -1);
502 return isoSurfacePoint::interpolate
507 interpolatedOldPoints_,
508 interpolationWeights_,
void size(const label n)
Older name for setAddressableSize.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void append(const T &val)
Append an element at the end of the list.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const scalar iso_
Iso value.
#define forAll(list, i)
Loop across all elements in list.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
Generic templated field type.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Triangle point storage. Default constructable (triangle is not)
#define WarningInFunction
Report a warning using Foam::Warning.
List< label > labelList
A List of labels.
A class for managing temporary objects.
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))
static bool test(const UList< face > &faces)
Test if given list of faces satisfies criteria for TET. (4 tri)
const volScalarField & p0