"README.md" did not exist on "e42af28f7dbca878a6edb6b2ff9bdf7ab6ff1937"
Newer
Older
/*---------------------------------------------------------------------------*\
========= |
Franjo
committed
\\ / A nd | Copyright held by the origina author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
cfMesh is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
cfMesh is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with cfMesh. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "polyMeshGenChecks.H"
#include "polyMeshGenAddressing.H"
#include "pyramidPointFaceRef.H"
#include "tetrahedron.H"
# ifdef USE_OMP
# endif
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace polyMeshGenChecks
{
bool checkClosedBoundary(const polyMeshGen& mesh, const bool report)
{
// Loop through all boundary faces and sum up the face area vectors.
// For a closed boundary, this should be zero in all vector components
vector sumClosed(vector::zero);
scalar sumMagClosedBoundary = 0;
const vectorField& areas = mesh.addressingData().faceAreas();
for(label faceI = mesh.nInternalFaces();faceI<areas.size();++faceI)
{
sumClosed += areas[faceI];
sumMagClosedBoundary += mag(areas[faceI]);
}
// Check the openness in the maximum direction (this is APPROXIMATE!)
scalar maxOpen = max
(
sumClosed.component(vector::X),
max
(
sumClosed.component(vector::Y),
sumClosed.component(vector::Z)
)
);
reduce(sumClosed, sumOp<vector>());
reduce(maxOpen, maxOp<scalar>());
if( maxOpen > SMALL*max(1.0, sumMagClosedBoundary) )
{
SeriousErrorIn
(
"bool checkClosedBoundary(const polyMeshGen&, const bool report)"
) << "Possible hole in boundary description" << endl;
Info<< "Boundary openness in x-direction = "
<< sumClosed.component(vector::X) << endl;
Info<< "Boundary openness in y-direction = "
<< sumClosed.component(vector::Y) << endl;
Info<< "Boundary openness in z-direction = "
<< sumClosed.component(vector::Z) << endl;
return true;
}
else
{
if( report )
{
Info<< "Boundary openness in x-direction = "
<< sumClosed.component(vector::X) << endl;
Info<< "Boundary openness in y-direction = "
<< sumClosed.component(vector::Y) << endl;
Info<< "Boundary openness in z-direction = "
<< sumClosed.component(vector::Z) << endl;
Info<< "Boundary closed (OK)." << endl;
}
return false;
}
}
bool checkClosedCells
(
const polyMeshGen& mesh,
const bool report,
const scalar aspectWarn,
labelHashSet* setPtr
)
{
// Check that all cells labels are valid
const cellListPMG& cells = mesh.cells();
const label nFaces = mesh.faces().size();
# ifdef USE_OMP
# pragma omp parallel for schedule(guided) reduction(+ : nErrorClosed)
# endif
forAll(cells, cI)
{
const cell& curCell = cells[cI];
if( min(curCell) < 0 || max(curCell) > nFaces )
{
WarningIn
(
"bool checkClosedCells("
"const polyMeshGen&, const bool, const scalar, labelHashSet*)"
) << "Cell " << cI << " contains face labels out of range: "
<< curCell << " Max face index = " << nFaces << endl;
if( setPtr )
{
# ifdef USE_OMP
# endif
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
setPtr->insert(cI);
}
++nErrorClosed;
}
}
if( nErrorClosed > 0 )
{
SeriousErrorIn
(
"bool checkClosedCells("
"const polyMeshGen&, const bool, const scalar, labelHashSet*)"
) << nErrorClosed << " cells with invalid face labels found"
<< endl;
return true;
}
// Loop through cell faces and sum up the face area vectors for each cell.
// This should be zero in all vector components
vectorField sumClosed(cells.size(), vector::zero);
scalarField sumMagClosed(cells.size(), 0.0);
const labelList& own = mesh.owner();
const labelList& nei = mesh.neighbour();
const vectorField& areas = mesh.addressingData().faceAreas();
forAll(own, faceI)
{
// Add to owner
sumClosed[own[faceI]] += areas[faceI];
sumMagClosed[own[faceI]] += mag(areas[faceI]);
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
// Subtract from neighbour
sumClosed[nei[faceI]] -= areas[faceI];
sumMagClosed[nei[faceI]] += mag(areas[faceI]);
}
label nOpen = 0;
scalar maxOpenCell = 0;
label nAspect = 0;
scalar maxAspectRatio = 0;
const scalarField& vols = mesh.addressingData().cellVolumes();
// Check the sums
forAll(sumClosed, cellI)
{
scalar maxOpen = max
(
sumClosed[cellI].component(vector::X),
max
(
sumClosed[cellI].component(vector::Y),
sumClosed[cellI].component(vector::Z)
)
);
maxOpenCell = max(maxOpenCell, maxOpen);
if( maxOpen > SMALL*max(1.0, sumMagClosed[cellI]) )
{
if( report )
{
Pout<< "Cell " << cellI << " is not closed. "
<< "Face area vectors sum up to " << mag(sumClosed[cellI])
<< " directionwise " << sumClosed[cellI] << " or "
<< mag(sumClosed[cellI])
/(mag(sumMagClosed[cellI]) + VSMALL)
<< " of the area of all faces of the cell. " << endl
<< " Area magnitudes sum up to "
<< sumMagClosed[cellI] << endl;
}
if( setPtr )
{
setPtr->insert(cellI);
}
++nOpen;
}
scalar aspectRatio =
1.0/6.0*sumMagClosed[cellI]/pow(vols[cellI], 2.0/3.0);
maxAspectRatio = max(maxAspectRatio, aspectRatio);
if( aspectRatio > aspectWarn )
{
if( report )
{
Pout<< "High aspect ratio for cell " << cellI
<< ": " << aspectRatio << endl;
}
++nAspect;
}
}
reduce(nOpen, sumOp<label>());
reduce(maxOpenCell, maxOp<scalar>());
reduce(nAspect, sumOp<label>());
reduce(maxAspectRatio, maxOp<scalar>());
if( nOpen > 0 )
{
SeriousErrorIn
(
"bool checkClosedCells("
"const polyMeshGen&, const bool, const scalar, labelHashSet*)"
) << nOpen << " open cells found. Max cell openness: "
<< maxOpenCell << endl;
return true;
}
if( nAspect > 0 )
{
SeriousErrorIn
(
"bool checkClosedCells("
"const polyMeshGen&, const bool, const scalar, labelHashSet*)"
) << nAspect << " high aspect ratio cells found. "
<< "Max aspect ratio: " << maxAspectRatio
<< endl;
return true;
}
if( report )
{
Info<< "Max cell openness = " << maxOpenCell
<< " Max aspect ratio = " << maxAspectRatio
<< ". All cells OK.\n" << endl;
}
return false;
}
bool checkCellVolumes
(
const polyMeshGen& mesh,
const bool report,
labelHashSet* setPtr
)
{
const scalarField& vols = mesh.addressingData().cellVolumes();
scalar minVolume = GREAT;
scalar maxVolume = -GREAT;
label nNegVolCells = 0;
forAll(vols, cellI)
{
if( vols[cellI] < VSMALL )
{
if( report )
SeriousErrorIn
(
"bool checkCellVolumes("
"const polyMeshGen&, const bool, labelHashSet*)"
) << "Zero or negative cell volume detected for cell "
<< cellI << ". Volume = " << vols[cellI] << endl;
if( setPtr )
setPtr->insert(cellI);
++nNegVolCells;
}
minVolume = min(minVolume, vols[cellI]);
maxVolume = max(maxVolume, vols[cellI]);
}
reduce(minVolume, minOp<scalar>());
reduce(maxVolume, maxOp<scalar>());
reduce(nNegVolCells, sumOp<label>());
if( minVolume < VSMALL )
{
SeriousErrorIn
(
"bool checkCellVolumes("
"const polyMeshGen&, const bool, labelHashSet*)"
) << "Zero or negative cell volume detected. "
<< "Minimum negative volume: "
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
<< minVolume << ".\nNumber of negative volume cells: "
<< nNegVolCells << ". This mesh is invalid"
<< endl;
return true;
}
else
{
if( report )
{
Info<< "Min volume = " << minVolume
<< ". Max volume = " << maxVolume
<< ". Total volume = " << sum(vols)
<< ". Cell volumes OK.\n" << endl;
}
return false;
}
}
bool checkFaceAreas
(
const polyMeshGen& mesh,
const bool report,
const scalar minFaceArea,
labelHashSet* setPtr,
const boolList* changedFacePtr
)
{
const vectorField& faceAreas = mesh.addressingData().faceAreas();
const labelList& own = mesh.owner();
const labelList& nei = mesh.neighbour();
scalar minArea = VGREAT;
scalar maxArea = -VGREAT;
# ifdef USE_OMP
# endif
# ifdef USE_OMP
# endif
{
if( changedFacePtr && !changedFacePtr->operator[](faceI) )
continue;
const scalar magFaceArea = mag(faceAreas[faceI]);
if( magFaceArea < minFaceArea )
{
if( report )
{
if( nei[faceI] != -1 )
{
Pout<< "Zero or negative face area detected for "
<< "internal face " << faceI << " between cells "
<< own[faceI] << " and " << nei[faceI]
<< ". Face area magnitude = "
}
else
{
Pout<< "Zero or negative face area detected for "
<< "boundary face " << faceI << " next to cell "
<< own[faceI] << ". Face area magnitude = "
# ifdef USE_OMP
# endif
localMinArea = Foam::min(localMinArea, magFaceArea);
localMaxArea = Foam::max(localMaxArea, magFaceArea);
# ifdef USE_OMP
# endif
{
minArea = Foam::min(minArea, localMinArea);
maxArea = Foam::max(maxArea, localMaxArea);
}
}
reduce(minArea, minOp<scalar>());
reduce(maxArea, maxOp<scalar>());
if( minArea < VSMALL )
{
SeriousErrorIn
(
"bool checkFaceAreas("
"const polyMeshGen&, const bool, const scalar,"
" labelHashSet*, const boolList*)"
) << "Zero or negative face area detected. Minimum negative area: "
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
<< minArea << ". This mesh is invalid"
<< endl;
return true;
}
else
{
if( report )
{
Info<< "Minumum face area = " << minArea
<< ". Maximum face area = " << maxArea
<< ". Face area magnitudes OK.\n" << endl;
}
return false;
}
}
bool checkCellPartTetrahedra
(
const polyMeshGen& mesh,
const bool report,
const scalar minPartTet,
labelHashSet* setPtr,
const boolList* changedFacePtr
)
{
const pointFieldPMG& points = mesh.points();
const faceListPMG& faces = mesh.faces();
const labelList& owner = mesh.owner();
const labelList& neighbour = mesh.neighbour();
const vectorField& fCentres = mesh.addressingData().faceCentres();
const vectorField& cCentres = mesh.addressingData().cellCentres();
label nNegVolCells = 0;
# ifdef USE_OMP
# pragma omp parallel for if( owner.size() > 100 ) \
schedule(guided) reduction(+ : nNegVolCells)
# endif
forAll(owner, faceI)
{
if( changedFacePtr && !(*changedFacePtr)[faceI] )
continue;
forAll(f, eI)
{
const tetrahedron<point, point> tetOwn
(
fCentres[faceI],
points[f.nextLabel(eI)],
points[f[eI]],
cCentres[owner[faceI]]
);
if( tetOwn.mag() < minPartTet )
{
if( report )
{
# ifdef USE_OMP
# endif
Pout<< "Zero or negative cell volume detected for cell "
<< owner[faceI] << "." << endl;
}
const tetrahedron<point, point> tetNei
(
fCentres[faceI],
points[f[eI]],
points[f.nextLabel(eI)],
cCentres[neighbour[faceI]]
);
if( tetNei.mag() < minPartTet )
{
if( report )
{
# ifdef USE_OMP
# endif
Pout<< "Zero or negative cell volume detected for cell "
<< neighbour[faceI] << "." << endl;
}
# ifdef USE_OMP
# endif
if( setPtr )
{
//- ensure that faces are selected at both sides
const PtrList<processorBoundaryPatch>& procBnd = mesh.procBoundaries();
forAll(procBnd, patchI)
{
const label start = procBnd[patchI].patchStart();
const label size = procBnd[patchI].patchSize();
franjo_j@hotmail.com
committed
labelLongList sendData;
for(label faceI=0;faceI<size;++faceI)
{
if( setPtr->found(faceI+start) )
sendData.append(faceI);
}
OPstream toOtherProc
(
Pstream::blocking,
procBnd[patchI].neiProcNo(),
sendData.byteSize()
);
IPstream fromOtherProc
(
Pstream::blocking,
procBnd[patchI].neiProcNo()
);
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
const label start = procBnd[patchI].patchStart();
forAll(receivedData, i)
setPtr->insert(start+receivedData[i]);
}
}
reduce(nNegVolCells, sumOp<label>());
if( nNegVolCells != 0 )
{
WarningIn
(
"bool checkCellPartTetrahedra("
"const polyMeshGen&, const bool, const scalar,"
" labelHashSet*, const boolList*)"
) << nNegVolCells << " zero or negative part tetrahedra detected."
<< endl;
return true;
}
else
{
if( report )
Info<< "Part cell tetrahedra OK.\n" << endl;
return false;
}
}
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
void checkFaceDotProduct
(
const polyMeshGen& mesh,
scalarField& faceDotProduct,
const boolList* changedFacePtr
)
{
//- for all internal faces check theat the d dot S product is positive
const vectorField& centres = mesh.addressingData().cellCentres();
const vectorField& areas = mesh.addressingData().faceAreas();
const labelList& own = mesh.owner();
const labelList& nei = mesh.neighbour();
const label nInternalFaces = mesh.nInternalFaces();
faceDotProduct.setSize(own.size());
faceDotProduct = 1.0;
# ifdef USE_OMP
# pragma omp parallel if( nInternalFaces > 1000 )
# endif
{
# ifdef USE_OMP
# pragma omp for schedule(dynamic, 100)
# endif
for(label faceI=0;faceI<nInternalFaces;++faceI)
{
if( changedFacePtr && !(*changedFacePtr)[faceI] )
continue;
const vector d = centres[nei[faceI]] - centres[own[faceI]];
const vector& s = areas[faceI];
faceDotProduct[faceI] = (d & s)/(mag(d)*mag(s) + VSMALL);
}
}
if( Pstream::parRun() )
{
const PtrList<processorBoundaryPatch>& procBoundaries =
mesh.procBoundaries();
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
vectorField cCentres(procBoundaries[patchI].patchSize());
forAll(cCentres, faceI)
cCentres[faceI] = centres[own[start+faceI]];
OPstream toOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo(),
cCentres.byteSize()
);
toOtherProc << cCentres;
}
forAll(procBoundaries, patchI)
{
vectorField otherCentres;
IPstream fromOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo()
);
fromOtherProc >> otherCentres;
//- calculate skewness at processor faces
const label start = procBoundaries[patchI].patchStart();
# ifdef USE_OMP
# pragma omp parallel
# endif
{
# ifdef USE_OMP
# pragma omp for schedule(dynamic, 100)
# endif
forAll(otherCentres, fI)
{
const label faceI = start + fI;
if( changedFacePtr && !(*changedFacePtr)[faceI] )
continue;
const point& cOwn = centres[own[faceI]];
const point& cNei = otherCentres[fI];
const vector d = cNei - cOwn;
const vector& s = areas[faceI];
faceDotProduct[faceI] = (d & s)/(mag(d)*mag(s) + VSMALL);
}
}
}
}
}
bool checkFaceDotProduct
(
const polyMeshGen& mesh,
const bool report,
const scalar nonOrthWarn,
labelHashSet* setPtr,
const boolList* changedFacePtr
)
{
//- for all internal faces check theat the d dot S product is positive
scalarField faceDotProduct;
checkFaceDotProduct(mesh, faceDotProduct, changedFacePtr);
const labelList& own = mesh.owner();
const labelList& nei = mesh.neighbour();
const label nInternalFaces = mesh.nInternalFaces();
// Severe nonorthogonality threshold
const scalar severeNonorthogonalityThreshold =
::cos(nonOrthWarn/180.0*M_PI);
scalar minDDotS = VGREAT;
scalar sumDDotS = 0.0;
label severeNonOrth = 0;
label errorNonOrth = 0;
label counter = 0;
# ifdef USE_OMP
# pragma omp parallel if( nInternalFaces > 1000 ) \
reduction(+ : severeNonOrth, errorNonOrth, sumDDotS)
# endif
# ifdef USE_OMP
# endif
if( changedFacePtr && !(*changedFacePtr)[faceI] )
continue;
const scalar dDotS = faceDotProduct[faceI];
if( dDotS < severeNonorthogonalityThreshold )
{
if( dDotS > SMALL )
{
if( report )
{
// Severe non-orthogonality but mesh still OK
# ifdef USE_OMP
# endif
Pout<< "Severe non-orthogonality for face " << faceI
<< " between cells " << own[faceI]
<< " and " << nei[faceI]
<< ": Angle = "
<< ::acos(dDotS)/M_PI*180.0
# ifdef USE_OMP
# endif
# ifdef USE_OMP
# endif
localMinDDotS = Foam::min(dDotS, localMinDDotS);
sumDDotS += dDotS;
}
# ifdef USE_OMP
# endif
const PtrList<processorBoundaryPatch>& procBoundaries =
const label start = procBoundaries[0].patchStart();
# ifdef USE_OMP
# pragma omp parallel \
reduction(+ : severeNonOrth, errorNonOrth, sumDDotS, counter)
# endif
scalar localMinDDotS(VGREAT);
# ifdef USE_OMP
# pragma omp for schedule(dynamic, 100)
# endif
for(label faceI=start;faceI<own.size();++faceI)
if( changedFacePtr && !(*changedFacePtr)[start+faceI] )
continue;
const scalar dDotS = faceDotProduct[faceI];
if( dDotS < severeNonorthogonalityThreshold )
{
if( dDotS > SMALL )
if( report )
// Severe non-orthogonality but mesh still OK
# ifdef USE_OMP
# pragma omp critical
# endif
const scalar angle
(
Foam::acos(dDotS) /
M_PI * 180.0
);
Pout<< "Severe non-orthogonality for face "
<< start+faceI
<< ": Angle = "
<< angle
<< " deg." << endl;
if( setPtr )
{
# ifdef USE_OMP
# pragma omp critical
# endif
setPtr->insert(start+faceI);
++severeNonOrth;
else
{
++errorNonOrth;
if( setPtr )
{
# ifdef USE_OMP
# pragma omp critical
# endif
setPtr->insert(start+faceI);
}
}
localMinDDotS = Foam::min(dDotS, localMinDDotS);
sumDDotS += 0.5 * dDotS;
++counter;
# ifdef USE_OMP
# pragma omp critical
# endif
minDDotS = Foam::min(minDDotS, localMinDDotS);
reduce(minDDotS, minOp<scalar>());
reduce(sumDDotS, sumOp<scalar>());
reduce(severeNonOrth, sumOp<label>());
reduce(errorNonOrth, sumOp<label>());
reduce(counter, sumOp<label>());
// Only report if there are some internal faces
if( counter > 0 )
{
if( minDDotS < severeNonorthogonalityThreshold )
{
Info<< "Number of non-orthogonality errors: " << errorNonOrth
<< ". Number of severely non-orthogonal faces: "
<< severeNonOrth << "." << endl;
}
}
if( report )
{
if( counter > 0 )
{
Info<< "Mesh non-orthogonality Max: "
<< ::acos(minDDotS)/M_PI*180.0
::acos(sumDDotS/counter)/M_PI*180.0
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
<< endl;
}
}
if( errorNonOrth > 0 )
{
WarningIn
(
"checkFaceDotProduct("
"const polyMeshGen&, const bool, const scalar,"
" labelHashSet*, const boolList*)"
) << "Error in non-orthogonality detected" << endl;
return true;
}
else
{
if( report )
Info<< "Non-orthogonality check OK.\n" << endl;
return false;
}
}
bool checkFacePyramids
(
const polyMeshGen& mesh,
const bool report,
const scalar minPyrVol,
labelHashSet* setPtr,
const boolList* changedFacePtr
)
{
// check whether face area vector points to the cell with higher label
const vectorField& ctrs = mesh.addressingData().cellCentres();
const labelList& owner = mesh.owner();
const labelList& neighbour = mesh.neighbour();
const faceListPMG& faces = mesh.faces();
const pointFieldPMG& points = mesh.points();
label nErrorPyrs = 0;
# ifdef USE_OMP
# pragma omp parallel for schedule(guided) reduction(+ : nErrorPyrs)
# endif
if( changedFacePtr && !(*changedFacePtr)[faceI] )
continue;
// Create the owner pyramid - it will have negative volume
const scalar pyrVol = pyramidPointFaceRef
(
faces[faceI],
ctrs[owner[faceI]]
).mag(points);