tetDecomposer.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-2022 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 "tetDecomposer.H"
27 #include "meshTools.H"
28 #include "polyMesh.H"
29 #include "polyTopoChange.H"
30 #include "polyTopoChangeMap.H"
31 #include "OFstream.H"
32 #include "EdgeMap.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(tetDecomposer, 0);
39 
40  template<>
42  {
43  "faceCentre",
44  "faceDiagonal"
45  };
46 
49 }
50 
51 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
52 
53 void Foam::tetDecomposer::modifyFace
54 (
55  polyTopoChange& meshMod,
56  const face& f,
57  const label facei,
58  const label own,
59  const label nei,
60  const label patchi,
61  const label zoneI,
62  const bool zoneFlip
63 ) const
64 {
65  // First usage of face. Modify.
66  if (nei == -1 || own < nei)
67  {
68  meshMod.modifyFace
69  (
70  f, // modified face
71  facei, // label of face
72  own, // owner
73  nei, // neighbour
74  false, // face flip
75  patchi, // patch for face
76  zoneI, // zone for face
77  zoneFlip // face flip in zone
78  );
79  }
80  else
81  {
82  meshMod.modifyFace
83  (
84  f.reverseFace(), // modified face
85  facei, // label of face
86  nei, // owner
87  own, // neighbour
88  true, // face flip
89  patchi, // patch for face
90  zoneI, // zone for face
91  !zoneFlip // face flip in zone
92  );
93  }
94 }
95 
96 
97 void Foam::tetDecomposer::addFace
98 (
99  polyTopoChange& meshMod,
100  const face& f,
101  const label own,
102  const label nei,
103  const label masterPointID,
104  const label masterEdgeID,
105  const label masterFaceID,
106  const label patchi,
107  const label zoneI,
108  const bool zoneFlip
109 ) const
110 {
111  // Second or more usage of face. Add.
112  if (nei == -1 || own < nei)
113  {
114  meshMod.addFace
115  (
116  f, // modified face
117  own, // owner
118  nei, // neighbour
119  masterPointID, // master point
120  masterEdgeID, // master edge
121  masterFaceID, // master face
122  false, // face flip
123  patchi, // patch for face
124  zoneI, // zone for face
125  zoneFlip // face flip in zone
126  );
127  }
128  else
129  {
130  meshMod.addFace
131  (
132  f.reverseFace(), // modified face
133  nei, // owner
134  own, // neighbour
135  masterPointID, // master point
136  masterEdgeID, // master edge
137  masterFaceID, // master face
138  true, // face flip
139  patchi, // patch for face
140  zoneI, // zone for face
141  !zoneFlip // face flip in zone
142  );
143  }
144 }
145 
146 
147 // Work out triangle index given the starting vertex in the face
148 Foam::label Foam::tetDecomposer::triIndex(const label facei, const label fp)
149 const
150 {
151  const face& f = mesh_.faces()[facei];
152  const label fp0 = mesh_.tetBasePtIs()[facei];
153 
154  // Work out triangle index on this face
155  label thisTriI;
156  if (fp == fp0)
157  {
158  thisTriI = 0;
159  }
160  else if (fp == f.rcIndex(fp0))
161  {
162  thisTriI = f.size()-3;
163  }
164  else
165  {
166  thisTriI = (fp-fp0-1) % (f.size()-2);
167  }
168  return thisTriI;
169 }
170 
171 
172 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
173 
175 :
176  mesh_(mesh)
177 {}
178 
179 
180 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
181 
183 (
184  const decompositionType decomposeType,
185  polyTopoChange& meshMod
186 )
187 {
188  cellToPoint_.setSize(mesh_.nCells());
189  forAll(mesh_.cellCentres(), celli)
190  {
191  // Any point on the cell
192  label masterPointi = mesh_.faces()[mesh_.cells()[celli][0]][0];
193 
194  cellToPoint_[celli] = meshMod.addPoint
195  (
196  mesh_.cellCentres()[celli],
197  masterPointi,
198  -1,
199  true
200  );
201  }
202 
203 
204  // Add face centre points
205  if (decomposeType == FACE_CENTRE_TRIS)
206  {
207  faceToPoint_.setSize(mesh_.nFaces());
208  forAll(mesh_.faceCentres(), facei)
209  {
210  // Any point on the face
211  const label masterPointi = mesh_.faces()[facei][0];
212 
213  faceToPoint_[facei] = meshMod.addPoint
214  (
215  mesh_.faceCentres()[facei],
216  masterPointi,
217  -1,
218  true
219  );
220  }
221  }
222 
223 
224  // Per face, per point (faceCentre) or triangle (faceDiag) the added cell
225  faceOwnerCells_.setSize(mesh_.nFaces());
226  faceNeighbourCells_.setSize(mesh_.nFaces());
227 
228  if (decomposeType == FACE_CENTRE_TRIS)
229  {
230  forAll(faceOwnerCells_, facei)
231  {
232  const face& f = mesh_.faces()[facei];
233  faceOwnerCells_[facei].setSize(f.size(), -1);
234  faceNeighbourCells_[facei].setSize(f.size(), -1);
235  }
236  }
237  else
238  {
239  // Force construction of diagonal decomposition
240  (void)mesh_.tetBasePtIs();
241 
242  forAll(faceOwnerCells_, facei)
243  {
244  const face& f = mesh_.faces()[facei];
245  faceOwnerCells_[facei].setSize(f.size()-2, -1);
246  faceNeighbourCells_[facei].setSize(f.size()-2, -1);
247  }
248  }
249 
250 
251  forAll(mesh_.cells(), celli)
252  {
253  const cell& cFaces = mesh_.cells()[celli];
254 
255  EdgeMap<label> edgeToFace(8*cFaces.size());
256 
257  forAll(cFaces, cFacei)
258  {
259  label facei = cFaces[cFacei];
260  const face& f = mesh_.faces()[facei];
261 
262  // Get reference to either owner or neighbour
263  labelList& added =
264  (
265  (mesh_.faceOwner()[facei] == celli)
266  ? faceOwnerCells_[facei]
267  : faceNeighbourCells_[facei]
268  );
269 
270  if (decomposeType == FACE_CENTRE_TRIS)
271  {
272  forAll(f, fp)
273  {
274  if (cFacei == 0 && fp == 0)
275  {
276  // Reuse cell itself
277  added[fp] = celli;
278  }
279  else
280  {
281  added[fp] = meshMod.addCell
282  (
283  -1, // masterPoint
284  -1, // masterEdge
285  -1, // masterFace
286  celli, // masterCell
287  mesh_.cellZones().whichZone(celli)
288  );
289  }
290  }
291  }
292  else
293  {
294  for (label triI = 0; triI < f.size()-2; triI++)
295  {
296  if (cFacei == 0 && triI == 0)
297  {
298  // Reuse cell itself
299  added[triI] = celli;
300  }
301  else
302  {
303  added[triI] = meshMod.addCell
304  (
305  -1, // masterPoint
306  -1, // masterEdge
307  -1, // masterFace
308  celli, // masterCell
309  mesh_.cellZones().whichZone(celli)
310  );
311  }
312  }
313  }
314  }
315  }
316 
317 
318 
319  // Add triangle faces
320  face triangle(3);
321 
322  forAll(mesh_.faces(), facei)
323  {
324  label own = mesh_.faceOwner()[facei];
325  const labelList& addedOwn = faceOwnerCells_[facei];
326  const labelList& addedNei = faceNeighbourCells_[facei];
327  const face& f = mesh_.faces()[facei];
328 
329  label patchi = -1;
330  if (facei >= mesh_.nInternalFaces())
331  {
332  patchi = mesh_.boundaryMesh().whichPatch(facei);
333  }
334 
335  label zoneI = mesh_.faceZones().whichZone(facei);
336  bool zoneFlip = false;
337  if (zoneI != -1)
338  {
339  const faceZone& fz = mesh_.faceZones()[zoneI];
340  zoneFlip = fz.flipMap()[fz.whichFace(facei)];
341  }
342 
343 
344  if (decomposeType == FACE_CENTRE_TRIS)
345  {
346  forAll(f, fp)
347  {
348  // 1. Front triangle (decomposition of face itself)
349  // (between owner and neighbour cell)
350  {
351  triangle[0] = f[fp];
352  triangle[1] = f[f.fcIndex(fp)];
353  triangle[2] = faceToPoint_[facei];
354 
355  if (fp == 0)
356  {
357  modifyFace
358  (
359  meshMod,
360  triangle,
361  facei,
362  addedOwn[fp],
363  addedNei[fp],
364  patchi,
365  zoneI,
366  zoneFlip
367  );
368  }
369  else
370  {
371  addFace
372  (
373  meshMod,
374  triangle,
375  addedOwn[fp],
376  addedNei[fp],
377  -1, // point
378  -1, // edge
379  facei, // face
380  patchi,
381  zoneI,
382  zoneFlip
383  );
384  }
385  }
386 
387 
388  // 2. Within owner cell - to cell centre
389  {
390  label newOwn = addedOwn[f.rcIndex(fp)];
391  label newNei = addedOwn[fp];
392 
393  triangle[0] = f[fp];
394  triangle[1] = cellToPoint_[own];
395  triangle[2] = faceToPoint_[facei];
396 
397  addFace
398  (
399  meshMod,
400  triangle,
401  newOwn,
402  newNei,
403  f[fp], // point
404  -1, // edge
405  -1, // face
406  -1, // patchi
407  zoneI,
408  zoneFlip
409  );
410  }
411  // 2b. Within neighbour cell - to cell centre
412  if (facei < mesh_.nInternalFaces())
413  {
414  label newOwn = addedNei[f.rcIndex(fp)];
415  label newNei = addedNei[fp];
416 
417  triangle[0] = f[fp];
418  triangle[1] = faceToPoint_[facei];
419  triangle[2] = cellToPoint_[mesh_.faceNeighbour()[facei]];
420 
421  addFace
422  (
423  meshMod,
424  triangle,
425  newOwn,
426  newNei,
427  f[fp], // point
428  -1, // edge
429  -1, // face
430  -1, // patchi
431  zoneI,
432  zoneFlip
433  );
434  }
435  }
436  }
437  else
438  {
439  label fp0 = mesh_.tetBasePtIs()[facei];
440  label fp = f.fcIndex(fp0);
441 
442  for (label triI = 0; triI < f.size()-2; triI++)
443  {
444  label nextTri = triI+1;
445  if (nextTri >= f.size()-2)
446  {
447  nextTri -= f.size()-2;
448  }
449  label nextFp = f.fcIndex(fp);
450 
451 
452  // Triangle triI consisting of f[fp0], f[fp], f[nextFp]
453 
454 
455  // 1. Front triangle (decomposition of face itself)
456  // (between owner and neighbour cell)
457  {
458  triangle[0] = f[fp0];
459  triangle[1] = f[fp];
460  triangle[2] = f[nextFp];
461 
462  if (triI == 0)
463  {
464  modifyFace
465  (
466  meshMod,
467  triangle,
468  facei,
469  addedOwn[triI],
470  addedNei[triI],
471  patchi,
472  zoneI,
473  zoneFlip
474  );
475  }
476  else
477  {
478  addFace
479  (
480  meshMod,
481  triangle,
482  addedOwn[triI],
483  addedNei[triI],
484  -1, // point
485  -1, // edge
486  facei, // face
487  patchi,
488  zoneI,
489  zoneFlip
490  );
491  }
492  }
493 
494 
495  // 2. Within owner cell - diagonal to cell centre
496  if (triI < f.size()-3)
497  {
498  label newOwn = addedOwn[triI];
499  label newNei = addedOwn[nextTri];
500 
501  triangle[0] = f[fp0];
502  triangle[1] = f[nextFp];
503  triangle[2] = cellToPoint_[own];
504 
505  addFace
506  (
507  meshMod,
508  triangle,
509  newOwn,
510  newNei,
511  f[fp], // point
512  -1, // edge
513  -1, // face
514  -1, // patchi
515  zoneI,
516  zoneFlip
517  );
518 
519  // 2b. Within neighbour cell - to cell centre
520  if (facei < mesh_.nInternalFaces())
521  {
522  label newOwn = addedNei[triI];
523  label newNei = addedNei[nextTri];
524 
525  triangle[0] = f[nextFp];
526  triangle[1] = f[fp0];
527  triangle[2] =
528  cellToPoint_[mesh_.faceNeighbour()[facei]];
529 
530  addFace
531  (
532  meshMod,
533  triangle,
534  newOwn,
535  newNei,
536  f[fp], // point
537  -1, // edge
538  -1, // face
539  -1, // patchi
540  zoneI,
541  zoneFlip
542  );
543  }
544  }
545 
546 
547  fp = nextFp;
548  }
549  }
550  }
551 
552 
553 
554  // Add triangles for all edges.
555  EdgeMap<label> edgeToFace;
556 
557  forAll(mesh_.cells(), celli)
558  {
559  const cell& cFaces = mesh_.cells()[celli];
560 
561  edgeToFace.clear();
562 
563  forAll(cFaces, cFacei)
564  {
565  label facei = cFaces[cFacei];
566 
567  label zoneI = mesh_.faceZones().whichZone(facei);
568  bool zoneFlip = false;
569  if (zoneI != -1)
570  {
571  const faceZone& fz = mesh_.faceZones()[zoneI];
572  zoneFlip = fz.flipMap()[fz.whichFace(facei)];
573  }
574 
575  const face& f = mesh_.faces()[facei];
576  // const labelList& fEdges = mesh_.faceEdges()[facei];
577  forAll(f, fp)
578  {
579  label p0 = f[fp];
580  label p1 = f[f.fcIndex(fp)];
581  const edge e(p0, p1);
582 
583  EdgeMap<label>::const_iterator edgeFnd = edgeToFace.find(e);
584  if (edgeFnd == edgeToFace.end())
585  {
586  edgeToFace.insert(e, facei);
587  }
588  else
589  {
590  // Found the other face on the edge.
591  label otherFacei = edgeFnd();
592  const face& otherF = mesh_.faces()[otherFacei];
593 
594  // Found the other face on the edge. Note that since
595  // we are looping in the same order the tets added for
596  // otherFacei will be before those of facei
597 
598  label otherFp = findIndex(otherF, p0);
599  if (otherF.nextLabel(otherFp) == p1)
600  {
601  // ok. otherFp is first vertex of edge.
602  }
603  else if (otherF.prevLabel(otherFp) == p1)
604  {
605  otherFp = otherF.rcIndex(otherFp);
606  }
607  else
608  {
610  << "problem." << abort(FatalError);
611  }
612 
613 
614  // Triangle from edge to cell centre
615  if (mesh_.faceOwner()[facei] == celli)
616  {
617  triangle[0] = p0;
618  triangle[1] = p1;
619  triangle[2] = cellToPoint_[celli];
620  }
621  else
622  {
623  triangle[0] = p1;
624  triangle[1] = p0;
625  triangle[2] = cellToPoint_[celli];
626  }
627 
628  // Determine tets on either side
629  label thisTet, otherTet;
630 
631  if (decomposeType == FACE_CENTRE_TRIS)
632  {
633  if (mesh_.faceOwner()[facei] == celli)
634  {
635  thisTet = faceOwnerCells_[facei][fp];
636  }
637  else
638  {
639  thisTet = faceNeighbourCells_[facei][fp];
640  }
641 
642  if (mesh_.faceOwner()[otherFacei] == celli)
643  {
644  otherTet = faceOwnerCells_[otherFacei][otherFp];
645  }
646  else
647  {
648  otherTet =
649  faceNeighbourCells_[otherFacei][otherFp];
650  }
651  }
652  else
653  {
654  label thisTriI = triIndex(facei, fp);
655  if (mesh_.faceOwner()[facei] == celli)
656  {
657  thisTet = faceOwnerCells_[facei][thisTriI];
658  }
659  else
660  {
661  thisTet = faceNeighbourCells_[facei][thisTriI];
662  }
663 
664  label otherTriI = triIndex(otherFacei, otherFp);
665  if (mesh_.faceOwner()[otherFacei] == celli)
666  {
667  otherTet = faceOwnerCells_[otherFacei][otherTriI];
668  }
669  else
670  {
671  otherTet =
672  faceNeighbourCells_[otherFacei][otherTriI];
673  }
674  }
675 
676 
677  addFace
678  (
679  meshMod,
680  triangle,
681  otherTet,
682  thisTet,
683  -1, // masterPoint
684  -1, // fEdges[fp], // masterEdge
685  facei, // masterFace
686  -1, // patchi
687  zoneI,
688  zoneFlip
689  );
690  }
691  }
692  }
693  }
694 }
695 
696 
698 {
699  inplaceRenumber(map.reversePointMap(), cellToPoint_);
700  inplaceRenumber(map.reversePointMap(), faceToPoint_);
701 
702  forAll(faceOwnerCells_, facei)
703  {
704  inplaceRenumber(map.reverseCellMap(), faceOwnerCells_[facei]);
705  }
706  forAll(faceNeighbourCells_, facei)
707  {
708  inplaceRenumber(map.reverseCellMap(), faceNeighbourCells_[facei]);
709  }
710 }
711 
712 
713 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A triangle primitive used to calculate face areas and swept volumes.
Definition: triangle.H:57
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:112
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const meshCellZones & cellZones() const
Return cell zones.
Definition: polyMesh.H:501
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
label nInternalFaces() const
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
Definition: UListI.H:65
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1243
const labelIOList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:928
label nFaces() const
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:249
void setRefinement(const decompositionType decomposeType, polyTopoChange &meshMod)
Insert all changes into meshMod to convert the polyMesh into.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
label nCells() const
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:100
const cellList & cells() const
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:307
fvMesh & mesh
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:51
label addCell(const label masterPointID, const label masterEdgeID, const label masterFaceID, const label masterCellID, const label zoneID)
Add cell. Return new cell label.
label prevLabel(const label i) const
Previous vertex on face.
Definition: faceI.H:106
bool insert(const edge &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
iterator find(const edge &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
face reverseFace() const
Return face with reverse direction.
Definition: face.C:256
label addFace(const face &f, const label own, const label nei, const label masterPointID, const label masterEdgeID, const label masterFaceID, const bool flipFaceFlux, const label patchID, const label zoneID, const bool zoneFlip)
Add face to cells. Return new face label.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
void clear()
Clear all entries from table.
Definition: HashTable.C:468
void topoChange(const polyTopoChangeMap &)
Force recalculation of locally stored data on topological change.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1237
const labelList & reversePointMap() const
Reverse point map.
label addPoint(const point &, const label masterPointID, const label zoneID, const bool inCell)
Add point. Return new point label.
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1224
tetDecomposer(const polyMesh &)
Construct from mesh.
const vectorField & cellCentres() const
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Map from edge (expressed as its endpoints) to value.
Definition: EdgeMap.H:47
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
const labelList & reverseCellMap() const
Reverse cell map.
const vectorField & faceCentres() const
void setSize(const label)
Reset size of List.
Definition: List.C:281
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: MeshZones.C:221
void modifyFace(const face &f, const label facei, const label own, const label nei, const bool flipFaceFlux, const label patchID, const label zoneID, const bool zoneFlip)
Modify vertices or cell of face.
static const NamedEnum< decompositionType, 2 > decompositionTypeNames
Definition: tetDecomposer.H:70
label patchi
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
const meshFaceZones & faceZones() const
Return face zones.
Definition: polyMesh.H:495
Direct mesh changes based on v1.3 polyTopoChange syntax.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:65
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Namespace for OpenFOAM.