59 for (i=ts.size()-1; i>=0; i--)
72 for (
label j=i; j>=0; j--)
78 Pout<<
" triSurface::triSurfInstance(const Time& d)" 79 <<
"reading " << foamName
80 <<
" from " << ts[j].
name()/typeName
91 Pout<<
" triSurface::triSurfInstance(const Time& d)" 92 <<
"reading " << foamName
93 <<
" from constant/" <<
endl;
102 const label defaultRegion
109 const face& f = faces[facei];
114 <<
"Face at position " << facei
115 <<
" does not have three vertices:" << f
124 tri.
region() = defaultRegion;
134 const label defaultRegion
141 const triFace& f = faces[facei];
148 tri.
region() = defaultRegion;
157 void Foam::triSurface::printTriangle
166 << pre.c_str() <<
"vertex numbers:" 167 << f[0] <<
' ' << f[1] <<
' ' << f[2] <<
endl 168 << pre.c_str() <<
"vertex coords :" 169 << points[f[0]] <<
' ' << points[f[1]] <<
' ' << points[f[2]]
170 << pre.c_str() <<
"region :" << f.
region() <<
endl 182 while ((line.empty() || line[0] ==
'#') && is.
good());
188 Foam::scalar Foam::triSurface::pointNormalWeight
204 const vector e1 = points[f[index]] - points[f[f.
fcIndex(index)]];
205 const vector e2 = points[f[index]] - points[f[f.rcIndex(index)]];
224 const labelList& meshPoints = this->meshPoints();
232 const label facei = pFaces[fi];
233 const triFace& f = operator[](facei);
236 const scalar weight =
237 pointNormalWeight(f, meshPoints[pi], fa, points);
239 pointNormals[
pi] += weight*fa;
242 pointNormals[
pi] /=
mag(pointNormals[pi]) + vSmall;
245 return tpointNormals;
258 const labelList& meshPoints = this->meshPoints();
262 const point& pt = points[meshPoints[
pi]];
263 const vector& normal = pointNormals[
pi];
265 if (
mag(normal) < small)
277 vector dir2 = dir1 ^ normal;
280 pointCoordSys[
pi] =
triad(dir1, dir2, normal);
283 return pointCoordSys;
287 bool Foam::triSurface::read(
Istream& is)
289 is >> patches_ >> storedPoints() >> storedFaces();
295 bool Foam::triSurface::read
302 if (check && !
exists(name))
313 return read(unzipName, unzipName.
ext(),
false);
315 else if (ext ==
"ftr")
319 else if (ext ==
"stl")
321 return readSTL(name);
323 else if (ext ==
"stlb")
325 return readSTLBINARY(name);
327 else if (ext ==
"gts")
329 return readGTS(name);
331 else if (ext ==
"obj")
333 return readOBJ(name);
335 else if (ext ==
"off")
337 return readOFF(name);
339 else if (ext ==
"tri")
341 return readTRI(name);
343 else if (ext ==
"ac")
347 else if (ext ==
"nas")
349 return readNAS(name);
351 else if (ext ==
"vtk")
353 return readVTK(name);
358 <<
"unknown file extension " << ext
359 <<
". Supported extensions are '.ftr', '.stl', '.stlb', '.gts'" 360 <<
", '.obj', '.ac', '.off', '.nas', '.tri' and '.vtk'" 368 void Foam::triSurface::write
379 else if (ext ==
"stl")
381 return writeSTLASCII(sort,
OFstream(name)());
383 else if (ext ==
"stlb")
385 ofstream outFile(name.c_str(), std::ios::binary);
387 writeSTLBINARY(outFile);
389 else if (ext ==
"gts")
391 return writeGTS(sort,
OFstream(name)());
393 else if (ext ==
"obj")
397 else if (ext ==
"off")
401 else if (ext ==
"vtk")
405 else if (ext ==
"tri")
409 else if (ext ==
"ac")
413 else if (ext ==
"smesh")
420 <<
"unknown file extension " << ext
421 <<
" for file " << name
422 <<
". Supported extensions are '.ftr', '.stl', '.stlb', " 423 <<
"'.gts', '.obj', '.vtk'" 424 <<
", '.off', '.dx', '.smesh', '.ac' and '.tri'" 435 forAll(sortedRegion, facei)
437 sortedRegion[facei] = operator[](facei).region();
441 faceMap = sortedRegion.
indices();
451 operator[](faceMap.
last()).region()
461 label region = operator[](facei).region();
463 newPatches[region].size()++;
469 label startFacei = 0;
470 forAll(newPatches, newPatchi)
474 newPatch.
index() = newPatchi;
476 label oldPatchi = newPatchi;
479 newPatch.
start() = startFacei;
483 if ((oldPatchi < patches_.size()) && (patches_[oldPatchi].
name() !=
""))
485 newPatch.
name() = patches_[oldPatchi].name();
494 (oldPatchi < patches_.size())
495 && (patches_[oldPatchi].geometricType() !=
"")
498 newPatch.
geometricType() = patches_[oldPatchi].geometricType();
505 startFacei += newPatch.
size();
512 void Foam::triSurface::setDefaultPatches()
520 patches_.setSize(newPatches.
size());
526 patches_[
patchi].geometricType() = newPatches[
patchi].geometricType();
537 sortedEdgeFacesPtr_(nullptr),
538 edgeOwnerPtr_(nullptr)
552 sortedEdgeFacesPtr_(
nullptr),
553 edgeOwnerPtr_(
nullptr)
567 sortedEdgeFacesPtr_(
nullptr),
568 edgeOwnerPtr_(
nullptr)
581 sortedEdgeFacesPtr_(
nullptr),
582 edgeOwnerPtr_(
nullptr)
594 sortedEdgeFacesPtr_(
nullptr),
595 edgeOwnerPtr_(
nullptr)
607 ParentType(convertToTri(triangles, 0), points),
609 sortedEdgeFacesPtr_(
nullptr),
610 edgeOwnerPtr_(
nullptr)
620 sortedEdgeFacesPtr_(nullptr),
621 edgeOwnerPtr_(nullptr)
635 sortedEdgeFacesPtr_(nullptr),
636 edgeOwnerPtr_(nullptr)
648 sortedEdgeFacesPtr_(nullptr),
649 edgeOwnerPtr_(nullptr)
666 patches_(ts.patches()),
667 sortedEdgeFacesPtr_(nullptr),
668 edgeOwnerPtr_(nullptr)
675 patches_(move(ts.patches())),
676 sortedEdgeFacesPtr_(nullptr),
677 edgeOwnerPtr_(nullptr)
716 if (!sortedEdgeFacesPtr_)
718 calcSortedEdgeFaces();
721 return *sortedEdgeFacesPtr_;
732 return *edgeOwnerPtr_;
752 if (scaleFactor > 0 && scaleFactor != 1.0)
776 if (f[fp] < 0 || f[fp] > maxPointi)
780 <<
" uses point indices outside point range 0.." 794 bool hasInvalid =
false;
800 if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
803 valid[facei] =
false;
809 <<
"triangle " << facei
810 <<
" does not have three unique vertices:\n";
828 label neighbour = eFaces[i];
830 if (neighbour > facei)
837 ((f[0] == n[0]) || (f[0] == n[1]) || (f[0] == n[2]))
838 && ((f[1] == n[0]) || (f[1] == n[1]) || (f[1] == n[2]))
839 && ((f[2] == n[0]) || (f[2] == n[1]) || (f[2] == n[2]))
842 valid[facei] =
false;
848 <<
"triangles share the same vertices:\n" 849 <<
" face 1 :" << facei <<
endl;
855 << neighbour <<
endl;
876 (*this)[newFacei++] =
f;
883 <<
"Removing " <<
size() - newFacei
884 <<
" illegal faces." <<
endl;
886 (*this).setSize(newFacei);
900 const labelList& myFaces = eFaces[edgeI];
905 <<
"Edge " << edgeI <<
" with vertices " <<
edges()[edgeI]
906 <<
" has no edgeFaces" 909 else if (myFaces.
size() > 2 && verbose)
912 <<
"Edge " << edgeI <<
" with vertices " <<
edges()[edgeI]
913 <<
" has more than 2 faces connected to it : " << myFaces
923 stitchTriangles(small, verbose);
938 const label currentZone,
952 label facei = changedFaces[i];
958 label edgeI = fEdges[i];
960 if (!borderEdge[edgeI])
966 label nbrFacei = eFaces[j];
968 if (faceZone[nbrFacei] == -1)
970 faceZone[nbrFacei] = currentZone;
971 newChangedFaces.
append(nbrFacei);
973 else if (faceZone[nbrFacei] != currentZone)
976 <<
"Zones " << faceZone[nbrFacei]
977 <<
" at face " << nbrFacei
978 <<
" connects to zone " << currentZone
979 <<
" at face " << facei
987 if (newChangedFaces.empty())
992 changedFaces.
transfer(newChangedFaces);
1009 <<
"borderEdge boolList not same size as number of edges" <<
endl 1010 <<
"borderEdge:" << borderEdge.
size() <<
endl 1011 <<
"nEdges :" <<
nEdges()
1017 for (
label startFacei = 0;; zoneI++)
1020 for (; startFacei <
size(); startFacei++)
1022 if (faceZone[startFacei] == -1)
1028 if (startFacei >=
size())
1033 faceZone[startFacei] = zoneI;
1035 markZone(borderEdge, startFacei, zoneI, faceZone);
1060 forAll(include, oldFacei)
1062 if (include[oldFacei])
1065 faceMap[facei++] = oldFacei;
1073 if (!pointHad[labI])
1075 pointHad[labI] =
true;
1076 pointMap[pointi++] = labI;
1107 newPoints[pointi] = locPoints[pointMap[pointi]];
1108 oldToNew[pointMap[pointi]] = pointi;
1117 const labelledTri& tri = locFaces[faceMap[facei]];
1119 newTriangles[facei][0] = oldToNew[tri[0]];
1120 newTriangles[facei][1] = oldToNew[tri[1]];
1121 newTriangles[facei][2] = oldToNew[tri[2]];
1122 newTriangles[facei].region() = tri.region();
1136 faces[facei] =
operator[](facei).triFaceFace();
1155 this->weightedPointNormals()
1160 this->pointCoordSys(pointNormals)
1180 const edge&
e = fEdges[feI];
1182 edgeVectors[feI] = e.
vec(points);
1183 normalDifferences[feI] =
1184 pointNormals[meshPointMap[e[0]]]
1185 - pointNormals[meshPointMap[e[1]]];
1189 const vector& e0 = edgeVectors[0];
1191 const vector e1 = (e0 ^ eN);
1193 if (
magSqr(eN) < rootVSmall)
1198 triad faceCoordSys(e0, e1, eN);
1206 for (
label i = 0; i < 3; ++i)
1208 const scalar
x = edgeVectors[i] & faceCoordSys[0];
1209 const scalar
y = edgeVectors[i] & faceCoordSys[1];
1217 const scalar dndx = normalDifferences[i] & faceCoordSys[0];
1218 const scalar dndy = normalDifferences[i] & faceCoordSys[1];
1221 Z[1] += dndx*y + dndy*
x;
1228 const symmTensor2D secondFundamentalTensor(Z[0], Z[1], Z[2]);
1234 const label patchPointIndex = meshPointMap[f[fpi]];
1236 const triad& ptCoordSys = pointCoordSys[patchPointIndex];
1238 if (!ptCoordSys.
set())
1249 const triad rotatedFaceCoordSys = rotTensor &
tensor(faceCoordSys);
1255 ptCoordSys[0] & rotatedFaceCoordSys[0],
1256 ptCoordSys[0] & rotatedFaceCoordSys[1]
1260 ptCoordSys[1] & rotatedFaceCoordSys[0],
1261 ptCoordSys[1] & rotatedFaceCoordSys[1]
1272 projTensor.
x() & (secondFundamentalTensor & projTensor.
x()),
1273 projTensor.
x() & (secondFundamentalTensor & projTensor.
y()),
1274 projTensor.
y() & (secondFundamentalTensor & projTensor.
y())
1278 const scalar weight = pointNormalWeight
1281 meshPoints[patchPointIndex],
1287 pointFundamentalTensors[patchPointIndex] +=
1288 weight*projectedFundamentalTensor;
1290 accumulatedWeights[patchPointIndex] += weight;
1294 forAll(curvaturePointField, pi)
1296 pointFundamentalTensors[
pi] /= (accumulatedWeights[
pi] + small);
1298 const vector2D principalCurvatures =
1305 mag(principalCurvatures[0]),
1306 mag(principalCurvatures[1])
1313 return tcurvaturePointField;
1317 void Foam::triSurface::write
1320 const bool sortByRegion
1323 write(name, name.
ext(), sortByRegion);
1336 os.
check(
"triSurface::write(Ostream&)");
1340 void Foam::triSurface::write(
const Time& d)
const 1367 label pointi = f[fp];
1368 if (pointIsUsed.
set(pointi, 1))
1377 os <<
"Triangles : " <<
size() <<
endl 1378 <<
"Vertices : " << nPoints <<
endl 1379 <<
"Bounding Box : " << bb <<
endl;
1399 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)
fileName path() const
Return path.
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.
void normalize()
Normalize each set axis vector to have a unit magnitude.
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.
const fileName & caseName() const
Return case name.
void checkTriangles(const bool verbose)
Check/remove duplicate/degenerate triangles.
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.
triSurface()
Construct null.
std::remove_reference< ::Foam::List< labelledTri > >::type::value_type FaceType
scalar timeOutputValue() const
Return current time value.
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.
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.
ISstream & getLine(string &)
Read line into a string.