81 packedLst[packedI++] = i;
84 packedLst.setSize(packedI);
99 packedElems[packedI++] = elems[i];
113 const scalar cosAngle,
114 const scalar sinAngle,
126 bool sameFaceOrder = !((own[f0] == celli) ^ (own[f1] == celli));
130 scalar normalCosAngle = n0 & n1;
134 normalCosAngle = -normalCosAngle;
144 mesh.faceCentres()[f1] -
mesh.faceCentres()[f0]
147 scalar fcCosAngle = n0 & c1c0;
149 if (own[f0] != celli)
151 fcCosAngle = -fcCosAngle;
164 if (normalCosAngle < cosAngle)
185 if (normalCosAngle > cosAngle)
216 const edge&
e = edges[edgeI];
226 const cell& cFaces =
mesh.cells()[celli];
228 for (
const label facei : cFaces)
230 const face&
f = faces[facei];
232 label fp0 =
f.find(
e[0]);
233 label fp1 =
f.find(
e[1]);
269 if (leftI == -1 || rightI == -1)
272 <<
" rightI:" << rightI <<
abort(FatalError);
277 const face& leftF = faces[leftI];
279 label leftV = leftF[(leftFp + 2) % leftF.
size()];
281 const face& rightF = faces[rightI];
283 label rightV = rightF[(rightFp + 2) % rightF.
size()];
286 loop[0] = ev.vertToEVert(
e[0]);
287 loop[1] = ev.vertToEVert(leftV);
288 loop[2] = ev.vertToEVert(rightV);
289 loop[3] = ev.vertToEVert(
e[1]);
292 loopWeights[0] = -GREAT;
293 loopWeights[1] = -GREAT;
294 loopWeights[2] = -GREAT;
295 loopWeights[3] = -GREAT;
299 cellEdgeWeights.
append(loopWeights);
328 vector planeN = eVec ^ halfNorm;
333 planeN += 0.01*halfNorm;
337 plane cutPlane(
mesh.points()[
e.start()], planeN);
359 cellEdgeWeights.
append(loopWeights);
401 for (
const label celli : cellsToCut)
403 const labelList& cEdges = cellEdges[celli];
405 for (
const label edgeI : cEdges)
408 meshTools::getEdgeFaces(
mesh, celli, edgeI, f0, f1);
429 bool splitOk =
false;
449 if ((own[f0] == celli) ^ (own[f1] == celli))
452 halfNorm = 0.5*(n0 - n1);
458 halfNorm = 0.5*(n0 + n1);
484 label index = cellLoops.
size() - 1;
485 const labelList& loop = cellLoops[index];
486 const scalarField& loopWeights = cellEdgeWeights[index];
494 edgeIsCut[ev.getEdge(
cut)] =
true;
495 edgeWeight[ev.getEdge(
cut)] = loopWeights[i];
499 vertIsCut[ev.getVertex(
cut)] =
true;
517int main(
int argc,
char *argv[])
521 "Split cells with flat faces"
524 argList::noParallel();
535 "Split cells from specified cellSet only"
537 argList::addBoolOption
540 "Use geometric cut for hexes as well"
546 "Edge snap tolerance (default 0.2)"
549 argList::noFunctionObjects();
555 const word oldInstance =
mesh.pointsInstance();
557 const scalar featureAngle =
args.
get<scalar>(1);
562 const bool geometry =
args.
found(
"geometry");
563 const bool overwrite =
args.
found(
"overwrite");
567 Info<<
"Trying to split cells with internal angles > feature angle\n" <<
nl
568 <<
"featureAngle : " << featureAngle <<
nl
569 <<
"edge snapping tol : " << edgeTol <<
nl;
572 Info<<
"candidate cells : cellSet " <<
args[
"set"] <<
nl;
576 Info<<
"candidate cells : all cells" <<
nl;
580 Info<<
"hex cuts : geometric; using edge tolerance" <<
nl;
584 Info<<
"hex cuts : topological; cut to opposite edge" <<
nl;
592 geomCellLooper::setSnapTol(edgeTol);
610 for (label celli = 0; celli <
mesh.nCells(); celli++)
612 cellsToCut.insert(celli);
637 cutSet.insert(cutCells);
640 Info<<
"Writing " << cutSet.size() <<
" cells to cut to cellSet "
641 << cutSet.instance()/cutSet.local()/cutSet.name()
654 Info<<
"Actually cut cells:" << cuts.nLoops() <<
nl <<
endl;
656 if (cuts.nLoops() == 0)
662 forAll(cuts.cellLoops(), celli)
664 if (cuts.cellLoops()[celli].size())
668 cellsToCut.erase(celli);
679 cutter.setRefinement(cuts, meshMod);
691 if (morphMap().hasMotionPoints())
693 mesh.movePoints(morphMap().preMotionPoints());
697 cutter.updateMesh(morphMap());
700 cellsToCut.updateMesh(morphMap());
702 Info<<
"Remaining:" << cellsToCut.size() <<
endl;
707 mesh.setInstance(oldInstance);
710 Info<<
"Writing refined morphMesh to time " <<
runTime.timeName()
Various functions to operate on Lists.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
void append(const T &val)
Copy append an element to the end of this list.
void setSize(const label n)
Alias for resize()
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
void size(const label n)
Older name for setAddressableSize.
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.
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Description of cuts across cells.
Maps a geometry to a set of cell primitives.
A collection of cell labels.
A cell is defined as a list of faces with extra functionality.
Combines edge or vertex in single label. Used to specify cuts across cell circumference.
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.
Implementation of cellLooper. Does pure geometric cut through cell.
virtual bool cut(const vector &refDir, const label celli, const boolList &vertIsCut, const boolList &edgeIsCut, const scalarField &edgeWeight, labelList &loop, scalarField &loopWeights) const
Create cut along circumference of celli. Gets current mesh cuts.
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
Mesh consisting of general polyhedral cells.
Direct mesh changes based on v1.3 polyTopoChange syntax.
Cell-face mesh analysis engine.
Patchify triangles based on orientation w.r.t other (triangulated or triangulatable) surfaces.
Description of cell after splitting. Contains cellLabel and pointers to cells it it split in....
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.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
dimensionedScalar sin(const dimensionedScalar &ds)
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.
errorManip< error > abort(error &err)
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.