splitCells.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 splitCells
29
30Group
31 grpMeshAdvancedUtilities
32
33Description
34 Utility to split cells with flat faces.
35
36 Uses a geometric cut with a plane dividing the edge angle into two so
37 might produce funny cells. For hexes it will use by default a cut from
38 edge onto opposite edge (i.e. purely topological).
39
40 Options:
41 - split cells from cellSet only
42 - use geometric cut for hexes as well
43
44 The angle is the angle between two faces sharing an edge as seen from
45 inside each cell. So a cube will have all angles 90. If you want
46 to split cells with cell centre outside use e.g. angle 200
47
48\*---------------------------------------------------------------------------*/
49
50#include "argList.H"
51#include "Time.H"
52#include "polyTopoChange.H"
53#include "polyTopoChanger.H"
54#include "mapPolyMesh.H"
55#include "polyMesh.H"
56#include "cellCuts.H"
57#include "cellSet.H"
58#include "cellModel.H"
59#include "meshCutter.H"
60#include "unitConversion.H"
61#include "geomCellLooper.H"
62#include "plane.H"
63#include "edgeVertex.H"
64#include "meshTools.H"
65#include "ListOps.H"
66
67using namespace Foam;
68
69// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
70
71
72labelList pack(const boolList& lst)
73{
74 labelList packedLst(lst.size());
75 label packedI = 0;
76
77 forAll(lst, i)
78 {
79 if (lst[i])
80 {
81 packedLst[packedI++] = i;
82 }
83 }
84 packedLst.setSize(packedI);
85
86 return packedLst;
87}
88
89
90scalarField pack(const boolList& lst, const scalarField& elems)
91{
92 scalarField packedElems(lst.size());
93 label packedI = 0;
94
95 forAll(lst, i)
96 {
97 if (lst[i])
98 {
99 packedElems[packedI++] = elems[i];
100 }
101 }
102 packedElems.setSize(packedI);
103
104 return packedElems;
105}
106
107
108// Given sin and cos of max angle between normals calculate whether f0 and f1
109// on celli make larger angle. Uses sinAngle only for quadrant detection.
110bool largerAngle
111(
112 const primitiveMesh& mesh,
113 const scalar cosAngle,
114 const scalar sinAngle,
115
116 const label celli,
117 const label f0, // face label
118 const label f1,
119
120 const vector& n0, // normal at f0
121 const vector& n1
122)
123{
124 const labelList& own = mesh.faceOwner();
125
126 bool sameFaceOrder = !((own[f0] == celli) ^ (own[f1] == celli));
127
128 // Get cos between faceArea vectors. Correct so flat angle (180 degrees)
129 // gives -1.
130 scalar normalCosAngle = n0 & n1;
131
132 if (sameFaceOrder)
133 {
134 normalCosAngle = -normalCosAngle;
135 }
136
137
138 // Get cos between faceCentre and normal vector to determine in
139 // which quadrant angle is. (Is correct for unwarped faces only!)
140 // Correct for non-outwards pointing normal.
141 const vector c1c0 =
143 (
144 mesh.faceCentres()[f1] - mesh.faceCentres()[f0]
145 );
146
147 scalar fcCosAngle = n0 & c1c0;
148
149 if (own[f0] != celli)
150 {
151 fcCosAngle = -fcCosAngle;
152 }
153
154 if (sinAngle < 0.0)
155 {
156 // Looking for concave angles (quadrant 3 or 4)
157 if (fcCosAngle <= 0)
158 {
159 // Angle is convex so smaller.
160 return false;
161 }
162 else
163 {
164 if (normalCosAngle < cosAngle)
165 {
166 return false;
167 }
168 else
169 {
170 return true;
171 }
172 }
173 }
174 else
175 {
176 // Looking for convex angles (quadrant 1 or 2)
177 if (fcCosAngle > 0)
178 {
179 // Concave angle
180 return true;
181 }
182 else
183 {
184 // Convex. Check cos of normal vectors.
185 if (normalCosAngle > cosAngle)
186 {
187 return false;
188 }
189 else
190 {
191 return true;
192 }
193 }
194 }
195}
196
197
198// Split hex (and hex only) along edgeI creating two prisms
199bool splitHex
200(
201 const polyMesh& mesh,
202 const label celli,
203 const label edgeI,
204
205 DynamicList<label>& cutCells,
206 DynamicList<labelList>& cellLoops,
207 DynamicList<scalarField>& cellEdgeWeights
208)
209{
210 // cut handling functions
211 edgeVertex ev(mesh);
212
213 const edgeList& edges = mesh.edges();
214 const faceList& faces = mesh.faces();
215
216 const edge& e = edges[edgeI];
217
218 // Get faces on the side, i.e. faces not using edge but still using one of
219 // the edge endpoints.
220
221 label leftI = -1;
222 label rightI = -1;
223 label leftFp = -1;
224 label rightFp = -1;
225
226 const cell& cFaces = mesh.cells()[celli];
227
228 for (const label facei : cFaces)
229 {
230 const face& f = faces[facei];
231
232 label fp0 = f.find(e[0]);
233 label fp1 = f.find(e[1]);
234
235 if (fp0 == -1)
236 {
237 if (fp1 != -1)
238 {
239 // Face uses e[1] but not e[0]
240 rightI = facei;
241 rightFp = fp1;
242
243 if (leftI != -1)
244 {
245 // Have both faces so exit
246 break;
247 }
248 }
249 }
250 else
251 {
252 if (fp1 != -1)
253 {
254 // Face uses both e[1] and e[0]
255 }
256 else
257 {
258 leftI = facei;
259 leftFp = fp0;
260
261 if (rightI != -1)
262 {
263 break;
264 }
265 }
266 }
267 }
268
269 if (leftI == -1 || rightI == -1)
270 {
272 << " rightI:" << rightI << abort(FatalError);
273 }
274
275 // Walk two vertices further on faces.
276
277 const face& leftF = faces[leftI];
278
279 label leftV = leftF[(leftFp + 2) % leftF.size()];
280
281 const face& rightF = faces[rightI];
282
283 label rightV = rightF[(rightFp + 2) % rightF.size()];
284
285 labelList loop(4);
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]);
290
291 scalarField loopWeights(4);
292 loopWeights[0] = -GREAT;
293 loopWeights[1] = -GREAT;
294 loopWeights[2] = -GREAT;
295 loopWeights[3] = -GREAT;
296
297 cutCells.append(celli);
298 cellLoops.append(loop);
299 cellEdgeWeights.append(loopWeights);
300
301 return true;
302}
303
304
305// Split celli along edgeI with a plane along halfNorm direction.
306bool splitCell
307(
308 const polyMesh& mesh,
309 const geomCellLooper& cellCutter,
310
311 const label celli,
312 const label edgeI,
313 const vector& halfNorm,
314
315 const boolList& vertIsCut,
316 const boolList& edgeIsCut,
317 const scalarField& edgeWeight,
318
319 DynamicList<label>& cutCells,
320 DynamicList<labelList>& cellLoops,
321 DynamicList<scalarField>& cellEdgeWeights
322)
323{
324 const edge& e = mesh.edges()[edgeI];
325
326 const vector eVec = e.unitVec(mesh.points());
327
328 vector planeN = eVec ^ halfNorm;
329
330 // Slightly tilt plane to make it not cut edges exactly
331 // halfway on fully regular meshes (since we want cuts
332 // to be snapped to vertices)
333 planeN += 0.01*halfNorm;
334 planeN.normalise();
335
336 // Define plane through edge
337 plane cutPlane(mesh.points()[e.start()], planeN);
338
339 labelList loop;
340 scalarField loopWeights;
341
342 if
343 (
344 cellCutter.cut
345 (
346 cutPlane,
347 celli,
348 vertIsCut,
349 edgeIsCut,
350 edgeWeight,
351 loop,
352 loopWeights
353 )
354 )
355 {
356 // Did manage to cut cell. Copy into overall list.
357 cutCells.append(celli);
358 cellLoops.append(loop);
359 cellEdgeWeights.append(loopWeights);
360
361 return true;
362 }
363
364 return false;
365}
366
367
368// Collects cuts for all cells in cellSet
369void collectCuts
370(
371 const polyMesh& mesh,
372 const geomCellLooper& cellCutter,
373 const bool geometry,
374 const scalar minCos,
375 const scalar minSin,
376 const cellSet& cellsToCut,
377
378 DynamicList<label>& cutCells,
379 DynamicList<labelList>& cellLoops,
380 DynamicList<scalarField>& cellEdgeWeights
381)
382{
383 // Get data from mesh
384 const cellShapeList& cellShapes = mesh.cellShapes();
385 const labelList& own = mesh.faceOwner();
386 const labelListList& cellEdges = mesh.cellEdges();
387 const vectorField& faceAreas = mesh.faceAreas();
388
389 // Hex shape
390 const cellModel& hex = cellModel::ref(cellModel::HEX);
391
392 // cut handling functions
393 edgeVertex ev(mesh);
394
395
396 // Cut information per mesh entity
397 boolList vertIsCut(mesh.nPoints(), false);
398 boolList edgeIsCut(mesh.nEdges(), false);
399 scalarField edgeWeight(mesh.nEdges(), -GREAT);
400
401 for (const label celli : cellsToCut)
402 {
403 const labelList& cEdges = cellEdges[celli];
404
405 for (const label edgeI : cEdges)
406 {
407 label f0, f1;
408 meshTools::getEdgeFaces(mesh, celli, edgeI, f0, f1);
409
410 const vector n0 = normalised(faceAreas[f0]);
411 const vector n1 = normalised(faceAreas[f1]);
412
413 if
414 (
415 largerAngle
416 (
417 mesh,
418 minCos,
419 minSin,
420
421 celli,
422 f0,
423 f1,
424 n0,
425 n1
426 )
427 )
428 {
429 bool splitOk = false;
430
431 if (!geometry && cellShapes[celli].model() == hex)
432 {
433 splitOk =
434 splitHex
435 (
436 mesh,
437 celli,
438 edgeI,
439
440 cutCells,
441 cellLoops,
442 cellEdgeWeights
443 );
444 }
445 else
446 {
447 vector halfNorm;
448
449 if ((own[f0] == celli) ^ (own[f1] == celli))
450 {
451 // Opposite owner orientation
452 halfNorm = 0.5*(n0 - n1);
453 }
454 else
455 {
456 // Faces have same owner or same neighbour so
457 // normals point in same direction
458 halfNorm = 0.5*(n0 + n1);
459 }
460
461 splitOk =
463 (
464 mesh,
465 cellCutter,
466 celli,
467 edgeI,
468 halfNorm,
469
470 vertIsCut,
471 edgeIsCut,
472 edgeWeight,
473
474 cutCells,
475 cellLoops,
476 cellEdgeWeights
477 );
478 }
479
480
481 if (splitOk)
482 {
483 // Update cell/edge/vertex wise info.
484 label index = cellLoops.size() - 1;
485 const labelList& loop = cellLoops[index];
486 const scalarField& loopWeights = cellEdgeWeights[index];
487
488 forAll(loop, i)
489 {
490 label cut = loop[i];
491
492 if (ev.isEdge(cut))
493 {
494 edgeIsCut[ev.getEdge(cut)] = true;
495 edgeWeight[ev.getEdge(cut)] = loopWeights[i];
496 }
497 else
498 {
499 vertIsCut[ev.getVertex(cut)] = true;
500 }
501 }
502
503 // Stop checking edges for this cell.
504 break;
505 }
506 }
507 }
508 }
509
510 cutCells.shrink();
511 cellLoops.shrink();
512 cellEdgeWeights.shrink();
513}
514
515
516
517int main(int argc, char *argv[])
518{
519 argList::addNote
520 (
521 "Split cells with flat faces"
522 );
523 #include "addOverwriteOption.H"
524 argList::noParallel();
525 argList::addArgument
526 (
527 "edgeAngle",
528 "in degrees [0-360]"
529 );
530
531 argList::addOption
532 (
533 "set",
534 "name",
535 "Split cells from specified cellSet only"
536 );
537 argList::addBoolOption
538 (
539 "geometry",
540 "Use geometric cut for hexes as well"
541 );
542 argList::addOption
543 (
544 "tol",
545 "scalar",
546 "Edge snap tolerance (default 0.2)"
547 );
548
549 argList::noFunctionObjects(); // Never use function objects
550
551 #include "setRootCase.H"
552 #include "createTime.H"
553 #include "createPolyMesh.H"
554
555 const word oldInstance = mesh.pointsInstance();
556
557 const scalar featureAngle = args.get<scalar>(1);
558 const scalar minCos = Foam::cos(degToRad(featureAngle));
559 const scalar minSin = Foam::sin(degToRad(featureAngle));
560
561 const bool readSet = args.found("set");
562 const bool geometry = args.found("geometry");
563 const bool overwrite = args.found("overwrite");
564
565 const scalar edgeTol = args.getOrDefault<scalar>("tol", 0.2);
566
567 Info<< "Trying to split cells with internal angles > feature angle\n" << nl
568 << "featureAngle : " << featureAngle << nl
569 << "edge snapping tol : " << edgeTol << nl;
570 if (readSet)
571 {
572 Info<< "candidate cells : cellSet " << args["set"] << nl;
573 }
574 else
575 {
576 Info<< "candidate cells : all cells" << nl;
577 }
578 if (geometry)
579 {
580 Info<< "hex cuts : geometric; using edge tolerance" << nl;
581 }
582 else
583 {
584 Info<< "hex cuts : topological; cut to opposite edge" << nl;
585 }
586 Info<< endl;
587
588
589 // Cell circumference cutter
590 geomCellLooper cellCutter(mesh);
591 // Snap all edge cuts close to endpoints to vertices.
592 geomCellLooper::setSnapTol(edgeTol);
593
594 // Candidate cells to cut
595 cellSet cellsToCut(mesh, "cellsToCut", mesh.nCells()/100);
596
597 if (readSet)
598 {
599 // Read cells to cut from cellSet
600 cellSet cells(mesh, args["set"]);
601
602 cellsToCut = cells;
603 }
604
605 while (true)
606 {
607 if (!readSet)
608 {
609 // Try all cells for cutting
610 for (label celli = 0; celli < mesh.nCells(); celli++)
611 {
612 cellsToCut.insert(celli);
613 }
614 }
615
616
617 // Cut information per cut cell
618 DynamicList<label> cutCells(mesh.nCells()/10 + 10);
619 DynamicList<labelList> cellLoops(mesh.nCells()/10 + 10);
620 DynamicList<scalarField> cellEdgeWeights(mesh.nCells()/10 + 10);
621
622 collectCuts
623 (
624 mesh,
625 cellCutter,
626 geometry,
627 minCos,
628 minSin,
629 cellsToCut,
630
631 cutCells,
632 cellLoops,
633 cellEdgeWeights
634 );
635
636 cellSet cutSet(mesh, "cutSet", cutCells.size());
637 cutSet.insert(cutCells);
638
639 // Gets cuts across cells from cuts through edges.
640 Info<< "Writing " << cutSet.size() << " cells to cut to cellSet "
641 << cutSet.instance()/cutSet.local()/cutSet.name()
642 << endl << endl;
643 cutSet.write();
644
645 // Analyze cuts for clashes.
646 cellCuts cuts
647 (
648 mesh,
649 cutCells, // cells candidate for cutting
650 cellLoops,
651 cellEdgeWeights
652 );
653
654 Info<< "Actually cut cells:" << cuts.nLoops() << nl << endl;
655
656 if (cuts.nLoops() == 0)
657 {
658 break;
659 }
660
661 // Remove cut cells from cellsToCut (Note:only relevant if -readSet)
662 forAll(cuts.cellLoops(), celli)
663 {
664 if (cuts.cellLoops()[celli].size())
665 {
666 //Info<< "Removing cut cell " << celli << " from wishlist"
667 // << endl;
668 cellsToCut.erase(celli);
669 }
670 }
671
672 // At least some cells are cut.
673 polyTopoChange meshMod(mesh);
674
675 // Cutting engine
676 meshCutter cutter(mesh);
677
678 // Insert mesh refinement into polyTopoChange.
679 cutter.setRefinement(cuts, meshMod);
680
681 // Do all changes
682 Info<< "Morphing ..." << endl;
683
684 if (!overwrite)
685 {
686 ++runTime;
687 }
688
689 autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
690
691 if (morphMap().hasMotionPoints())
692 {
693 mesh.movePoints(morphMap().preMotionPoints());
694 }
695
696 // Update stored labels on meshCutter
697 cutter.updateMesh(morphMap());
698
699 // Update cellSet
700 cellsToCut.updateMesh(morphMap());
701
702 Info<< "Remaining:" << cellsToCut.size() << endl;
703
704 // Write resulting mesh
705 if (overwrite)
706 {
707 mesh.setInstance(oldInstance);
708 }
709
710 Info<< "Writing refined morphMesh to time " << runTime.timeName()
711 << endl;
712
713 mesh.write();
714 }
715
716 Info<< "End\n" << endl;
717
718 return 0;
719}
720
721
722// ************************************************************************* //
Various functions to operate on Lists.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:434
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
void setSize(const label n)
Alias for resize()
Definition: List.H:218
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
Description of cuts across cells.
Definition: cellCuts.H:113
Maps a geometry to a set of cell primitives.
Definition: cellModel.H:73
A collection of cell labels.
Definition: cellSet.H:54
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
Combines edge or vertex in single label. Used to specify cuts across cell circumference.
Definition: edgeVertex.H:56
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
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.
Cuts (splits) cells.
Definition: meshCutter.H:141
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
Definition: plane.H:95
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
Direct mesh changes based on v1.3 polyTopoChange syntax.
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:79
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....
Definition: splitCell.H:52
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.
Definition: word.H:68
const cellModel & hex
cellShapeList cellShapes
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const cellShapeList & cells
Namespace for OpenFOAM.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:680
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.
Definition: Ostream.H:372
errorManip< error > abort(error &err)
Definition: errorManip.H:144
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
Foam::argList args(argc, argv)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Unit conversion functions.