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