83 scalar minDistSqr = GREAT;
86 scalar almostMinDistSqr = GREAT;
87 label almostMinI = -1;
89 for (
const label pointi : meshPoints)
93 if (distSqr < minDistSqr)
95 almostMinDistSqr = minDistSqr;
101 else if (distSqr < almostMinDistSqr)
103 almostMinDistSqr = distSqr;
110 Info<<
"Found to point " << nearPoint <<
nl
111 <<
" nearest point : " << minI
114 <<
" next nearest point : " << almostMinI
115 <<
" distance " <<
Foam::sqrt(almostMinDistSqr)
118 if (almostMinDistSqr < 4*minDistSqr)
120 Info<<
"Next nearest too close to nearest. Aborting" <<
endl;
136 const point& nearPoint
145 scalar minDist = GREAT;
148 scalar almostMinDist = GREAT;
149 label almostMinI = -1;
151 for (
const edge&
e : edges)
153 pointHit pHit(
e.line(localPoints).nearestDist(nearPoint));
157 if (pHit.distance() < minDist)
159 almostMinDist = minDist;
162 minDist = pHit.distance();
163 minI = meshTools::findEdge
170 else if (pHit.distance() < almostMinDist)
172 almostMinDist = pHit.distance();
173 almostMinI = meshTools::findEdge
185 Info<<
"Did not find edge close to point " << nearPoint <<
endl;
192 Info<<
"Found to point " << nearPoint <<
nl
193 <<
" nearest edge : " << minI
194 <<
" distance " << minDist <<
" endpoints "
196 <<
" next nearest edge : " << almostMinI
197 <<
" distance " << almostMinDist <<
" endpoints "
201 if (almostMinDist < 2*minDist)
203 Info<<
"Next nearest too close to nearest. Aborting" <<
endl;
219 const point& nearPoint
225 scalar minDist = GREAT;
228 scalar almostMinDist = GREAT;
229 label almostMinI = -1;
237 if (pHit.distance() < minDist)
239 almostMinDist = minDist;
242 minDist = pHit.distance();
243 minI = patchFacei +
mesh.nInternalFaces();
245 else if (pHit.distance() < almostMinDist)
247 almostMinDist = pHit.distance();
248 almostMinI = patchFacei +
mesh.nInternalFaces();
255 Info<<
"Did not find face close to point " << nearPoint <<
endl;
262 Info<<
"Found to point " << nearPoint <<
nl
263 <<
" nearest face : " << minI
264 <<
" distance " << minDist
265 <<
" to face centre " <<
mesh.faceCentres()[minI] <<
nl
266 <<
" next nearest face : " << almostMinI
267 <<
" distance " << almostMinDist
268 <<
" to face centre " <<
mesh.faceCentres()[almostMinI] <<
nl
271 if (almostMinDist < 2*minDist)
273 Info<<
"Next nearest too close to nearest. Aborting" <<
endl;
287 label celli =
mesh.findCell(nearPoint);
291 scalar distToCcSqr =
magSqr(nearPoint -
mesh.cellCentres()[celli]);
296 scalar minDistSqr = GREAT;
298 for (
const label pointi : cPoints)
300 scalar distSqr =
magSqr(nearPoint -
mesh.points()[pointi]);
302 if (distSqr < minDistSqr)
304 minDistSqr = distSqr;
310 Info<<
"Found to point " << nearPoint <<
nl
311 <<
" nearest cell : " << celli
313 <<
" to cell centre " <<
mesh.cellCentres()[celli] <<
nl
314 <<
" nearest mesh point : " << minI
316 <<
" to " <<
mesh.points()[minI] <<
nl
319 if (minDistSqr < 4*distToCcSqr)
321 Info<<
"Mesh point too close to nearest cell centre. Aborting"
334int main(
int argc,
char *argv[])
338 "Manipulate mesh elements.\n"
339 "For example, moving points, splitting/collapsing edges etc."
342 argList::addOption(
"dict",
"file",
"Alternative modifyMeshDict");
344 argList::noFunctionObjects();
350 const word oldInstance =
mesh.pointsInstance();
352 const bool overwrite =
args.
found(
"overwrite");
354 Info<<
"Reading modifyMeshDict\n" <<
endl;
366 dict.lookup(
"facesToTriangulate")
370 dict.readIfPresent(
"facesToSplit", facesToSplit);
375 || edgesToSplit.
size()
376 || facesToTriangulate.
size()
377 || facesToSplit.
size()
386 bool cellsToSplit = cellsToPyramidise.
size();
392 <<
" Boundary cutting module:" <<
nl
393 <<
" points to move :" << pointsToMove.
size() <<
nl
394 <<
" edges to split :" << edgesToSplit.
size() <<
nl
395 <<
" faces to split :" << facesToSplit.
size() <<
nl
396 <<
" faces to triangulate:" << facesToTriangulate.
size() <<
nl
397 <<
" Cell splitting module:" <<
nl
398 <<
" cells to split :" << cellsToPyramidise.
size() <<
nl
399 <<
" Edge collapsing module:" <<
nl
400 <<
" edges to collapse :" << edgesToCollapse.
size() <<
nl
407 || (cutBoundary && cellsToSplit)
412 <<
"Used more than one mesh modifying module "
413 <<
"(boundary cutting, cell splitting, edge collapsing)" <<
nl
414 <<
"Please do them in separate passes." <<
exit(FatalError);
423 mesh.nBoundaryFaces(),
424 mesh.nInternalFaces()
432 bool validInputs =
true;
435 Info<<
nl <<
"Looking up points to move ..." <<
nl <<
endl;
439 const label pointi = findPoint(allBoundary, pts.first());
441 if (pointi == -1 || !pointToPos.insert(pointi, pts.second()))
443 Info<<
"Could not insert mesh point " << pointi
444 <<
" for input point " << pts.first() <<
nl
445 <<
"Perhaps the point is already marked for moving?" <<
endl;
451 Info<<
nl <<
"Looking up edges to split ..." <<
nl <<
endl;
460 || !edgeToCuts.insert(edgeI,
List<point>(1, pts.second()))
463 Info<<
"Could not insert mesh edge " << edgeI
464 <<
" for input point " << pts.first() <<
nl
465 <<
"Perhaps the edge is already marked for cutting?" <<
endl;
472 Info<<
nl <<
"Looking up faces to triangulate ..." <<
nl <<
endl;
476 label facei = findFace(
mesh, allBoundary, pts.first());
478 if (facei == -1 || !faceToDecompose.insert(facei, pts.second()))
480 Info<<
"Could not insert mesh face " << facei
481 <<
" for input point " << pts.first() <<
nl
482 <<
"Perhaps the face is already marked for splitting?" <<
endl;
489 Info<<
nl <<
"Looking up faces to split ..." <<
nl <<
endl;
493 label facei = findFace(
mesh, allBoundary, pts.first());
496 Info<<
"Could not insert mesh face " << facei
497 <<
" for input point " << pts.first() <<
nl
498 <<
"Perhaps the face is already marked for splitting?" <<
endl;
511 const label
p0 = findPoint(pp, pts.first());
512 const label p1 = findPoint(pp, pts.second());
520 && faceToSplit.insert(facei,
labelPair(
f.find(
p0),
f.find(p1)))
525 Info<<
"Could not insert mesh face " << facei
526 <<
" for input coordinates " << pts
527 <<
" with vertices " <<
p0 <<
" and " << p1 <<
nl
528 <<
"Perhaps the face is already marked for splitting?"
536 Info<<
nl <<
"Looking up cells to convert to pyramids around"
537 <<
" cell centre ..." <<
nl <<
endl;
541 label celli = findCell(
mesh, pts.first());
543 if (celli == -1 || !cellToPyrCentre.insert(celli, pts.second()))
545 Info<<
"Could not insert mesh cell " << celli
546 <<
" for input point " << pts.first() <<
nl
547 <<
"Perhaps the cell is already marked for splitting?" <<
endl;
554 Info<<
nl <<
"Looking up edges to collapse ..." <<
nl <<
endl;
560 if (edgeI == -1 || !edgeToPos.insert(edgeI, pts.second()))
562 Info<<
"Could not insert mesh edge " << edgeI
563 <<
" for input point " << pts.first() <<
nl
564 <<
"Perhaps the edge is already marked for collaping?" <<
endl;
574 Info<<
nl <<
"There was a problem in one of the inputs in the"
575 <<
" dictionary. Not modifying mesh." <<
endl;
577 else if (cellToPyrCentre.size())
579 Info<<
nl <<
"All input cells located. Modifying mesh." <<
endl;
588 cutter.setRefinement(cellToPyrCentre, meshMod);
593 if (morphMap().hasMotionPoints())
595 mesh.movePoints(morphMap().preMotionPoints());
598 cutter.updateMesh(morphMap());
606 mesh.setInstance(oldInstance);
612 topoSet::removeFiles(
mesh);
613 processorMeshes::removeFiles(
mesh);
615 else if (edgeToPos.size())
617 Info<<
nl <<
"All input edges located. Modifying mesh." <<
endl;
633 label edgeI = iter.key();
634 const edge&
e = edges[edgeI];
637 collapsePointToLocation.set(
e[1],
points[
e[0]]);
639 newPoints[
e[0]] = iter.val();
643 mesh.movePoints(newPoints);
649 cutter.consistentCollapse
653 collapsePointToLocation,
662 cutter.setRefinement(allPointInfo, meshMod);
667 if (morphMap().hasMotionPoints())
669 mesh.movePoints(morphMap().preMotionPoints());
682 mesh.setInstance(oldInstance);
688 topoSet::removeFiles(
mesh);
689 processorMeshes::removeFiles(
mesh);
693 Info<<
nl <<
"All input points located. Modifying mesh." <<
endl;
714 if (morphMap().hasMotionPoints())
716 mesh.movePoints(morphMap().preMotionPoints());
719 cutter.updateMesh(morphMap());
727 mesh.setInstance(oldInstance);
733 topoSet::removeFiles(
mesh);
734 processorMeshes::removeFiles(
mesh);
static constexpr label size() noexcept
Return the number of elements in the FixedList.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
A HashTable to objects of type <T> with a label key.
An ordered pair of two objects of type <T> with first() and second() elements.
Describes the interaction of a face and a point. It carries the info of a successful hit and (if succ...
A list of faces which address into the list of points.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
A List obtained as a section of another List.
void size(const label n)
Older name for setAddressableSize.
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...
Does modifications to boundary faces.
Does pyramidal decomposition of selected cells. So all faces will become base of pyramid with as top ...
Does polyTopoChanges to remove edges. Can remove faces due to edge collapse but can not remove cells ...
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
A face is a list of labels corresponding to mesh vertices.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Calculates points shared by more than two processor patches or cyclic patches.
Direct mesh changes based on v1.3 polyTopoChange syntax.
Cell-face mesh analysis engine.
A class for handling words, derived from Foam::string.
label collapseEdge(triSurface &surf, const scalar minLen)
Keep collapsing all edges < minLen.
const volScalarField & p0
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const word dictName("faMeshDefinition")
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar sqrt(const dimensionedScalar &ds)
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.