60 for (i=ts.
size()-1; i>=0; i--)
73 for (
label j=i; j>=0; j--)
75 if (
isFile(d.
path()/ts[j].name()/typeName/foamName))
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
120 labelledTri& tri = triFaces[facei];
125 tri.region() = defaultRegion;
135 const label defaultRegion
138 List<labelledTri> triFaces(faces.size());
144 labelledTri& tri = triFaces[facei];
149 tri.region() = defaultRegion;
158 void Foam::triSurface::printTriangle
162 const labelledTri&
f,
167 << pre.c_str() <<
"vertex numbers:"
168 <<
f[0] <<
' ' <<
f[1] <<
' ' <<
f[2] <<
endl
169 << pre.c_str() <<
"vertex coords :"
171 << pre.c_str() <<
"region :" <<
f.region() <<
endl
176 Foam::string Foam::triSurface::getLineNoComment(IFstream& is)
183 while ((line.empty() || line[0] ==
'#') && is.good());
189 Foam::scalar Foam::triSurface::pointNormalWeight
217 tmp<vectorField> tpointNormals
225 const labelList& meshPoints = this->meshPoints();
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;
256 triadField& pointCoordSys = tpointCoordSys.ref();
259 const labelList& meshPoints = this->meshPoints();
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
298 const fileName&
name,
311 fileName unzipName =
name.lessExt();
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
371 const fileName&
name,
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")
412 writeAC(OFstream(
name)());
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'"
434 SortableList<label> sortedRegion(size());
436 forAll(sortedRegion, facei)
438 sortedRegion[facei] = operator[](facei).region();
442 faceMap = sortedRegion.indices();
445 label maxRegion = patches_.size()-1;
452 operator[](
faceMap.last()).region()
462 label region = operator[](facei).region();
464 newPatches[region].size()++;
470 label startFacei = 0;
471 forAll(newPatches, newPatchi)
473 surfacePatch& newPatch = 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();
490 newPatch.name() = word(
"patch") +
name(newPatchi);
495 (oldPatchi < patches_.size())
496 && (patches_[oldPatchi].geometricType() !=
"")
499 newPatch.geometricType() = patches_[oldPatchi].geometricType();
503 newPatch.geometricType() =
"empty";
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)
552 sortedEdgeFacesPtr_(nullptr),
553 edgeOwnerPtr_(nullptr)
567 sortedEdgeFacesPtr_(nullptr),
568 edgeOwnerPtr_(nullptr)
581 sortedEdgeFacesPtr_(nullptr),
582 edgeOwnerPtr_(nullptr)
594 sortedEdgeFacesPtr_(nullptr),
595 edgeOwnerPtr_(nullptr)
609 sortedEdgeFacesPtr_(nullptr),
610 edgeOwnerPtr_(nullptr)
620 sortedEdgeFacesPtr_(nullptr),
621 edgeOwnerPtr_(nullptr)
635 sortedEdgeFacesPtr_(nullptr),
636 edgeOwnerPtr_(nullptr)
648 sortedEdgeFacesPtr_(nullptr),
649 edgeOwnerPtr_(nullptr)
667 sortedEdgeFacesPtr_(nullptr),
668 edgeOwnerPtr_(nullptr)
676 sortedEdgeFacesPtr_(nullptr),
677 edgeOwnerPtr_(nullptr)
693 ParentType::clearTopology();
701 ParentType::clearPatchMeshAddr();
707 ParentType::clearOut();
710 clearPatchMeshAddr();
716 if (!sortedEdgeFacesPtr_)
718 calcSortedEdgeFaces();
721 return *sortedEdgeFacesPtr_;
732 return *edgeOwnerPtr_;
742 ParentType::clearGeom();
745 storedPoints() = newPoints;
752 if (scaleFactor > 0 && scaleFactor != 1.0)
758 ParentType::clearGeom();
760 storedPoints() *= scaleFactor;
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";
817 const labelList& fEdges = faceEdges()[facei];
824 const labelList& eFaces = edgeFaces()[fEdges[fp]];
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);
928 checkTriangles(verbose);
938 const label currentZone,
952 label facei = changedFaces[i];
954 const labelList& fEdges = faceEdges()[facei];
958 label edgeI = fEdges[i];
960 if (!borderEdge[edgeI])
962 const labelList& eFaces = edgeFaces()[edgeI];
966 label nbrFacei = eFaces[j];
971 newChangedFaces.
append(nbrFacei);
973 else if (
faceZone[nbrFacei] != currentZone)
977 <<
" at face " << nbrFacei
978 <<
" connects to zone " << currentZone
979 <<
" at face " << facei
987 if (newChangedFaces.
empty())
992 changedFaces.
transfer(newChangedFaces);
1006 if (borderEdge.
size() != nEdges())
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++)
1028 if (startFacei >= size())
1035 markZone(borderEdge, startFacei, zoneI,
faceZone);
1060 forAll(include, oldFacei)
1062 if (include[oldFacei])
1073 if (!pointHad[labI])
1075 pointHad[labI] =
true;
1076 pointMap[pointi++] = labI;
1099 subsetMeshMap(include, pointMap,
faceMap);
1107 newPoints[pointi] = locPoints[pointMap[pointi]];
1108 oldToNew[pointMap[pointi]] = pointi;
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();
1150 const labelList& meshPoints = this->meshPoints();
1151 const Map<label>& meshPointMap = this->meshPointMap();
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;
1296 pointFundamentalTensors[
pi] /= (accumulatedWeights[
pi] + small);
1298 const vector2D principalCurvatures =
1303 const scalar curvature =
max
1305 mag(principalCurvatures[0]),
1306 mag(principalCurvatures[1])
1310 curvaturePointField[meshPoints[
pi]] = curvature;
1313 return tcurvaturePointField;
1317 void Foam::triSurface::write
1320 const bool sortByRegion
1336 os.
check(
"triSurface::write(Ostream&)");
1340 void Foam::triSurface::write(
const Time& d)
const
1344 fileName foamPath(d.
path()/triSurfInstance(d)/typeName/foamFile);
1368 if (pointIsUsed.
set(pointi, 1))
1377 os <<
"Triangles : " << size() <<
endl
1379 <<
"Bounding Box : " << bb <<
endl;
1389 storedPoints() = ts.
points();
1398 storedPoints() = move(ts.points());
1399 patches_ = move(ts.patches());
Functions used by OpenFOAM that are specific to POSIX compliant operating systems and need to be repl...
#define forAll(list, i)
Loop across all elements in list.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
virtual bool check(const char *operation) const
Check IOstream status for given operation.
An Istream is an abstract base class for all input systems (streams, files, token lists etc)....
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
void size(const label)
Override size to be inconsistent with allocated storage.
void operator=(const UList< T > &)
Assignment to UList operator. Takes linear time.
void setSize(const label)
Reset size of List.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
void set(const PackedList< 1 > &)
Set specified bits.
const Field< PointType > & points() const
Return reference to global points.
std::remove_reference< ::Foam::List< labelledTri > >::type::value_type FaceType
Templated 2D symmetric tensor derived from VectorSpace adding construction from 4 components,...
A templated 2D square symmetric matrix of objects of <T>, where the n x n matrix dimension is known a...
Templated 2D tensor derived from VectorSpace adding construction from 4 components,...
Vector2D< Cmpt > y() const
Vector2D< Cmpt > x() const
static const word & constant()
Return constant name.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
instantList times() const
Search the case for valid time directories.
fileName path() const
Explicitly inherit path from TimePaths to disambiguate from.
scalar userTimeValue() const
Return current user time value.
const fileName & caseName() const
Explicitly inherit caseName from TimePaths to disambiguate from.
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
A bounding box defined in terms of the points at its extremities.
const point & min() const
Minimum point defining the bounding box.
const point & max() const
Maximum point defining the bounding box.
static const boundBox invertedBox
A very large inverted boundBox: min/max == +/- vGreat.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
A subset of mesh faces organised as a primitive patch.
A face is a list of labels corresponding to mesh vertices.
A class for handling file names.
Triangle with additional region number.
label region() const
Return region label.
A class for handling character strings derived from std::string.
A class for managing temporary objects.
T & ref() const
Return non-const reference or generate a fatal error.
A triangular face using a FixedList of labels corresponding to mesh vertices.
Triangulated surface description with patch information.
virtual void scalePoints(const scalar)
Scale points. A non-positive factor is ignored.
static fileName triSurfInstance(const Time &)
Name of triSurface directory to use.
triSurface()
Construct null.
triSurface subsetMesh(const boolList &include, labelList &pointMap, labelList &faceMap) const
Return new surface. Returns pointMap, faceMap from.
label markZones(const boolList &borderEdge, labelList &faceZone) const
(size and) fills faceZone with zone of face. Zone is area
tmp< scalarField > curvature() const
Return the curvature of surface at the points.
void operator=(const triSurface &)
void cleanup(const bool verbose)
Remove non-valid triangles.
const geometricSurfacePatchList & patches() const
void checkTriangles(const bool verbose)
Check/remove duplicate/degenerate triangles.
void checkEdges(const bool verbose)
Check triply (or more) connected edges.
void markZone(const boolList &borderEdge, const label facei, const label currentZone, labelList &faceZone) const
Fill faceZone with currentZone for every face reachable.
void writeStats(Ostream &) const
Write some statistics.
faceList faces() const
Return the list of triangles as a faceList.
const labelList & edgeOwner() const
If 2 face neighbours: label of face where ordering of edge.
virtual void movePoints(const pointField &)
Move points.
void subsetMeshMap(const boolList &include, labelList &pointMap, labelList &faceMap) const
'Create' sub mesh, including only faces for which
void clearPatchMeshAddr()
const labelListList & sortedEdgeFaces() const
Return edge-face addressing sorted (for edges with more than.
virtual ~triSurface()
Destructor.
Representation of a 3D Cartesian coordinate system as a Vector of vectors.
bool set(const direction d) const
Is the vector in the direction d set.
void normalise()
Normalise each set axis vector to have a unit magnitude.
A class for handling words, derived from string.
Template functions to aid in the implementation of demand driven data.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const fvPatchList & patches
#define WarningInFunction
Report a warning using Foam::Warning.
bool valid(const PtrList< ModelType > &l)
void write(std::ostream &os, const bool binary, List< floatScalar > &fField)
Write floats ascii or binary.
bool isFile(const fileName &, const bool checkVariants=true, const bool followLink=true)
Does the name exist as a file in the file system?
errorManipArg< error, int > exit(error &err, const int errNo=1)
bool read(const char *, int32_t &)
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from unit vector n1 to n2.
List< label > labelList
A List of labels.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Tensor< scalar > tensor
Tensor of scalars.
Ostream & endl(Ostream &os)
Add newline and flush stream.
word name(const bool)
Return a word representation of a bool.
errorManip< error > abort(error &err)
bool exists(const fileName &, const bool checkVariants=true, const bool followLink=true)
Does the name exist (as directory or file) in the file system?
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
void deleteDemandDrivenData(DataPtr &dataPtr)
Field< triad > triadField
Specialisation of Field<T> for triad.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
List< triFace > triFaceList
list of triFaces
Vector< scalar > vector
A scalar version of the templated Vector.
List< labelList > labelListList
A List of labelList.
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
dimensionedVector eigenValues(const dimensionedTensor &dt)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
Field< vector > vectorField
Specialisation of Field<T> for vector.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
void writeVTK(OFstream &os, const Type &value)
prefixOSstream Pout(cout, "Pout")
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
Ostream & operator<<(Ostream &, const ensightPart &)
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
Solve the matrix using LU decomposition with pivoting.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
List< surfacePatch > surfacePatchList
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]