shortEdgeFilter2D.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) 2013-2016 OpenFOAM Foundation
9 Copyright (C) 2020-2021 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
27\*---------------------------------------------------------------------------*/
28
29#include "shortEdgeFilter2D.H"
30
31namespace Foam
32{
33 defineTypeNameAndDebug(shortEdgeFilter2D, 0);
34}
35
36
37// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
38
39void Foam::shortEdgeFilter2D::assignBoundaryPointRegions
40(
41 List<DynamicList<label>>& boundaryPointRegions
42) const
43{
44 forAllConstIters(mapEdgesRegion_, iter)
45 {
46 const edge& e = iter.key();
47 const label regi = iter.val();
48
49 boundaryPointRegions[e.start()].appendUniq(regi);
50 boundaryPointRegions[e.end()].appendUniq(regi);
51 }
52}
53
54
55void Foam::shortEdgeFilter2D::updateEdgeRegionMap
56(
57 const MeshedSurface<face>& surfMesh,
58 const List<DynamicList<label>>& boundaryPtRegions,
59 const labelList& surfPtToBoundaryPt,
60 EdgeMap<label>& mapEdgesRegion,
61 labelList& patchSizes
62) const
63{
64 EdgeMap<label> newMapEdgesRegion(mapEdgesRegion.size());
65
66 const edgeList& edges = surfMesh.edges();
67 const labelList& meshPoints = surfMesh.meshPoints();
68
69 patchSizes.setSize(patchNames_.size(), 0);
70 patchSizes = 0;
71
72 forAll(edges, edgeI)
73 {
74 if (surfMesh.isInternalEdge(edgeI))
75 {
76 continue;
77 }
78
79 const edge& e = edges[edgeI];
80
81 const label startI = meshPoints[e[0]];
82 const label endI = meshPoints[e[1]];
83
84 label region = -1;
85
86 const DynamicList<label> startPtRegions =
87 boundaryPtRegions[surfPtToBoundaryPt[startI]];
88 const DynamicList<label> endPtRegions =
89 boundaryPtRegions[surfPtToBoundaryPt[endI]];
90
91 if (startPtRegions.size() > 1 && endPtRegions.size() > 1)
92 {
93 region = startPtRegions[0];
94
96 << "Both points in edge are in different regions."
97 << " Assigning edge to region " << region
98 << endl;
99 }
100 else if (startPtRegions.size() > 1 || endPtRegions.size() > 1)
101 {
102 region =
103 (
104 startPtRegions.size() > 1
105 ? endPtRegions[0]
106 : startPtRegions[0]
107 );
108 }
109 else if
110 (
111 startPtRegions[0] == endPtRegions[0]
112 && startPtRegions[0] != -1
113 )
114 {
115 region = startPtRegions[0];
116 }
117
118 if (region != -1)
119 {
120 newMapEdgesRegion.insert(e, region);
121 patchSizes[region]++;
122 }
123 }
124
125 mapEdgesRegion.transfer(newMapEdgesRegion);
126}
127
128
129// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
130
132(
133 const Foam::CV2D& cv2Dmesh,
134 const dictionary& dict
135)
136:
137 shortEdgeFilterFactor_(dict.get<scalar>("shortEdgeFilterFactor")),
138 edgeAttachedToBoundaryFactor_
139 (
140 dict.getOrDefault<scalar>("edgeAttachedToBoundaryFactor", 2.0)
141 ),
142 patchNames_(),
143 patchSizes_(),
144 mapEdgesRegion_(),
145 indirectPatchEdge_()
146{
147 point2DField points2D;
148 faceList faces;
149
150 cv2Dmesh.calcDual
151 (
152 points2D,
153 faces,
154 patchNames_,
155 patchSizes_,
156 mapEdgesRegion_,
157 indirectPatchEdge_
158 );
159
160 pointField points(points2D.size());
161 forAll(points, ip)
162 {
163 points[ip] = cv2Dmesh.toPoint3D(points2D[ip]);
164 }
165
166 if (debug)
167 {
168 OFstream str("indirectPatchEdges.obj");
169 label count = 0;
170
171 Info<< "Writing indirectPatchEdges to " << str.name() << endl;
172
173 forAllConstIters(indirectPatchEdge_, iter)
174 {
175 const edge& e = iter.key();
176
177 meshTools::writeOBJ
178 (
179 str,
180 points[e.start()],
181 points[e.end()],
182 count
183 );
184 }
185 }
186
187 points2D.clear();
188
189 ms_ = MeshedSurface<face>(std::move(points), std::move(faces));
190
191 Info<< "Meshed surface stats before edge filtering :" << endl;
192 ms_.writeStats(Info);
193
194 if (debug)
195 {
196 writeInfo(Info);
197
198 ms_.write("MeshedSurface_preFilter.obj");
199 }
200}
201
202
203// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
204
206{
207 // These are global indices.
208 const pointField& points = ms_.points();
209 const edgeList& edges = ms_.edges();
210 const faceList& faces = ms_.surfFaces();
211 const labelList& meshPoints = ms_.meshPoints();
212 const labelList& boundaryPoints = ms_.boundaryPoints();
213
214 label maxChain = 0;
215 label nPointsToRemove = 0;
216
217 labelList pointsToRemove(ms_.points().size(), -1);
218
219 // List of number of vertices in a face.
220 labelList newFaceVertexCount(faces.size(), -1);
221 forAll(faces, facei)
222 {
223 newFaceVertexCount[facei] = faces[facei].size();
224 }
225
226 // Check if the point is a boundary point. Flag if it is so that
227 // it will not be deleted.
228 List<DynamicList<label>> boundaryPointRegions
229 (
230 points.size(),
231 DynamicList<label>()
232 );
233 assignBoundaryPointRegions(boundaryPointRegions);
234
235 // Check if an edge has a boundary point. It it does the edge length
236 // will be doubled when working out its length.
237 Info<< " Marking edges attached to boundaries." << endl;
238 boolList edgeAttachedToBoundary(edges.size(), false);
239 forAll(edges, edgeI)
240 {
241 const edge& e = edges[edgeI];
242 const label startVertex = e.start();
243 const label endVertex = e.end();
244
245 forAll(boundaryPoints, bPoint)
246 {
247 if
248 (
249 boundaryPoints[bPoint] == startVertex
250 || boundaryPoints[bPoint] == endVertex
251 )
252 {
253 edgeAttachedToBoundary[edgeI] = true;
254 }
255 }
256 }
257
258 forAll(edges, edgeI)
259 {
260 const edge& e = edges[edgeI];
261
262 // get the vertices of that edge.
263 const label startVertex = e.start();
264 const label endVertex = e.end();
265
266 scalar edgeLength =
267 mag
268 (
269 points[meshPoints[startVertex]]
270 - points[meshPoints[endVertex]]
271 );
272
273 if (edgeAttachedToBoundary[edgeI])
274 {
275 edgeLength *= edgeAttachedToBoundaryFactor_;
276 }
277
278 scalar shortEdgeFilterValue = 0.0;
279
280 const labelList& psEdges = ms_.pointEdges()[startVertex];
281 const labelList& peEdges = ms_.pointEdges()[endVertex];
282
283 forAll(psEdges, psEdgeI)
284 {
285 const edge& psE = edges[psEdges[psEdgeI]];
286 if (edgeI != psEdges[psEdgeI])
287 {
288 shortEdgeFilterValue +=
289 mag
290 (
291 points[meshPoints[psE.start()]]
292 - points[meshPoints[psE.end()]]
293 );
294 }
295 }
296
297 forAll(peEdges, peEdgeI)
298 {
299 const edge& peE = edges[peEdges[peEdgeI]];
300 if (edgeI != peEdges[peEdgeI])
301 {
302 shortEdgeFilterValue +=
303 mag
304 (
305 points[meshPoints[peE.start()]]
306 - points[meshPoints[peE.end()]]
307 );
308 }
309 }
310
311 shortEdgeFilterValue *=
312 shortEdgeFilterFactor_
313 /(psEdges.size() + peEdges.size() - 2);
314
315 edge lookupInPatchEdge
316 (
317 meshPoints[startVertex],
318 meshPoints[endVertex]
319 );
320
321 if
322 (
323 edgeLength < shortEdgeFilterValue
324 || indirectPatchEdge_.found(lookupInPatchEdge)
325 )
326 {
327 bool flagDegenerateFace = false;
328 const labelList& pFaces = ms_.pointFaces()[startVertex];
329
330 forAll(pFaces, pFacei)
331 {
332 const face& f = ms_.localFaces()[pFaces[pFacei]];
333 forAll(f, fp)
334 {
335 // If the edge is part of this face...
336 if (f[fp] == endVertex)
337 {
338 // If deleting vertex would create a triangle, don't!
339 if (newFaceVertexCount[pFaces[pFacei]] < 4)
340 {
341 flagDegenerateFace = true;
342 }
343 else
344 {
345 newFaceVertexCount[pFaces[pFacei]]--;
346 }
347 }
348 // If the edge is not part of this face...
349 else
350 {
351 // Deleting vertex of a triangle...
352 if (newFaceVertexCount[pFaces[pFacei]] < 3)
353 {
354 flagDegenerateFace = true;
355 }
356 }
357 }
358 }
359
360 // This if statement determines whether a point should be deleted.
361 if
362 (
363 pointsToRemove[meshPoints[startVertex]] == -1
364 && pointsToRemove[meshPoints[endVertex]] == -1
365 && !flagDegenerateFace
366 )
367 {
368 const DynamicList<label>& startVertexRegions =
369 boundaryPointRegions[meshPoints[startVertex]];
370 const DynamicList<label>& endVertexRegions =
371 boundaryPointRegions[meshPoints[endVertex]];
372
373 if (startVertexRegions.size() && endVertexRegions.size())
374 {
375 if (startVertexRegions.size() > 1)
376 {
377 pointsToRemove[meshPoints[endVertex]] =
378 meshPoints[startVertex];
379 }
380 else
381 {
382 pointsToRemove[meshPoints[startVertex]] =
383 meshPoints[endVertex];
384 }
385 }
386 else if (startVertexRegions.size())
387 {
388 pointsToRemove[meshPoints[endVertex]] =
389 meshPoints[startVertex];
390 }
391 else
392 {
393 pointsToRemove[meshPoints[startVertex]] =
394 meshPoints[endVertex];
395 }
396
397 ++nPointsToRemove;
398 }
399 }
400 }
401
402 label totalNewPoints = points.size() - nPointsToRemove;
403
404 pointField newPoints(totalNewPoints, Zero);
405 labelList newPointNumbers(points.size(), -1);
406 label numberRemoved = 0;
407
408 // Maintain addressing from new to old point field
409 labelList newPtToOldPt(totalNewPoints, -1);
410
411 forAll(points, pointi)
412 {
413 // If the point is NOT going to be removed.
414 if (pointsToRemove[pointi] == -1)
415 {
416 newPoints[pointi - numberRemoved] = points[pointi];
417 newPointNumbers[pointi] = pointi - numberRemoved;
418 newPtToOldPt[pointi - numberRemoved] = pointi;
419 }
420 else
421 {
422 numberRemoved++;
423 }
424 }
425
426 // Need a new faceList
427 faceList newFaces(faces.size());
428 label newFacei = 0;
429
430 labelList newFace;
431 label newFaceSize = 0;
432
433 // Now need to iterate over the faces and remove points. Global index.
434 forAll(faces, facei)
435 {
436 const face& f = faces[facei];
437
438 newFace.clear();
439 newFace.setSize(f.size());
440 newFaceSize = 0;
441
442 forAll(f, fp)
443 {
444 label pointi = f[fp];
445 // If not removing the point, then add it to the new face.
446 if (pointsToRemove[pointi] == -1)
447 {
448 newFace[newFaceSize++] = newPointNumbers[pointi];
449 }
450 else
451 {
452 label newPointi = pointsToRemove[pointi];
453 // Replace deleted point with point that it is being
454 // collapsed to.
455 if
456 (
457 f.nextLabel(fp) != newPointi
458 && f.prevLabel(fp) != newPointi
459 )
460 {
461 label pChain = newPointi;
462 label totalChain = 0;
463 for (label nChain = 0; nChain <= totalChain; ++nChain)
464 {
465 if (newPointNumbers[pChain] != -1)
466 {
467 newFace[newFaceSize++] = newPointNumbers[pChain];
468 newPointNumbers[pointi] = newPointNumbers[pChain];
469 maxChain = max(totalChain, maxChain);
470 }
471 else
472 {
474 << "Point " << pChain
475 << " marked for deletion as well as point "
476 << pointi << nl
477 << " Incrementing maxChain by 1 from "
478 << totalChain << " to " << totalChain + 1
479 << endl;
480 totalChain++;
481 }
482 pChain = pointsToRemove[pChain];
483 }
484 }
485 else
486 {
487 if (newPointNumbers[newPointi] != -1)
488 {
489 newPointNumbers[pointi] = newPointNumbers[newPointi];
490 }
491 }
492 }
493 }
494
495 newFace.setSize(newFaceSize);
496
497 if (newFace.size() > 2)
498 {
499 newFaces[newFacei++] = face(newFace);
500 }
501 else
502 {
504 << "Only " << newFace.size() << " in face " << facei
505 << exit(FatalError);
506 }
507 }
508
509 newFaces.setSize(newFacei);
510
511 MeshedSurface<face> fMesh(std::move(newPoints), std::move(newFaces));
512
513 updateEdgeRegionMap
514 (
515 fMesh,
516 boundaryPointRegions,
517 newPtToOldPt,
518 mapEdgesRegion_,
519 patchSizes_
520 );
521
522 forAll(newPointNumbers, pointi)
523 {
524 if (newPointNumbers[pointi] == -1)
525 {
527 << pointi << " will be deleted and " << newPointNumbers[pointi]
528 << ", so it will not be replaced. "
529 << "This will cause edges to be deleted." << endl;
530 }
531 }
532
533 ms_.transfer(fMesh);
534
535 if (debug)
536 {
537 Info<< " Maximum number of chained collapses = " << maxChain << endl;
538
540 }
541}
542
543
545{
546 os << "Short Edge Filtering Information:" << nl
547 << " shortEdgeFilterFactor: " << shortEdgeFilterFactor_ << nl
548 << " edgeAttachedToBoundaryFactor: " << edgeAttachedToBoundaryFactor_
549 << endl;
550
551 forAll(patchNames_, patchi)
552 {
553 os << " Patch " << patchNames_[patchi]
554 << ", size " << patchSizes_[patchi] << endl;
555 }
556
557 os << " There are " << mapEdgesRegion_.size()
558 << " boundary edges." << endl;
559
560 os << " Mesh Info:" << nl
561 << " Points: " << ms_.nPoints() << nl
562 << " Faces: " << ms_.size() << nl
563 << " Edges: " << ms_.nEdges() << nl
564 << " Internal: " << ms_.nInternalEdges() << nl
565 << " External: " << ms_.nEdges() - ms_.nInternalEdges()
566 << endl;
567}
568
569
570// ************************************************************************* //
Conformal-Voronoi 2D automatic mesher with grid or read initial points and point position relaxation ...
Definition: CV2D.H:149
void calcDual(point2DField &dualPoints, faceList &dualFaces, wordList &patchNames, labelList &patchSizes, EdgeMap< label > &mapEdgesRegion, EdgeMap< label > &indirectPatchEdge) const
Calculates dual points (circumcentres of tets) and faces.
Foam::point toPoint3D(const point2D &) const
Definition: CV2DI.H:142
bool found(const Key &key) const
Return true if hashed entry is found in table.
Definition: HashTableI.H:100
void setSize(const label n)
Alias for resize()
Definition: List.H:218
label appendUniq(const T &val)
Append an element if not already in the list.
Definition: ListI.H:232
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
void transfer(pointField &pointLst, List< Face > &faceLst)
Transfer the components.
const List< Face > & surfFaces() const
Return const access to the faces.
const labelListList & pointEdges() const
Return point-edge addressing.
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 > & points() const noexcept
Return reference to global points.
const labelList & boundaryPoints() const
Return list of boundary points, address into LOCAL point list.
const labelListList & pointFaces() const
Return point-face addressing.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
iterator end() noexcept
Return an iterator to end traversing the UList.
Definition: UListI.H:350
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
This class filters short edges generated by the CV2D mesher.
const MeshedSurface< face > & fMesh() const
void writeInfo(Ostream &os)
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
const pointField & points
#define WarningInFunction
Report a warning using Foam::Warning.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
Definition: BitOps.H:78
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
Namespace for OpenFOAM.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
List< label > labelList
A List of labels.
Definition: List.H:66
vector2DField point2DField
point2DField is a vector2DField.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< bool > boolList
A List of bools.
Definition: List.H:64
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::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].reset(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, &oldCyclicPolyPatch::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]
Definition: readKivaGrid.H:235
label newPointi
Definition: readKivaGrid.H:496
labelList f(nPoints)
dictionary dict
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278