60 for (i=ts.size()-1; i>=0; i--)
73 for (
label j=i; j>=0; j--)
79 Pout<<
" triSurface::triSurfInstance(const Time& d)" 80 <<
"reading " << foamName
81 <<
" from " << ts[j].
name()/typeName
92 Pout<<
" triSurface::triSurfInstance(const Time& d)" 93 <<
"reading " << foamName
94 <<
" from constant/" <<
endl;
103 const label defaultRegion
110 const face& f = faces[facei];
115 <<
"Face at position " << facei
116 <<
" does not have three vertices:" << f
125 tri.
region() = defaultRegion;
135 const label defaultRegion
142 const triFace& f = faces[facei];
149 tri.
region() = defaultRegion;
158 void Foam::triSurface::printTriangle
167 << pre.c_str() <<
"vertex numbers:" 168 << f[0] <<
' ' << f[1] <<
' ' << f[2] <<
endl 169 << pre.c_str() <<
"vertex coords :" 170 << points[f[0]] <<
' ' << points[f[1]] <<
' ' << points[f[2]]
171 << pre.c_str() <<
"region :" << f.
region() <<
endl 183 while ((line.empty() || line[0] ==
'#') && is.
good());
189 Foam::scalar Foam::triSurface::pointNormalWeight
205 const vector e1 = points[f[index]] - points[f[f.
fcIndex(index)]];
206 const vector e2 = points[f[index]] - points[f[f.rcIndex(index)]];
225 const labelList& meshPoints = this->meshPoints();
233 const label facei = pFaces[fi];
234 const triFace& f = operator[](facei);
237 const scalar weight =
238 pointNormalWeight(f, meshPoints[pi], fa, points);
240 pointNormals[
pi] += weight*fa;
243 pointNormals[
pi] /=
mag(pointNormals[pi]) + vSmall;
246 return tpointNormals;
259 const labelList& meshPoints = this->meshPoints();
263 const point& pt = points[meshPoints[
pi]];
264 const vector& normal = pointNormals[
pi];
266 if (
mag(normal) < small)
278 vector dir2 = dir1 ^ normal;
281 pointCoordSys[
pi] =
triad(dir1, dir2, normal);
284 return pointCoordSys;
288 bool Foam::triSurface::read(
Istream& is)
290 is >> patches_ >> storedPoints() >> storedFaces();
296 bool Foam::triSurface::read
303 if (check && !
exists(name))
314 return read(unzipName, unzipName.
ext(),
false);
316 else if (ext ==
"ftr")
320 else if (ext ==
"stl")
322 return readSTL(name);
324 else if (ext ==
"stlb")
326 return readSTLBINARY(name);
328 else if (ext ==
"gts")
330 return readGTS(name);
332 else if (ext ==
"obj")
334 return readOBJ(name);
336 else if (ext ==
"off")
338 return readOFF(name);
340 else if (ext ==
"tri")
342 return readTRI(name);
344 else if (ext ==
"ac")
348 else if (ext ==
"nas")
350 return readNAS(name);
352 else if (ext ==
"vtk")
354 return readVTK(name);
359 <<
"unknown file extension " << ext
360 <<
". Supported extensions are '.ftr', '.stl', '.stlb', '.gts'" 361 <<
", '.obj', '.ac', '.off', '.nas', '.tri' and '.vtk'" 369 void Foam::triSurface::write
380 else if (ext ==
"stl")
382 return writeSTLASCII(sort,
OFstream(name)());
384 else if (ext ==
"stlb")
386 ofstream outFile(name.c_str(), std::ios::binary);
388 writeSTLBINARY(outFile);
390 else if (ext ==
"gts")
392 return writeGTS(sort,
OFstream(name)());
394 else if (ext ==
"obj")
398 else if (ext ==
"off")
402 else if (ext ==
"vtk")
406 else if (ext ==
"tri")
410 else if (ext ==
"ac")
414 else if (ext ==
"smesh")
421 <<
"unknown file extension " << ext
422 <<
" for file " << name
423 <<
". Supported extensions are '.ftr', '.stl', '.stlb', " 424 <<
"'.gts', '.obj', '.vtk'" 425 <<
", '.off', '.dx', '.smesh', '.ac' and '.tri'" 436 forAll(sortedRegion, facei)
438 sortedRegion[facei] = operator[](facei).region();
442 faceMap = sortedRegion.
indices();
452 operator[](faceMap.
last()).region()
462 label region = operator[](facei).region();
464 newPatches[region].size()++;
470 label startFacei = 0;
471 forAll(newPatches, newPatchi)
475 newPatch.
index() = newPatchi;
477 label oldPatchi = newPatchi;
480 newPatch.
start() = startFacei;
484 if ((oldPatchi < patches_.size()) && (patches_[oldPatchi].
name() !=
""))
486 newPatch.
name() = patches_[oldPatchi].name();
495 (oldPatchi < patches_.size())
496 && (patches_[oldPatchi].geometricType() !=
"")
499 newPatch.
geometricType() = patches_[oldPatchi].geometricType();
506 startFacei += newPatch.
size();
513 void Foam::triSurface::setDefaultPatches()
521 patches_.setSize(newPatches.
size());
527 patches_[
patchi].geometricType() = newPatches[
patchi].geometricType();
538 sortedEdgeFacesPtr_(nullptr),
539 edgeOwnerPtr_(nullptr)
553 sortedEdgeFacesPtr_(
nullptr),
554 edgeOwnerPtr_(
nullptr)
568 sortedEdgeFacesPtr_(
nullptr),
569 edgeOwnerPtr_(
nullptr)
582 sortedEdgeFacesPtr_(
nullptr),
583 edgeOwnerPtr_(
nullptr)
595 sortedEdgeFacesPtr_(
nullptr),
596 edgeOwnerPtr_(
nullptr)
608 ParentType(convertToTri(triangles, 0), points),
610 sortedEdgeFacesPtr_(
nullptr),
611 edgeOwnerPtr_(
nullptr)
621 sortedEdgeFacesPtr_(nullptr),
622 edgeOwnerPtr_(nullptr)
636 sortedEdgeFacesPtr_(nullptr),
637 edgeOwnerPtr_(nullptr)
649 sortedEdgeFacesPtr_(nullptr),
650 edgeOwnerPtr_(nullptr)
667 patches_(ts.patches()),
668 sortedEdgeFacesPtr_(nullptr),
669 edgeOwnerPtr_(nullptr)
676 patches_(move(ts.patches())),
677 sortedEdgeFacesPtr_(nullptr),
678 edgeOwnerPtr_(nullptr)
717 if (!sortedEdgeFacesPtr_)
719 calcSortedEdgeFaces();
722 return *sortedEdgeFacesPtr_;
733 return *edgeOwnerPtr_;
753 if (scaleFactor > 0 && scaleFactor != 1.0)
777 if (f[fp] < 0 || f[fp] > maxPointi)
781 <<
" uses point indices outside point range 0.." 795 bool hasInvalid =
false;
801 if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
804 valid[facei] =
false;
810 <<
"triangle " << facei
811 <<
" does not have three unique vertices:\n";
829 label neighbour = eFaces[i];
831 if (neighbour > facei)
838 ((f[0] == n[0]) || (f[0] == n[1]) || (f[0] == n[2]))
839 && ((f[1] == n[0]) || (f[1] == n[1]) || (f[1] == n[2]))
840 && ((f[2] == n[0]) || (f[2] == n[1]) || (f[2] == n[2]))
843 valid[facei] =
false;
849 <<
"triangles share the same vertices:\n" 850 <<
" face 1 :" << facei <<
endl;
856 << neighbour <<
endl;
877 (*this)[newFacei++] =
f;
884 <<
"Removing " <<
size() - newFacei
885 <<
" illegal faces." <<
endl;
887 (*this).setSize(newFacei);
901 const labelList& myFaces = eFaces[edgeI];
906 <<
"Edge " << edgeI <<
" with vertices " <<
edges()[edgeI]
907 <<
" has no edgeFaces" 910 else if (myFaces.
size() > 2 && verbose)
913 <<
"Edge " << edgeI <<
" with vertices " <<
edges()[edgeI]
914 <<
" has more than 2 faces connected to it : " << myFaces
924 stitchTriangles(small, verbose);
939 const label currentZone,
953 label facei = changedFaces[i];
959 label edgeI = fEdges[i];
961 if (!borderEdge[edgeI])
967 label nbrFacei = eFaces[j];
969 if (faceZone[nbrFacei] == -1)
971 faceZone[nbrFacei] = currentZone;
972 newChangedFaces.
append(nbrFacei);
974 else if (faceZone[nbrFacei] != currentZone)
977 <<
"Zones " << faceZone[nbrFacei]
978 <<
" at face " << nbrFacei
979 <<
" connects to zone " << currentZone
980 <<
" at face " << facei
988 if (newChangedFaces.empty())
993 changedFaces.
transfer(newChangedFaces);
1010 <<
"borderEdge boolList not same size as number of edges" <<
endl 1011 <<
"borderEdge:" << borderEdge.
size() <<
endl 1012 <<
"nEdges :" <<
nEdges()
1018 for (
label startFacei = 0;; zoneI++)
1021 for (; startFacei <
size(); startFacei++)
1023 if (faceZone[startFacei] == -1)
1029 if (startFacei >=
size())
1034 faceZone[startFacei] = zoneI;
1036 markZone(borderEdge, startFacei, zoneI, faceZone);
1061 forAll(include, oldFacei)
1063 if (include[oldFacei])
1066 faceMap[facei++] = oldFacei;
1074 if (!pointHad[labI])
1076 pointHad[labI] =
true;
1077 pointMap[pointi++] = labI;
1108 newPoints[pointi] = locPoints[pointMap[pointi]];
1109 oldToNew[pointMap[pointi]] = pointi;
1118 const labelledTri& tri = locFaces[faceMap[facei]];
1120 newTriangles[facei][0] = oldToNew[tri[0]];
1121 newTriangles[facei][1] = oldToNew[tri[1]];
1122 newTriangles[facei][2] = oldToNew[tri[2]];
1123 newTriangles[facei].region() = tri.region();
1137 faces[facei] =
operator[](facei).triFaceFace();
1156 this->weightedPointNormals()
1161 this->pointCoordSys(pointNormals)
1181 const edge&
e = fEdges[feI];
1183 edgeVectors[feI] = e.
vec(points);
1184 normalDifferences[feI] =
1185 pointNormals[meshPointMap[e[0]]]
1186 - pointNormals[meshPointMap[e[1]]];
1190 const vector& e0 = edgeVectors[0];
1192 const vector e1 = (e0 ^ eN);
1194 if (
magSqr(eN) < rootVSmall)
1199 triad faceCoordSys(e0, e1, eN);
1207 for (
label i = 0; i < 3; ++i)
1209 const scalar
x = edgeVectors[i] & faceCoordSys[0];
1210 const scalar
y = edgeVectors[i] & faceCoordSys[1];
1218 const scalar dndx = normalDifferences[i] & faceCoordSys[0];
1219 const scalar dndy = normalDifferences[i] & faceCoordSys[1];
1222 Z[1] += dndx*y + dndy*
x;
1229 const symmTensor2D secondFundamentalTensor(Z[0], Z[1], Z[2]);
1235 const label patchPointIndex = meshPointMap[f[fpi]];
1237 const triad& ptCoordSys = pointCoordSys[patchPointIndex];
1239 if (!ptCoordSys.
set())
1250 const triad rotatedFaceCoordSys = rotTensor &
tensor(faceCoordSys);
1256 ptCoordSys[0] & rotatedFaceCoordSys[0],
1257 ptCoordSys[0] & rotatedFaceCoordSys[1]
1261 ptCoordSys[1] & rotatedFaceCoordSys[0],
1262 ptCoordSys[1] & rotatedFaceCoordSys[1]
1273 projTensor.
x() & (secondFundamentalTensor & projTensor.
x()),
1274 projTensor.
x() & (secondFundamentalTensor & projTensor.
y()),
1275 projTensor.
y() & (secondFundamentalTensor & projTensor.
y())
1279 const scalar weight = pointNormalWeight
1282 meshPoints[patchPointIndex],
1288 pointFundamentalTensors[patchPointIndex] +=
1289 weight*projectedFundamentalTensor;
1291 accumulatedWeights[patchPointIndex] += weight;
1295 forAll(curvaturePointField, pi)
1297 pointFundamentalTensors[
pi] /= (accumulatedWeights[
pi] + small);
1299 const vector2D principalCurvatures =
1306 mag(principalCurvatures[0]),
1307 mag(principalCurvatures[1])
1314 return tcurvaturePointField;
1318 void Foam::triSurface::write
1321 const bool sortByRegion
1324 write(name, name.
ext(), sortByRegion);
1337 os.
check(
"triSurface::write(Ostream&)");
1341 void Foam::triSurface::write(
const Time& d)
const 1368 label pointi = f[fp];
1369 if (pointIsUsed.
set(pointi, 1))
1378 os <<
"Triangles : " <<
size() <<
endl 1379 <<
"Vertices : " << nPoints <<
endl 1380 <<
"Bounding Box : " << bb <<
endl;
1400 patches_ = move(ts.patches());
virtual const fileName & name() const
Return the name of the stream.
label nPoints() const
Return number of points supporting patch faces.
void cleanup(const bool verbose)
Remove non-valid triangles.
bool exists(const fileName &, const bool checkVariants=true, const bool followLink=true)
Does the name exist (as directory or file) in the file system?
#define forAll(list, i)
Loop across all elements in list.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
label start() const
Return start label of this patch in the polyMesh face list.
virtual void movePoints(const Field< PointType > &)
Correct patch after moving points.
void sort()
(stable) sort the list (if changed after construction time)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
A class for handling file names.
const word & name() const
Return name.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A face is a list of labels corresponding to mesh vertices.
void operator=(const triSurface &)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
T & operator[](const label)
Return element of UList.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual bool check(const char *operation) const
Check IOstream status for given operation.
A list that is sorted upon construction or when explicitly requested with the sort() method...
T & ref() const
Return non-const reference or generate a fatal error.
bool isFile(const fileName &, const bool checkVariants=true, const bool followLink=true)
Does the name exist as a file in the file system?
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Vector2D< Cmpt > x() const
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
void size(const label)
Override size to be inconsistent with allocated storage.
const word & geometricType() const
Return the type of the patch.
dimensionedVector eigenValues(const dimensionedTensor &dt)
tmp< scalarField > curvature() const
Return the curvature of surface at the points.
Ostream & endl(Ostream &os)
Add newline and flush stream.
A bounding box defined in terms of the points at its extremities.
edgeList edges() const
Return edges in face point ordering,.
virtual void scalePoints(const scalar)
Scale points. A non-positive factor is ignored.
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
const Field< PointType > & localPoints() const
Return pointField of points in patch.
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from unit vector n1 to n2.
void clearPatchMeshAddr()
void clearPatchMeshAddr()
vector area(const pointField &) const
Return vector area.
Vector2D< Cmpt > y() const
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
bool good() const
Return true if next operation might succeed.
word ext() const
Return file name extension (part after last .)
Templated 2D tensor derived from VectorSpace adding construction from 4 components, element access using xx(), xy(), yx() and yy() member functions and the iner-product (dot-product) and outer-product of two Vector2Ds (tensor-product) operators.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
bool set(const direction d) const
Is the vector in the direction d set.
virtual void movePoints(const pointField &)
Move points.
Functions used by OpenFOAM that are specific to POSIX compliant operating systems and need to be repl...
label region() const
Return region label.
bool read(const char *, int32_t &)
static fileName triSurfInstance(const Time &)
Name of triSurface directory to use.
const Map< label > & meshPointMap() const
Mesh point map. Given the global point index find its.
vectorField pointField
pointField is a vectorField.
const labelListList & sortedEdgeFaces() const
Return edge-face addressing sorted (for edges with more than.
'Patch' on surface as subset of triSurface.
void set(const PackedList< 1 > &)
Set specified bits.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
void writeStats(Ostream &) const
Write some statistics.
point aPoint() const
Return a point on the plane.
A class for handling words, derived from string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
A triangular face using a FixedList of labels corresponding to mesh vertices.
virtual const fileName & name() const
Return the name of the stream.
void checkTriangles(const bool verbose)
Check/remove duplicate/degenerate triangles.
const fileName & caseName() const
Return case name.
void append(const T &)
Append an element at the end of the list.
const word & constant() const
Return constant name.
void checkEdges(const bool verbose)
Check triply (or more) connected edges.
word name() const
Return file name (part beyond last /)
const Field< PointType > & points() const
Return reference to global points.
void write(std::ostream &os, const bool binary, List< floatScalar > &fField)
Write floats ascii or binary.
const labelListList & edgeFaces() const
Return edge-face addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
Triangle with additional region number.
errorManip< error > abort(error &err)
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< small) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &mergedCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
pointField & storedPoints()
Non-const access to global points.
Field< triad > triadField
Specialisation of Field<T> for triad.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
vector vec(const pointField &) const
Return the vector (end - start)
label nEdges() const
Return number of edges in patch.
label size() const
Return the number of elements in the FixedList.
void markZone(const boolList &borderEdge, const label facei, const label currentZone, labelList &faceZone) const
Fill faceZone with currentZone for every face reachable.
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
A templated 2D square symmetric matrix of objects of <T>, where the n x n matrix dimension is known a...
Representation of a 3D Cartesian coordinate system as a Vector of vectors.
const Field< PointType > & pointNormals() const
Return point normals for patch.
void operator=(const UList< T > &)
Assignment to UList operator. Takes linear time.
void setSize(const label)
Reset size of List.
Template functions to aid in the implementation of demand driven data.
const point & max() const
Maximum describing the bounding box.
const geometricSurfacePatchList & patches() const
static const boundBox invertedBox
A very large inverted boundBox: min/max == +/- vGreat.
ISstream & getLine(string &, const bool continuation=true)
Read line into a string.
triSurface()
Construct null.
std::remove_reference< ::Foam::List< labelledTri > >::type::value_type FaceType
scalar timeOutputValue() const
Return current time value.
void normalise()
Normalise each set axis vector to have a unit magnitude.
label size() const
Return size of this patch in the polyMesh face list.
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
void subsetMeshMap(const boolList &include, labelList &pointMap, labelList &faceMap) const
'Create' sub mesh, including only faces for which
instantList times() const
Search the case for valid time directories.
Ostream & operator<<(Ostream &, const ensightPart &)
fileName lessExt() const
Return file name without extension (part before last .)
Templated 2D symmetric tensor derived from VectorSpace adding construction from 4 components...
const labelListList & faceEdges() const
Return face-edge addressing.
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
faceList faces() const
Return the list of triangles as a faceList.
const doubleScalar e
Elementary charge.
label index() const
Return the index of this patch in the boundaryMesh.
const point & min() const
Minimum describing the bounding box.
label markZones(const boolList &borderEdge, labelList &faceZone) const
(size and) fills faceZone with zone of face. Zone is area
void writeVTK(OFstream &os, const Type &value)
A class for managing temporary objects.
triSurface subsetMesh(const boolList &include, labelList &pointMap, labelList &faceMap) const
Return new surface. Returns pointMap, faceMap from.
T & last()
Return the last element of the list.
Triangulated surface description with patch information.
fileName path() const
Return path.
void deleteDemandDrivenData(DataPtr &dataPtr)
const labelList & edgeOwner() const
If 2 face neighbours: label of face where ordering of edge.
A class for handling character strings derived from std::string.
label size() const
Return the number of elements in the UList.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Tensor< scalar > tensor
Tensor of scalars.
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
Solve the matrix using LU decomposition with pivoting.