DelaunayMeshIO.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) 2012-2017 OpenFOAM Foundation
9 Copyright (C) 2019 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 "DelaunayMesh.H"
30#include "fvMesh.H"
31#include "pointConversion.H"
32#include "wallPolyPatch.H"
33#include "processorPolyPatch.H"
34#include "labelIOField.H"
35
36// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37
38template<class Triangulation>
40(
41 faceList& faces,
42 labelList& owner,
43 labelList& neighbour
44) const
45{
46 // Upper triangular order:
47 // + owner is sorted in ascending cell order
48 // + within each block of equal value for owner, neighbour is sorted in
49 // ascending cell order.
50 // + faces sorted to correspond
51 // e.g.
52 // owner | neighbour
53 // 0 | 2
54 // 0 | 23
55 // 0 | 71
56 // 1 | 23
57 // 1 | 24
58 // 1 | 91
59
60 List<labelPair> ownerNeighbourPair(owner.size());
61
62 forAll(ownerNeighbourPair, oNI)
63 {
64 ownerNeighbourPair[oNI] = labelPair(owner[oNI], neighbour[oNI]);
65 }
66
67 Info<< nl
68 << "Sorting faces, owner and neighbour into upper triangular order"
69 << endl;
70
71 labelList oldToNew(sortedOrder(ownerNeighbourPair));
72
73 oldToNew = invert(oldToNew.size(), oldToNew);
74
75 inplaceReorder(oldToNew, faces);
76 inplaceReorder(oldToNew, owner);
77 inplaceReorder(oldToNew, neighbour);
78}
79
80
81template<class Triangulation>
83(
84 const label nInternalFaces,
85 faceList& faces,
86 labelList& owner,
87 PtrList<dictionary>& patchDicts,
88 const List<DynamicList<face>>& patchFaces,
89 const List<DynamicList<label>>& patchOwners
90) const
91{
92 label nPatches = patchFaces.size();
93
94 patchDicts.setSize(nPatches);
95 forAll(patchDicts, patchi)
96 {
97 patchDicts.set(patchi, new dictionary());
98 }
99
100 label nBoundaryFaces = 0;
101
102 forAll(patchFaces, p)
103 {
104 patchDicts[p].set("nFaces", patchFaces[p].size());
105 patchDicts[p].set("startFace", nInternalFaces + nBoundaryFaces);
106
107 nBoundaryFaces += patchFaces[p].size();
108 }
109
110 faces.setSize(nInternalFaces + nBoundaryFaces);
111 owner.setSize(nInternalFaces + nBoundaryFaces);
112
113 label facei = nInternalFaces;
114
115 forAll(patchFaces, p)
116 {
117 forAll(patchFaces[p], f)
118 {
119 faces[facei] = patchFaces[p][f];
120 owner[facei] = patchOwners[p][f];
121
122 facei++;
123 }
124 }
125}
126
127
128// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
129
130template<class Triangulation>
132{
133 PrintTable<word, label> triInfoTable("Mesh Statistics");
134
135 triInfoTable.add("Points", Triangulation::number_of_vertices());
136 triInfoTable.add("Edges", Triangulation::number_of_finite_edges());
137 triInfoTable.add("Faces", Triangulation::number_of_finite_facets());
138 triInfoTable.add("Cells", Triangulation::number_of_finite_cells());
139
140 scalar minSize = GREAT;
141 scalar maxSize = 0;
142
143 for
144 (
145 Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
146 vit != Triangulation::finite_vertices_end();
147 ++vit
148 )
149 {
150 // Only internal or boundary vertices have a size
151 if (vit->internalOrBoundaryPoint())
152 {
153 minSize = min(vit->targetCellSize(), minSize);
154 maxSize = max(vit->targetCellSize(), maxSize);
155 }
156 }
157
159 triInfoTable.print(Info, true, true);
160
161 Info<< "Size (Min/Max) = "
162 << returnReduce(minSize, minOp<scalar>()) << " "
163 << returnReduce(maxSize, maxOp<scalar>()) << endl;
164
166}
167
168
169template<class Triangulation>
171{
172 label nInternal = 0;
173 label nInternalRef = 0;
174 label nUnassigned = 0;
175 label nUnassignedRef = 0;
176 label nInternalNearBoundary = 0;
177 label nInternalNearBoundaryRef = 0;
178 label nInternalSurface = 0;
179 label nInternalSurfaceRef = 0;
180 label nInternalFeatureEdge = 0;
181 label nInternalFeatureEdgeRef = 0;
182 label nInternalFeaturePoint = 0;
183 label nInternalFeaturePointRef = 0;
184 label nExternalSurface = 0;
185 label nExternalSurfaceRef = 0;
186 label nExternalFeatureEdge = 0;
187 label nExternalFeatureEdgeRef = 0;
188 label nExternalFeaturePoint = 0;
189 label nExternalFeaturePointRef = 0;
190 label nFar = 0;
191 label nReferred = 0;
192
193 for
194 (
195 Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
196 vit != Triangulation::finite_vertices_end();
197 ++vit
198 )
199 {
200 if (vit->type() == Vb::vtInternal)
201 {
202 if (vit->referred())
203 {
204 nReferred++;
205 nInternalRef++;
206 }
207
208 nInternal++;
209 }
210 else if (vit->type() == Vb::vtUnassigned)
211 {
212 if (vit->referred())
213 {
214 nReferred++;
215 nUnassignedRef++;
216 }
217
218 nUnassigned++;
219 }
220 else if (vit->type() == Vb::vtInternalNearBoundary)
221 {
222 if (vit->referred())
223 {
224 nReferred++;
225 nInternalNearBoundaryRef++;
226 }
227
228 nInternalNearBoundary++;
229 }
230 else if (vit->type() == Vb::vtInternalSurface)
231 {
232 if (vit->referred())
233 {
234 nReferred++;
235 nInternalSurfaceRef++;
236 }
237
238 nInternalSurface++;
239 }
240 else if (vit->type() == Vb::vtInternalFeatureEdge)
241 {
242 if (vit->referred())
243 {
244 nReferred++;
245 nInternalFeatureEdgeRef++;
246 }
247
248 nInternalFeatureEdge++;
249 }
250 else if (vit->type() == Vb::vtInternalFeaturePoint)
251 {
252 if (vit->referred())
253 {
254 nReferred++;
255 nInternalFeaturePointRef++;
256 }
257
258 nInternalFeaturePoint++;
259 }
260 else if (vit->type() == Vb::vtExternalSurface)
261 {
262 if (vit->referred())
263 {
264 nReferred++;
265 nExternalSurfaceRef++;
266 }
267
268 nExternalSurface++;
269 }
270 else if (vit->type() == Vb::vtExternalFeatureEdge)
271 {
272 if (vit->referred())
273 {
274 nReferred++;
275 nExternalFeatureEdgeRef++;
276 }
277
278 nExternalFeatureEdge++;
279 }
280 else if (vit->type() == Vb::vtExternalFeaturePoint)
281 {
282 if (vit->referred())
283 {
284 nReferred++;
285 nExternalFeaturePointRef++;
286 }
287
288 nExternalFeaturePoint++;
289 }
290 else if (vit->type() == Vb::vtFar)
291 {
292 nFar++;
293 }
294 }
295
296 label nTotalVertices =
297 nUnassigned
298 + nInternal
299 + nInternalNearBoundary
300 + nInternalSurface
301 + nInternalFeatureEdge
302 + nInternalFeaturePoint
303 + nExternalSurface
304 + nExternalFeatureEdge
305 + nExternalFeaturePoint
306 + nFar;
307
308 if (nTotalVertices != label(Triangulation::number_of_vertices()))
309 {
311 << nTotalVertices << " does not equal "
312 << Triangulation::number_of_vertices()
313 << endl;
314 }
315
316 PrintTable<word, label> vertexTable("Vertex Type Information");
317
318 vertexTable.add("Total", nTotalVertices);
319 vertexTable.add("Unassigned", nUnassigned);
320 vertexTable.add("nInternal", nInternal);
321 vertexTable.add("nInternalNearBoundary", nInternalNearBoundary);
322 vertexTable.add("nInternalSurface", nInternalSurface);
323 vertexTable.add("nInternalFeatureEdge", nInternalFeatureEdge);
324 vertexTable.add("nInternalFeaturePoint", nInternalFeaturePoint);
325 vertexTable.add("nExternalSurface", nExternalSurface);
326 vertexTable.add("nExternalFeatureEdge", nExternalFeatureEdge);
327 vertexTable.add("nExternalFeaturePoint", nExternalFeaturePoint);
328 vertexTable.add("nFar", nFar);
329 vertexTable.add("nReferred", nReferred);
330
331 os << endl;
332 vertexTable.print(os);
333}
334
335
336template<class Triangulation>
339(
340 const fileName& name,
341 labelPairLookup& vertexMap,
342 labelList& cellMap,
343 const bool writeDelaunayData
344) const
345{
346 pointField points(Triangulation::number_of_vertices());
347 faceList faces(Triangulation::number_of_finite_facets());
348 labelList owner(Triangulation::number_of_finite_facets());
349 labelList neighbour(Triangulation::number_of_finite_facets());
350
351 wordList patchNames(1, "foamyHexMesh_defaultPatch");
352 wordList patchTypes(1, wallPolyPatch::typeName);
353
354 PtrList<dictionary> patchDicts(1);
355 patchDicts.set(0, new dictionary());
356
357 List<DynamicList<face>> patchFaces(1, DynamicList<face>());
358 List<DynamicList<label>> patchOwners(1, DynamicList<label>());
359
360 vertexMap.resize(vertexCount());
361 cellMap.setSize(Triangulation::number_of_finite_cells(), -1);
362
363 // Calculate pts and a map of point index to location in pts.
364 label vertI = 0;
365
366// labelIOField indices
367// (
368// IOobject
369// (
370// "indices",
371// time().timeName(),
372// name,
373// time(),
374// IOobject::NO_READ,
375// IOobject::AUTO_WRITE
376// ),
377// Triangulation::number_of_vertices()
378// );
379
380 labelIOField types
381 (
382 IOobject
383 (
384 "types",
385 time().timeName(),
386 name,
387 time(),
388 IOobject::NO_READ,
389 IOobject::AUTO_WRITE
390 ),
391 label(Triangulation::number_of_vertices())
392 );
393
394 labelIOField processorIndices
395 (
396 IOobject
397 (
398 "processorIndices",
399 time().timeName(),
400 name,
401 time(),
402 IOobject::NO_READ,
403 IOobject::AUTO_WRITE
404 ),
405 label(Triangulation::number_of_vertices())
406 );
407
408 for
409 (
410 Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
411 vit != Triangulation::finite_vertices_end();
412 ++vit
413 )
414 {
415 if (!vit->farPoint())
416 {
417 vertexMap(labelPair(vit->index(), vit->procIndex())) = vertI;
418 points[vertI] = topoint(vit->point());
419// indices[vertI] = vit->index();
420 types[vertI] = static_cast<label>(vit->type());
421 processorIndices[vertI] = vit->procIndex();
422 vertI++;
423 }
424 }
425
426 points.setSize(vertI);
427// indices.setSize(vertI);
428 types.setSize(vertI);
429 processorIndices.setSize(vertI);
430
431
432 // Index the cells
433 label celli = 0;
434
435 for
436 (
437 Finite_cells_iterator cit = Triangulation::finite_cells_begin();
438 cit != Triangulation::finite_cells_end();
439 ++cit
440 )
441 {
442 if
443 (
444 !cit->hasFarPoint()
445 && !Triangulation::is_infinite(cit)
446 && cit->real()
447 )
448 {
449 cellMap[cit->cellIndex()] = celli++;
450 }
451 }
452
453 label facei = 0;
454 labelList verticesOnTriFace(3, label(-1));
455 face newFace(verticesOnTriFace);
456
457 for
458 (
459 Finite_facets_iterator fit = Triangulation::finite_facets_begin();
460 fit != Triangulation::finite_facets_end();
461 ++fit
462 )
463 {
464 const Cell_handle c1(fit->first);
465 const label oppositeVertex = fit->second;
466 const Cell_handle c2(c1->neighbor(oppositeVertex));
467
468 // Do not output if face has neither opposite vertex as an internal
469// if
470// (
471// !c1->vertex(oppositeVertex)->internalPoint()
472// || !Triangulation::mirror_vertex(c1, oppositeVertex)->internalPoint()
473// )
474// {
475// continue;
476// }
477
478 label c1I = Cb::ctFar;
479 bool c1Real = false;
480 if
481 (
482 !Triangulation::is_infinite(c1)
483 && !c1->hasFarPoint()
484 && c1->real()
485 )
486 {
487 c1I = cellMap[c1->cellIndex()];
488 c1Real = true;
489 }
490
491 label c2I = Cb::ctFar;
492 bool c2Real = false;
493 if
494 (
495 !Triangulation::is_infinite(c2)
496 && !c2->hasFarPoint()
497 && c2->real()
498 )
499 {
500 c2I = cellMap[c2->cellIndex()];
501 c2Real = true;
502 }
503
504 if (!c1Real && !c2Real)
505 {
506 // Both tets are outside, skip
507 continue;
508 }
509
510 label ownerCell = -1;
511 label neighbourCell = -1;
512
513 for (label i = 0; i < 3; i++)
514 {
515 verticesOnTriFace[i] = vertexMap
516 [
518 (
519 c1->vertex
520 (
521 Triangulation::vertex_triple_index(oppositeVertex, i)
522 )->index(),
523 c1->vertex
524 (
525 Triangulation::vertex_triple_index(oppositeVertex, i)
526 )->procIndex()
527 )
528 ];
529 }
530
531 newFace = face(verticesOnTriFace);
532
533 if (!c1Real || !c2Real)
534 {
535 // Boundary face...
536 if (!c1Real)
537 {
538 //... with c1 outside
539 ownerCell = c2I;
540 }
541 else
542 {
543 // ... with c2 outside
544 ownerCell = c1I;
545
546 reverse(newFace);
547 }
548
549 patchFaces[0].append(newFace);
550 patchOwners[0].append(ownerCell);
551 }
552 else
553 {
554 // Internal face...
555 if (c1I < c2I)
556 {
557 // ...with c1 as the ownerCell
558 ownerCell = c1I;
559 neighbourCell = c2I;
560
561 reverse(newFace);
562 }
563 else
564 {
565 // ...with c2 as the ownerCell
566 ownerCell = c2I;
567 neighbourCell = c1I;
568 }
569
570 faces[facei] = newFace;
571 owner[facei] = ownerCell;
572 neighbour[facei] = neighbourCell;
573 facei++;
574 }
575 }
576
577 faces.setSize(facei);
578 owner.setSize(facei);
579 neighbour.setSize(facei);
580
581 sortFaces(faces, owner, neighbour);
582
583 Info<< "Creating patches" << endl;
584
585 addPatches
586 (
587 facei,
588 faces,
589 owner,
591 patchFaces,
592 patchOwners
593 );
594
595 Info<< "Creating mesh" << endl;
596
597 auto meshPtr = autoPtr<polyMesh>::New
598 (
599 IOobject
600 (
601 name,
602 time().timeName(),
603 time(),
604 IOobject::NO_READ,
605 IOobject::NO_WRITE
606 ),
607 std::move(points),
608 std::move(faces),
609 std::move(owner),
610 std::move(neighbour)
611 );
612
613 Info<< "Adding patches" << endl;
614
615 List<polyPatch*> patches(patchNames.size());
616
617 label nValidPatches = 0;
618
620 {
621 patches[nValidPatches] = polyPatch::New
622 (
623 patchTypes[p],
624 patchNames[p],
625 patchDicts[p],
626 nValidPatches,
627 meshPtr().boundaryMesh()
628 ).ptr();
629
630 nValidPatches++;
631 }
632
633 patches.setSize(nValidPatches);
634
635 meshPtr().addPatches(patches);
636
637 if (writeDelaunayData)
638 {
639// indices.write();
640 types.write();
641 processorIndices.write();
642 }
643
644 Info<< "Mesh created" << endl;
645
646 return meshPtr;
647}
648
649
650// ************************************************************************* //
Y[inertIndex] max(0.0)
The vertex and cell classes must have an index defined.
Definition: DelaunayMesh.H:63
virtual void print(Ostream &os) const
Print stream description to Ostream.
Definition: IOstream.C:80
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
void printInfo() const
Print general information about the mesh.
Definition: ccmReader.C:522
Switch printVertexInfo() const
Return the printVertexInfo Switch.
Definition: cvControlsI.H:163
volScalarField & p
const polyBoundaryMesh & patches
Foam::autoPtr< Foam::dynamicFvMesh > meshPtr
OBJstream os(runTime.globalPath()/outputName)
const label nPatches
const pointField & points
word timeName
Definition: getTimeIndex.H:3
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:57
List< word > wordList
A List of words.
Definition: fileName.H:63
List< label > labelList
A List of labels.
Definition: List.H:66
pointFromPoint topoint(const Point &P)
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
Definition: UListI.H:449
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
IOField< label > labelIOField
labelField with IO.
Definition: labelIOField.H:44
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:349
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
Definition: ListOps.C:36
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:356
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
wordList patchTypes(nPatches)
wordList patchNames(nPatches)
PtrList< dictionary > patchDicts
Definition: readKivaGrid.H:532
labelList f(nPoints)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333