cellShapeControlMesh.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) 2016 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
31#include "pointIOField.H"
32#include "scalarIOField.H"
33#include "triadIOField.H"
34#include "tetPointRef.H"
35#include "plane.H"
36#include "transform.H"
37#include "meshTools.H"
38
39// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40
41namespace Foam
42{
43defineTypeNameAndDebug(cellShapeControlMesh, 0);
44}
45
47
48
49// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
50
51//Foam::tensor Foam::cellShapeControlMesh::requiredAlignment
52//(
53// const Foam::point& pt,
54// const searchableSurfaces& allGeometry,
55// const conformationSurfaces& geometryToConformTo
56//) const
57//{
58// pointIndexHit surfHit;
59// label hitSurface;
60//
61// geometryToConformTo.findSurfaceNearest
62// (
63// pt,
64// sqr(GREAT),
65// surfHit,
66// hitSurface
67// );
68//
69// if (!surfHit.hit())
70// {
71// FatalErrorInFunction
72// << "findSurfaceNearest did not find a hit across the surfaces."
73// << exit(FatalError) << endl;
74// }
75//
76// // Primary alignment
77//
78// vectorField norm(1);
79//
80// allGeometry[hitSurface].getNormal
81// (
82// List<pointIndexHit>(1, surfHit),
83// norm
84// );
85//
86// const vector np = norm[0];
87//
88// // Generate equally spaced 'spokes' in a circle normal to the
89// // direction from the vertex to the closest point on the surface
90// // and look for a secondary intersection.
91//
92// const vector d = surfHit.hitPoint() - pt;
93//
94// const tensor Rp = rotationTensor(vector(0,0,1), np);
95//
96// const label s = 36;//foamyHexMeshControls().alignmentSearchSpokes();
97//
98// scalar closestSpokeHitDistance = GREAT;
99//
100// pointIndexHit closestSpokeHit;
101//
102// label closestSpokeSurface = -1;
103//
104// const scalar spanMag = geometryToConformTo.globalBounds().mag();
105//
106// for (label i = 0; i < s; i++)
107// {
108// vector spoke
109// (
110// Foam::cos(i*constant::mathematical::twoPi/s),
111// Foam::sin(i*constant::mathematical::twoPi/s),
112// 0
113// );
114//
115// spoke *= spanMag;
116//
117// spoke = Rp & spoke;
118//
119// pointIndexHit spokeHit;
120//
121// label spokeSurface = -1;
122//
123// // internal spoke
124//
125// geometryToConformTo.findSurfaceNearestIntersection
126// (
127// pt,
128// pt + spoke,
129// spokeHit,
130// spokeSurface
131// );
132//
133// if (spokeHit.hit())
134// {
135// scalar spokeHitDistance = mag
136// (
137// spokeHit.hitPoint() - pt
138// );
139//
140// if (spokeHitDistance < closestSpokeHitDistance)
141// {
142// closestSpokeHit = spokeHit;
143// closestSpokeSurface = spokeSurface;
144// closestSpokeHitDistance = spokeHitDistance;
145// }
146// }
147//
148// //external spoke
149//
150// Foam::point mirrorPt = pt + 2*d;
151//
152// geometryToConformTo.findSurfaceNearestIntersection
153// (
154// mirrorPt,
155// mirrorPt + spoke,
156// spokeHit,
157// spokeSurface
158// );
159//
160// if (spokeHit.hit())
161// {
162// scalar spokeHitDistance = mag
163// (
164// spokeHit.hitPoint() - mirrorPt
165// );
166//
167// if (spokeHitDistance < closestSpokeHitDistance)
168// {
169// closestSpokeHit = spokeHit;
170// closestSpokeSurface = spokeSurface;
171// closestSpokeHitDistance = spokeHitDistance;
172// }
173// }
174// }
175//
176// if (closestSpokeSurface == -1)
177// {
183//
184// return I;
185// }
186//
187// // Auxiliary alignment generated by spoke intersection normal.
188//
189// allGeometry[closestSpokeSurface].getNormal
190// (
191// List<pointIndexHit>(1, closestSpokeHit),
192// norm
193// );
194//
195// const vector& na = norm[0];
196//
197// // Secondary alignment
198// vector ns = np ^ na;
199//
200// if (mag(ns) < SMALL)
201// {
202// FatalErrorInFunction
203// << "Parallel normals detected in spoke search." << nl
204// << "point: " << pt << nl
205// << "closest surface point: " << surfHit.hitPoint() << nl
206// << "closest spoke hit: " << closestSpokeHit.hitPoint() << nl
207// << "np: " << surfHit.hitPoint() + np << nl
208// << "ns: " << closestSpokeHit.hitPoint() + na << nl
209// << exit(FatalError);
210// }
211//
212// ns /= mag(ns);
213//
214// tensor Rs = rotationTensor((Rp & vector(0,1,0)), ns);
215//
216// return (Rs & Rp);
217//}
218
219
221{
222 label nRemoved = 0;
223 for
224 (
225 CellSizeDelaunay::Finite_vertices_iterator vit =
226 finite_vertices_begin();
227 vit != finite_vertices_end();
228 ++vit
229 )
230 {
231 std::list<Vertex_handle> verts;
232 adjacent_vertices(vit, std::back_inserter(verts));
233
234 bool removePt = true;
235 for
236 (
237 std::list<Vertex_handle>::iterator aVit = verts.begin();
238 aVit != verts.end();
239 ++aVit
240 )
241 {
242 Vertex_handle avh = *aVit;
243
244 scalar diff =
245 mag(avh->targetCellSize() - vit->targetCellSize())
246 /max(vit->targetCellSize(), 1e-6);
247
248 if (diff > 0.05)
249 {
250 removePt = false;
251 }
252 }
253
254 if (removePt)
255 {
256 remove(vit);
257 nRemoved++;
258 }
259 }
260
261 return nRemoved;
262}
263
264
266{
267 tmp<pointField> tcellCentres(new pointField(number_of_finite_cells()));
268 pointField& cellCentres = tcellCentres.ref();
269
270 label count = 0;
271 for
272 (
273 CellSizeDelaunay::Finite_cells_iterator c = finite_cells_begin();
274 c != finite_cells_end();
275 ++c
276 )
277 {
278 if (c->hasFarPoint())
279 {
280 continue;
281 }
282
283 scalarList bary;
285
286 const Foam::point centre = topoint
287 (
288 CGAL::centroid<baseK>
289 (
290 c->vertex(0)->point(),
291 c->vertex(1)->point(),
292 c->vertex(2)->point(),
293 c->vertex(3)->point()
294 )
295 );
296
297 cellCentres[count++] = centre;
298 }
299
300 cellCentres.resize(count);
301
302 return tcellCentres;
303}
304
305
307{
308 OFstream str
309 (
310 "refinementTriangulation_"
312 + ".obj"
313 );
314
315 label count = 0;
316
317 Info<< "Write refinementTriangulation" << endl;
318
319 for
320 (
321 CellSizeDelaunay::Finite_edges_iterator e = finite_edges_begin();
322 e != finite_edges_end();
323 ++e
324 )
325 {
326 Cell_handle c = e->first;
327 Vertex_handle vA = c->vertex(e->second);
328 Vertex_handle vB = c->vertex(e->third);
329
330 // Don't write far edges
331 if (vA->farPoint() || vB->farPoint())
332 {
333 continue;
334 }
335
336 // Don't write unowned edges
337 if (vA->referred() && vB->referred())
338 {
339 continue;
340 }
341
342 pointFromPoint p1 = topoint(vA->point());
343 pointFromPoint p2 = topoint(vB->point());
344
345 meshTools::writeOBJ(str, p1, p2, count);
346 }
347
348 if (is_valid())
349 {
350 Info<< " Triangulation is valid" << endl;
351 }
352 else
353 {
355 << "Triangulation is not valid"
356 << abort(FatalError);
357 }
358}
359
360
361// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
362
364:
365 DistributedDelaunayMesh<CellSizeDelaunay>
366 (
367 runTime,
368 meshSubDir
369 ),
370 runTime_(runTime)
371{
372 if (this->vertexCount())
373 {
374 fvMesh mesh
375 (
376 IOobject
377 (
378 meshSubDir,
379 runTime.timeName(),
380 runTime,
381 IOobject::READ_IF_PRESENT,
382 IOobject::NO_WRITE
383 )
384 );
385
386 if (mesh.nPoints() == this->vertexCount())
387 {
388 IOobject io
389 (
390 "sizes",
391 runTime.timeName(),
392 meshSubDir,
393 runTime,
394 IOobject::MUST_READ,
395 IOobject::NO_WRITE,
396 false
397 );
398
399 if (io.typeHeaderOk<pointScalarField>(true))
400 {
401 pointScalarField sizes(io, pointMesh::New(mesh));
402
403 triadIOField alignments
404 (
405 IOobject
406 (
407 "alignments",
408 mesh.time().timeName(),
409 meshSubDir,
410 mesh.time(),
411 IOobject::MUST_READ,
412 IOobject::NO_WRITE,
413 false
414 )
415 );
416
417 if (alignments.size() == this->vertexCount())
418 {
419 for
420 (
421 Finite_vertices_iterator vit = finite_vertices_begin();
422 vit != finite_vertices_end();
423 ++vit
424 )
425 {
426 vit->targetCellSize() = sizes[vit->index()];
427 vit->alignment() = alignments[vit->index()];
428 }
429 }
430 else
431 {
433 << "Cell alignments point field " << alignments.size()
434 << " is not the same size as the number of vertices"
435 << " in the mesh " << this->vertexCount()
436 << abort(FatalError);
437 }
438 }
439 }
440 }
441}
442
443
444// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
445
447(
448 const Foam::point& pt,
449 barycentric& bary,
450 Cell_handle& ch
451) const
452{
453 // Use the previous cell handle as a hint on where to start searching
454 // Giving a hint causes strange errors...
455 ch = locate(toPoint(pt));
456
457 if (dimension() > 2 && !is_infinite(ch))
458 {
459 oldCellHandle_ = ch;
460
461 tetPointRef tet
462 (
463 topoint(ch->vertex(0)->point()),
464 topoint(ch->vertex(1)->point()),
465 topoint(ch->vertex(2)->point()),
466 topoint(ch->vertex(3)->point())
467 );
468
469 bary = tet.pointToBarycentric(pt);
470 }
471}
472
473
475{
476 DynamicList<Foam::point> pts(number_of_vertices());
477
478 for
479 (
480 Finite_vertices_iterator vit = finite_vertices_begin();
481 vit != finite_vertices_end();
482 ++vit
483 )
484 {
485 if (vit->real())
486 {
487 pts.append(topoint(vit->point()));
488 }
489 }
490
491 boundBox bb(pts);
492
493 return bb;
494}
495
496
498(
499 const backgroundMeshDecomposition& decomposition
500)
501{
502 DynamicList<Foam::point> points(number_of_vertices());
503 DynamicList<scalar> sizes(number_of_vertices());
504 DynamicList<tensor> alignments(number_of_vertices());
505
506 DynamicList<Vb> farPts(8);
507
508 for
509 (
510 Finite_vertices_iterator vit = finite_vertices_begin();
511 vit != finite_vertices_end();
512 ++vit
513 )
514 {
515 if (vit->real())
516 {
517 points.append(topoint(vit->point()));
518 sizes.append(vit->targetCellSize());
519 alignments.append(vit->alignment());
520 }
521 else if (vit->farPoint())
522 {
523 farPts.append
524 (
525 Vb
526 (
527 vit->point(),
528 -1,
529 Vb::vtFar,
531 )
532 );
533
534 farPts.last().targetCellSize() = vit->targetCellSize();
535 farPts.last().alignment() = vit->alignment();
536 }
537 }
538
539 autoPtr<mapDistribute> mapDist =
540 DistributedDelaunayMesh<CellSizeDelaunay>::distribute
541 (
542 decomposition,
543 points
544 );
545
546 mapDist().distribute(sizes);
547 mapDist().distribute(alignments);
548
549 // Reset the entire tessellation
551
552
553 // Internal points have to be inserted first
554 DynamicList<Vb> verticesToInsert(points.size());
555
556
557 forAll(farPts, ptI)
558 {
559 verticesToInsert.append(farPts[ptI]);
560 }
561
562
563 forAll(points, pI)
564 {
565 verticesToInsert.append
566 (
567 Vb
568 (
569 toPoint(points[pI]),
570 -1,
573 )
574 );
575
576 verticesToInsert.last().targetCellSize() = sizes[pI];
577 verticesToInsert.last().alignment() = alignments[pI];
578 }
579
580 Info<< nl << " Inserting distributed background tessellation..." << endl;
581
582 this->rangeInsertWithInfo
583 (
584 verticesToInsert.begin(),
585 verticesToInsert.end(),
586 true
587 );
588
589 sync(decomposition.procBounds());
590
591 Info<< " Total number of vertices after redistribution "
592 << returnReduce(label(number_of_vertices()), sumOp<label>()) << endl;
593}
594
595
597{
598 tensorField alignmentsTmp(number_of_vertices(), Zero);
599
600 label count = 0;
601 for
602 (
603 Finite_vertices_iterator vit = finite_vertices_begin();
604 vit != finite_vertices_end();
605 ++vit
606 )
607 {
608 alignmentsTmp[count++] = vit->alignment();
609 }
610
611 return alignmentsTmp;
612}
613
614
616{
617 Info<< "Writing " << meshSubDir << endl;
618
619 // Reindex the cells
620 label cellCount = 0;
621 for
622 (
623 Finite_cells_iterator cit = finite_cells_begin();
624 cit != finite_cells_end();
625 ++cit
626 )
627 {
628 if (!cit->hasFarPoint() && !is_infinite(cit))
629 {
630 cit->cellIndex() = cellCount++;
631 }
632 }
633
634 labelPairLookup vertexMap;
635 labelList cellMap;
636
637 autoPtr<polyMesh> meshPtr = DelaunayMesh<CellSizeDelaunay>::createMesh
638 (
639 meshSubDir,
640 vertexMap,
641 cellMap
642 );
643 const polyMesh& mesh = meshPtr();
644
645 pointScalarField sizes
646 (
647 IOobject
648 (
649 "sizes",
650 mesh.time().timeName(),
651 meshSubDir,
652 mesh.time(),
655 ),
658 );
659
660 triadIOField alignments
661 (
662 IOobject
663 (
664 "alignments",
665 mesh.time().timeName(),
666 meshSubDir,
667 mesh.time(),
670 ),
671 sizes.size()
672 );
673
674 // Write alignments
675// OFstream str(runTime_.path()/"alignments.obj");
676
677 for
678 (
679 Finite_vertices_iterator vit = finite_vertices_begin();
680 vit != finite_vertices_end();
681 ++vit
682 )
683 {
684 if (!vit->farPoint())
685 {
686 // Populate sizes
687 sizes[vertexMap[labelPair(vit->index(), vit->procIndex())]] =
688 vit->targetCellSize();
689
690 alignments[vertexMap[labelPair(vit->index(), vit->procIndex())]] =
691 vit->alignment();
692
693// // Write alignments
694// const tensor& alignment = vit->alignment();
695// pointFromPoint pt = topoint(vit->point());
696//
697// if
698// (
699// alignment.x() == triad::unset[0]
700// || alignment.y() == triad::unset[0]
701// || alignment.z() == triad::unset[0]
702// )
703// {
704// Info<< "Bad alignment = " << vit->info();
705//
706// vit->alignment() = tensor::I;
707//
708// Info<< "New alignment = " << vit->info();
709//
710// continue;
711// }
712//
713// meshTools::writeOBJ(str, pt, alignment.x() + pt);
714// meshTools::writeOBJ(str, pt, alignment.y() + pt);
715// meshTools::writeOBJ(str, pt, alignment.z() + pt);
716 }
717 }
718
719 mesh.write();
720 sizes.write();
721 alignments.write();
722}
723
724
726(
727 const autoPtr<backgroundMeshDecomposition>& decomposition
728) const
729{
730 // Loop over all the tets and estimate the cell count in each one
731
732 scalar cellCount = 0;
733
734 for
735 (
736 Finite_cells_iterator cit = finite_cells_begin();
737 cit != finite_cells_end();
738 ++cit
739 )
740 {
741 if (!cit->hasFarPoint() && !is_infinite(cit))
742 {
743 // TODO: Check if tet centre is on the processor..
744 CGAL::Tetrahedron_3<baseK> tet
745 (
746 cit->vertex(0)->point(),
747 cit->vertex(1)->point(),
748 cit->vertex(2)->point(),
749 cit->vertex(3)->point()
750 );
751
752 pointFromPoint centre = topoint(CGAL::centroid(tet));
753
754 if
755 (
757 && !decomposition().positionOnThisProcessor(centre)
758 )
759 {
760 continue;
761 }
762
763 scalar volume = CGAL::to_double(tet.volume());
764
765 scalar averagedPointCellSize = 0;
766 //scalar averagedPointCellSize = 1;
767
768 // Get an average volume by averaging the cell size of the vertices
769 for (label vI = 0; vI < 4; ++vI)
770 {
771 averagedPointCellSize += cit->vertex(vI)->targetCellSize();
772 //averagedPointCellSize *= cit->vertex(vI)->targetCellSize();
773 }
774
775 averagedPointCellSize /= 4;
776 //averagedPointCellSize = ::sqrt(averagedPointCellSize);
777
778// if (averagedPointCellSize < SMALL)
779// {
780// Pout<< "Volume = " << volume << endl;
781//
782// for (label vI = 0; vI < 4; ++vI)
783// {
784// Pout<< "Point " << vI
785// << ", point = " << topoint(cit->vertex(vI)->point())
786// << ", size = " << cit->vertex(vI)->targetCellSize()
787// << endl;
788// }
789// }
790
791 cellCount += volume/pow(averagedPointCellSize, 3);
792 }
793 }
794
795 return cellCount;
796}
797
798
799// ************************************************************************* //
CGAL::Delaunay_triangulation_3< K, Tds, FastLocator > CellSizeDelaunay
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Definition: indexedVertex.H:88
void reset()
Clear the entire triangulation.
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
static word timeName(const scalar t, const int precision=precision_)
Definition: Time.C:780
T & first()
Return the first element of the list.
Definition: UListI.H:202
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
boundBox bounds() const
tmp< pointField > cellCentres() const
Get the centres of all the tets.
static word meshSubDir
Return the mesh sub-directory name (usually "cellShapeControlMesh")
label estimateCellCount(const autoPtr< backgroundMeshDecomposition > &decomposition) const
CellSizeDelaunay::Cell_handle Cell_handle
void barycentricCoords(const Foam::point &pt, barycentric &bary, Cell_handle &ch) const
Calculate and return the barycentric coordinates for.
CellSizeDelaunay::Vertex_handle Vertex_handle
tensorField dumpAlignments() const
void distribute(const backgroundMeshDecomposition &decomposition)
virtual bool write(const bool valid=true) const
Write mesh using IO settings from time.
Definition: fvMesh.C:1079
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
int myProcNo() const noexcept
Return processor number.
A class for managing temporary objects.
Definition: tmp.H:65
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
dynamicFvMesh & mesh
engineTime & runTime
Foam::autoPtr< Foam::dynamicFvMesh > meshPtr
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
Definition: BitOps.H:78
const dimensionedScalar c
Speed of light in a vacuum.
const wordList volume
Standard volume field types (scalar, vector, tensor, etc)
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
Namespace for OpenFOAM.
IOField< triad > triadIOField
triadField with IO.
Definition: triadIOField.H:44
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:57
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
tetrahedron< point, const point & > tetPointRef
A tetrahedron using referred points.
Definition: tetPointRef.H:47
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
List< label > labelList
A List of labels.
Definition: List.H:66
PointFrompoint toPoint(const Foam::point &p)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
pointFromPoint topoint(const Point &P)
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:64
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
HashTable< label, labelPair, Foam::Hash< labelPair > > labelPairLookup
This is a Map of a labelPair to a label. Used for e.g. for face1, face2 to shared edge....
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
Definition: triad.C:378
errorManip< error > abort(error &err)
Definition: errorManip.H:144
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
3D tensor transformation operations.