89void simpleMarkFeatures
92 const bitSet& isBoundaryEdge,
93 const scalar featureAngle,
94 const bool concaveMultiCells,
95 const bool doNotPreserveFaceZones,
124 label meshEdgeI = meshEdges[edgeI];
125 featureEdgeSet.
insert(meshEdgeI);
126 singleCellFeaturePointSet.
insert(
mesh.edges()[meshEdgeI][0]);
127 singleCellFeaturePointSet.
insert(
mesh.edges()[meshEdgeI][1]);
145 mesh.nBoundaryFaces(),
146 mesh.nInternalFaces()
152 allBoundary.checkPointManifold(
false, &singleCellFeaturePointSet);
156 const labelList& meshPoints = allBoundary.meshPoints();
160 const labelList& eFaces = edgeFaces[edgeI];
162 if (eFaces.
size() > 2)
164 const edge&
e = allBoundary.edges()[edgeI];
171 singleCellFeaturePointSet.
insert(meshPoints[
e[0]]);
172 singleCellFeaturePointSet.
insert(meshPoints[
e[1]]);
179 const labelList& eFaces = edgeFaces[edgeI];
181 if (eFaces.
size() == 2)
183 label f0 = eFaces[0];
184 label f1 = eFaces[1];
187 const vector& n0 = allBoundary.faceNormals()[f0];
188 const vector& n1 = allBoundary.faceNormals()[f1];
190 if ((n0 & n1) < minCos)
192 const edge&
e = allBoundary.edges()[edgeI];
193 label v0 = meshPoints[
e[0]];
194 label v1 = meshPoints[
e[1]];
196 label meshEdgeI = meshTools::findEdge(
mesh, v0, v1);
197 featureEdgeSet.
insert(meshEdgeI);
203 allBoundary[f1].centre(allBoundary.points())
204 - allBoundary[f0].centre(allBoundary.points())
207 if (concaveMultiCells && (c1c0 & n0) > SMALL)
210 Info<<
"Detected concave feature edge:" << edgeI
211 <<
" cos:" << (c1c0 & n0)
213 << allBoundary.points()[v0]
214 << allBoundary.points()[v1]
217 singleCellFeaturePointSet.
erase(v0);
218 multiCellFeaturePointSet.
insert(v0);
219 singleCellFeaturePointSet.
erase(v1);
220 multiCellFeaturePointSet.
insert(v1);
225 if (!multiCellFeaturePointSet.
found(v0))
227 singleCellFeaturePointSet.
insert(v0);
229 if (!multiCellFeaturePointSet.
found(v1))
231 singleCellFeaturePointSet.
insert(v1);
245 for (label facei =
mesh.nInternalFaces(); facei <
mesh.nFaces(); facei++)
247 featureFaceSet.insert(facei);
253 if (doNotPreserveFaceZones)
255 if (faceZones.
size() > 0)
258 <<
"Detected " << faceZones.
size()
259 <<
" faceZones. These will not be preserved."
265 if (faceZones.
size() > 0)
267 Info<<
"Detected " << faceZones.
size()
268 <<
" faceZones. Preserving these by marking their"
269 <<
" points, edges and faces as features." <<
endl;
274 const faceZone& fz = faceZones[zoneI];
276 Info<<
"Inserting all faces in faceZone " << fz.
name()
277 <<
" as features." <<
endl;
285 featureFaceSet.insert(facei);
290 singleCellFeaturePointSet.
erase(
f[fp]);
291 multiCellFeaturePointSet.
insert(
f[fp]);
294 featureEdgeSet.
insert(fEdges[fp]);
301 featureFaces = featureFaceSet.toc();
302 featureEdges = featureEdgeSet.
toc();
303 singleCellFeaturePoints = singleCellFeaturePointSet.
toc();
304 multiCellFeaturePoints = multiCellFeaturePointSet.
toc();
314 const labelList& singleCellFeaturePoints,
320 Info<<
"Dumping centres of featureFaces to obj file " << str.name()
324 meshTools::writeOBJ(str,
mesh.faceCentres()[featureFaces[i]]);
329 Info<<
"Dumping featureEdges to obj file " << str.name() <<
endl;
334 const edge&
e =
mesh.edges()[featureEdges[i]];
335 meshTools::writeOBJ(str,
mesh.points()[
e[0]]);
337 meshTools::writeOBJ(str,
mesh.points()[
e[1]]);
339 str<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
343 OFstream str(
"singleCellFeaturePoints.obj");
344 Info<<
"Dumping featurePoints that become a single cell to obj file "
345 << str.name() <<
endl;
346 forAll(singleCellFeaturePoints, i)
348 meshTools::writeOBJ(str,
mesh.points()[singleCellFeaturePoints[i]]);
352 OFstream str(
"multiCellFeaturePoints.obj");
353 Info<<
"Dumping featurePoints that become multiple cells to obj file "
354 << str.name() <<
endl;
355 forAll(multiCellFeaturePoints, i)
357 meshTools::writeOBJ(str,
mesh.points()[multiCellFeaturePoints[i]]);
363int main(
int argc,
char *argv[])
367 "Creates the dual of a polyMesh,"
368 " adhering to all the feature and patch edges."
372 argList::noParallel();
380 argList::addBoolOption
383 "Have multiple faces in between cells"
385 argList::addBoolOption
388 "Split cells on concave boundary edges into multiple cells"
390 argList::addBoolOption
392 "doNotPreserveFaceZones",
393 "Disable the default behaviour of preserving faceZones by having"
394 " multiple faces in between cells"
401 const word oldInstance =
mesh.pointsInstance();
407 for (label facei =
mesh.nInternalFaces(); facei <
mesh.nFaces(); facei++)
413 isBoundaryEdge.
set(fEdges[i]);
417 const scalar featureAngle =
args.
get<scalar>(1);
420 Info<<
"Feature:" << featureAngle <<
endl
421 <<
"minCos :" << minCos <<
endl
425 const bool splitAllFaces =
args.
found(
"splitAllFaces");
428 Info<<
"Splitting all internal faces to create multiple faces"
429 <<
" between two cells." <<
nl
433 const bool overwrite =
args.
found(
"overwrite");
434 const bool doNotPreserveFaceZones =
args.
found(
"doNotPreserveFaceZones");
435 const bool concaveMultiCells =
args.
found(
"concaveMultiCells");
436 if (concaveMultiCells)
438 Info<<
"Generating multiple cells for points on concave feature edges."
459 doNotPreserveFaceZones,
463 singleCellFeaturePoints,
464 multiCellFeaturePoints
483 singleCellFeaturePoints,
484 multiCellFeaturePoints
533 dualMaker.setRefinement
538 singleCellFeaturePoints,
539 multiCellFeaturePoints,
547 mesh.updateMesh(map());
550 if (map().hasMotionPoints())
552 mesh.movePoints(map().preMotionPoints());
561 mesh.setInstance(oldInstance);
567 topoSet::removeFiles(
mesh);
568 processorMeshes::removeFiles(
mesh);
Field reading functions for post-processing utilities.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
List< Key > toc() const
The table of contents (the keys) in unsorted order.
bool found(const Key &key) const
Return true if hashed entry is found in table.
bool erase(const iterator &iter)
Erase an entry specified by given iterator.
List of IOobjects with searching and retrieving facilities.
Output to file stream, using an OSstream.
A list of faces which address into the list of points.
label nEdges() const
Number of edges in patch.
label nInternalEdges() const
Number of internal edges.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
A List obtained as a section of another List.
void size(const label n)
Older name for setAddressableSize.
label size() const noexcept
The number of elements in the list.
T get(const label index) const
Get a value from the argument at index.
bool found(const word &optName) const
Return true if the named option is found.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
void set(const bitSet &bitset)
Set specified bits from another bitset.
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.
Creates dual of polyMesh. Every point becomes a cell (or multiple cells for feature points),...
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
A patch is a list of labels that address the faces in the global face list.
const labelList & meshEdges() const
Return global edge index for local edges.
Direct mesh changes based on v1.3 polyTopoChange syntax.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
const word & name() const noexcept
The zone name.
const polyBoundaryMesh & patches
#define WarningInFunction
Report a warning using Foam::Warning.
wordList ReadFields(const typename GeoMesh::Mesh &mesh, const IOobjectList &objects, PtrList< GeometricField< Type, PatchField, GeoMesh > > &fields, const bool syncPar=true, const bool readOldTime=false)
Read Geometric fields of templated type.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
messageStream Info
Information stream (stdout output on master, null elsewhere)
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Unit conversion functions.