cellSizeAndAlignmentGrid.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-2016 OpenFOAM Foundation
9 Copyright (C) 2020 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
27Application
28 Test-distributedDelaunayMesh
29
30Description
31
32\*---------------------------------------------------------------------------*/
33
35
36#include "indexedVertex.H"
37#include "indexedCell.H"
38
39#include "argList.H"
40#include "Time.H"
43#include "searchableSurfaces.H"
45#include "PrintTable.H"
46#include "Random.H"
47#include "boundBox.H"
48#include "point.H"
50#include "triadField.H"
51#include "scalarIOField.H"
52#include "pointIOField.H"
53#include "triadIOField.H"
54
55using namespace Foam;
56
57// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
58// Main program:
59
60template<class Triangulation, class Type>
61Foam::tmp<Foam::Field<Type>> filterFarPoints
62(
63 const Triangulation& mesh,
64 const Field<Type>& field
65)
66{
67 tmp<Field<Type>> tNewField(new Field<Type>(field.size()));
68 Field<Type>& newField = tNewField.ref();
69
70 label added = 0;
71 label count = 0;
72
73 for
74 (
76 mesh.finite_vertices_begin();
77 vit != mesh.finite_vertices_end();
78 ++vit
79 )
80 {
81 if (vit->real())
82 {
83 newField[added++] = field[count];
84 }
85
86 count++;
87 }
88
89 newField.resize(added);
90
91 return tNewField;
92}
93
94
95template<class T>
97(
98 const T& mesh,
99 labelListList& pointPoints
100)
101{
102 pointPoints.setSize(mesh.vertexCount());
103
104 globalIndex globalIndexing(mesh.vertexCount());
105
106 for
107 (
108 typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
109 vit != mesh.finite_vertices_end();
110 ++vit
111 )
112 {
113 if (!vit->real())
114 {
115 continue;
116 }
117
118 std::list<typename T::Vertex_handle> adjVerts;
119 mesh.finite_adjacent_vertices(vit, std::back_inserter(adjVerts));
120
121 DynamicList<label> indices(adjVerts.size());
122
123 for
124 (
125 typename std::list<typename T::Vertex_handle>::const_iterator
126 adjVertI = adjVerts.begin();
127 adjVertI != adjVerts.end();
128 ++adjVertI
129 )
130 {
131 typename T::Vertex_handle vh = *adjVertI;
132
133 if (!vh->farPoint())
134 {
135 indices.append
136 (
137 globalIndexing.toGlobal(vh->procIndex(), vh->index())
138 );
139 }
140 }
141
142 pointPoints[vit->index()].transfer(indices);
143 }
144
145 List<Map<label>> compactMap;
146
148 (
149 globalIndexing,
150 pointPoints,
151 compactMap
152 );
153}
154
155
156template<class T>
157Foam::tmp<Foam::triadField> buildAlignmentField(const T& mesh)
158{
159 tmp<triadField> tAlignments
160 (
161 new triadField(mesh.vertexCount(), triad::unset)
162 );
163 triadField& alignments = tAlignments.ref();
164
165 for
166 (
167 typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
168 vit != mesh.finite_vertices_end();
169 ++vit
170 )
171 {
172 if (!vit->real())
173 {
174 continue;
175 }
176
177 alignments[vit->index()] = vit->alignment();
178 }
179
180 return tAlignments;
181}
182
183
184template<class T>
185Foam::tmp<Foam::pointField> buildPointField(const T& mesh)
186{
187 tmp<pointField> tPoints
188 (
189 new pointField(mesh.vertexCount(), point(GREAT, GREAT, GREAT))
190 );
191 pointField& points = tPoints.ref();
192
193 for
194 (
195 typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
196 vit != mesh.finite_vertices_end();
197 ++vit
198 )
199 {
200 if (!vit->real())
201 {
202 continue;
203 }
204
205 points[vit->index()] = topoint(vit->point());
206 }
207
208 return tPoints;
209}
210
211
212void refine
213(
215 const conformationSurfaces& geometryToConformTo,
216 const label maxRefinementIterations,
217 const scalar defaultCellSize
218)
219{
220 for (label iter = 0; iter < maxRefinementIterations; ++iter)
221 {
222 DynamicList<point> ptsToInsert;
223
224 for
225 (
226 CellSizeDelaunay::Finite_cells_iterator cit =
227 mesh.finite_cells_begin();
228 cit != mesh.finite_cells_end();
229 ++cit
230 )
231 {
232 const point newPoint =
233 topoint
234 (
235 CGAL::centroid
236 (
237 cit->vertex(0)->point(),
238 cit->vertex(1)->point(),
239 cit->vertex(2)->point(),
240 cit->vertex(3)->point()
241 )
242 );
243
244 if (geometryToConformTo.inside(newPoint))
245 {
246 ptsToInsert.append(newPoint);
247 }
248 }
249
250 Info<< " Adding " << returnReduce(ptsToInsert.size(), sumOp<label>())
251 << endl;
252
253 forAll(ptsToInsert, ptI)
254 {
255 mesh.insert
256 (
257 ptsToInsert[ptI],
258 defaultCellSize,
259 triad::unset,
261 );
262 }
263 }
264}
265
266
267int main(int argc, char *argv[])
268{
269 #include "setRootCase.H"
270 #include "createTime.H"
271
272 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
273
274 label maxRefinementIterations = 2;
275 label maxSmoothingIterations = 200;
276 scalar minResidual = 0;
277 scalar defaultCellSize = 0.001;
278 scalar nearFeatDistSqrCoeff = 1e-8;
279
280
281 // Need to decouple vertex and cell type from this class?
282 // Vertex must have:
283 // + index
284 // + procIndex
285 // - type should be optional
287
288 IOdictionary foamyHexMeshDict
289 (
291 (
292 "foamyHexMeshDict",
293 runTime.system(),
294 runTime,
295 IOobject::MUST_READ,
296 IOobject::NO_WRITE
297 )
298 );
299
300 Random rndGen(64293*Pstream::myProcNo());
301
302 searchableSurfaces allGeometry
303 (
305 (
306 "cvSearchableSurfaces",
307 runTime.constant(),
308 "triSurface",
309 runTime,
310 IOobject::MUST_READ,
311 IOobject::NO_WRITE
312 ),
313 foamyHexMeshDict.subDict("geometry"),
314 foamyHexMeshDict.getOrDefault("singleRegionName", true)
315 );
316
317 conformationSurfaces geometryToConformTo
318 (
319 runTime,
320 rndGen,
321 allGeometry,
322 foamyHexMeshDict.subDict("surfaceConformation")
323 );
324
326 if (Pstream::parRun())
327 {
328 bMesh.set
329 (
331 (
332 runTime,
333 rndGen,
334 geometryToConformTo,
335 foamyHexMeshDict.subDict("backgroundMeshDecomposition")
336 )
337 );
338 }
339
340 // Nice to have IO for the delaunay mesh
341 // IO depend on vertex type.
342 //
343 // Define a delaunay mesh as:
344 // + list of points of the triangulation
345 // + optionally a list of cells
346
347 Info<< nl << "Loop over surfaces" << endl;
348
349 forAll(geometryToConformTo.surfaces(), sI)
350 {
351 const label surfI = geometryToConformTo.surfaces()[sI];
352
354 geometryToConformTo.geometry()[surfI];
355
356 Info<< nl << "Inserting points from surface " << surface.name()
357 << " (" << surface.type() << ")" << endl;
358
359 const tmp<pointField> tpoints(surface.points());
360 const pointField& points = tpoints();
361
362 Info<< " Number of points = " << points.size() << endl;
363
364 forAll(points, pI)
365 {
366 // Is the point in the extendedFeatureEdgeMesh? If so get the
367 // point normal, otherwise get the surface normal from
368 // searchableSurface
369
370 pointIndexHit info;
371 label infoFeature;
372 geometryToConformTo.findFeaturePointNearest
373 (
374 points[pI],
375 nearFeatDistSqrCoeff,
376 info,
377 infoFeature
378 );
379
380
381 autoPtr<triad> pointAlignment;
382
383 if (info.hit())
384 {
385 const extendedFeatureEdgeMesh& features =
386 geometryToConformTo.features()[infoFeature];
387
388 vectorField norms = features.featurePointNormals(info.index());
389
390 // Create a triad from these norms.
391 pointAlignment.set(new triad());
392 forAll(norms, nI)
393 {
394 pointAlignment() += norms[nI];
395 }
396
397 pointAlignment().normalize();
398 pointAlignment().orthogonalize();
399 }
400 else
401 {
402 geometryToConformTo.findEdgeNearest
403 (
404 points[pI],
405 nearFeatDistSqrCoeff,
406 info,
407 infoFeature
408 );
409
410 if (info.hit())
411 {
412 const extendedFeatureEdgeMesh& features =
413 geometryToConformTo.features()[infoFeature];
414
415 vectorField norms = features.edgeNormals(info.index());
416
417 // Create a triad from these norms.
418 pointAlignment.set(new triad());
419 forAll(norms, nI)
420 {
421 pointAlignment() += norms[nI];
422 }
423
424 pointAlignment().normalize();
425 pointAlignment().orthogonalize();
426 }
427 else
428 {
429 pointField ptField(1, points[pI]);
430 scalarField distField(1, nearFeatDistSqrCoeff);
431 List<pointIndexHit> infoList(1, pointIndexHit());
432
433 surface.findNearest(ptField, distField, infoList);
434
435 vectorField normals(1);
436 surface.getNormal(infoList, normals);
437
438 pointAlignment.set(new triad(normals[0]));
439 }
440 }
441
442 if (Pstream::parRun())
443 {
444 if (bMesh().positionOnThisProcessor(points[pI]))
445 {
446 CellSizeDelaunay::Vertex_handle vh = mesh.insert
447 (
448 points[pI],
449 defaultCellSize,
450 pointAlignment(),
452 );
453 }
454 }
455 else
456 {
457 CellSizeDelaunay::Vertex_handle vh = mesh.insert
458 (
459 points[pI],
460 defaultCellSize,
461 pointAlignment(),
463 );
464 }
465 }
466 }
467
468
469 // Refine the mesh
470 refine
471 (
472 mesh,
473 geometryToConformTo,
474 maxRefinementIterations,
475 defaultCellSize
476 );
477
478
479 if (Pstream::parRun())
480 {
481 mesh.distribute(bMesh);
482 }
483
484
485 labelListList pointPoints;
486 autoPtr<mapDistribute> meshDistributor = buildMap(mesh, pointPoints);
487
488
489 triadField alignments(buildAlignmentField(mesh));
490 pointField points(buildPointField(mesh));
491
492 mesh.printInfo(Info);
493
494
495 // Setup the sizes and alignments on each point
496 triadField fixedAlignments(mesh.vertexCount(), triad::unset);
497
498 for
499 (
500 CellSizeDelaunay::Finite_vertices_iterator vit =
501 mesh.finite_vertices_begin();
502 vit != mesh.finite_vertices_end();
503 ++vit
504 )
505 {
506 if (vit->nearBoundary())
507 {
508 fixedAlignments[vit->index()] = vit->alignment();
509 }
510 }
511
512 Info<< nl << "Smoothing alignments" << endl;
513
514 for (label iter = 0; iter < maxSmoothingIterations; iter++)
515 {
516 Info<< "Iteration " << iter;
517
518 meshDistributor().distribute(points);
519 meshDistributor().distribute(alignments);
520
521 scalar residual = 0;
522
523 triadField triadAv(alignments.size(), triad::unset);
524
525 forAll(pointPoints, pI)
526 {
527 const labelList& pPoints = pointPoints[pI];
528
529 if (pPoints.empty())
530 {
531 continue;
532 }
533
534 const triad& oldTriad = alignments[pI];
535 triad& newTriad = triadAv[pI];
536
537 // Enforce the boundary conditions
538 const triad& fixedAlignment = fixedAlignments[pI];
539
540 forAll(pPoints, adjPointi)
541 {
542 const label adjPointIndex = pPoints[adjPointi];
543
544 scalar dist = mag(points[pI] - points[adjPointIndex]);
545
546// dist = max(dist, SMALL);
547
548 triad tmpTriad = alignments[adjPointIndex];
549
550 for (direction dir = 0; dir < 3; dir++)
551 {
552 if (tmpTriad.set(dir))
553 {
554 tmpTriad[dir] *= (1.0/dist);
555 }
556 }
557
558 newTriad += tmpTriad;
559 }
560
561 newTriad.normalize();
562 newTriad.orthogonalize();
563// newTriad = newTriad.sortxyz();
564
565 label nFixed = 0;
566
567 forAll(fixedAlignment, dirI)
568 {
569 if (fixedAlignment.set(dirI))
570 {
571 nFixed++;
572 }
573 }
574
575 if (nFixed == 1)
576 {
577 forAll(fixedAlignment, dirI)
578 {
579 if (fixedAlignment.set(dirI))
580 {
581 newTriad.align(fixedAlignment[dirI]);
582 }
583 }
584 }
585 else if (nFixed == 2)
586 {
587 forAll(fixedAlignment, dirI)
588 {
589 if (fixedAlignment.set(dirI))
590 {
591 newTriad[dirI] = fixedAlignment[dirI];
592 }
593 else
594 {
595 newTriad[dirI] = triad::unset[dirI];
596 }
597 }
598
599 newTriad.orthogonalize();
600 }
601 else if (nFixed == 3)
602 {
603 forAll(fixedAlignment, dirI)
604 {
605 if (fixedAlignment.set(dirI))
606 {
607 newTriad[dirI] = fixedAlignment[dirI];
608 }
609 }
610 }
611
612 for (direction dir = 0; dir < 3; ++dir)
613 {
614 if
615 (
616 newTriad.set(dir)
617 && oldTriad.set(dir)
618 //&& !fixedAlignment.set(dir)
619 )
620 {
621 scalar dotProd = (oldTriad[dir] & newTriad[dir]);
622 scalar diff = mag(dotProd) - 1.0;
623
624 residual += mag(diff);
625 }
626 }
627 }
628
629 forAll(alignments, pI)
630 {
631 alignments[pI] = triadAv[pI].sortxyz();
632 }
633
634 reduce(residual, sumOp<scalar>());
635
636 Info<< ", Residual = " << residual << endl;
637
638 if (residual <= minResidual)
639 {
640 break;
641 }
642 }
643
644
645 // Write alignments to a .obj file
646 OFstream str(runTime.path()/"alignments.obj");
647
648 forAll(alignments, pI)
649 {
650 const triad& tri = alignments[pI];
651
652 if (tri.set())
653 {
654 forAll(tri, dirI)
655 {
656 meshTools::writeOBJ(str, points[pI], tri[dirI] + points[pI]);
657 }
658 }
659 }
660
661
662 // Remove the far points
663 pointIOField pointsIO
664 (
666 (
667 "points",
668 runTime.constant(),
669 runTime,
670 IOobject::NO_READ,
671 IOobject::AUTO_WRITE
672 ),
673 filterFarPoints(mesh, points)
674 );
675
676 scalarField sizes(points.size(), defaultCellSize);
677 scalarIOField sizesIO
678 (
680 (
681 "sizes",
682 runTime.constant(),
683 runTime,
684 IOobject::NO_READ,
685 IOobject::AUTO_WRITE
686 ),
687 filterFarPoints(mesh, sizes)
688 );
689
690 triadIOField alignmentsIO
691 (
693 (
694 "alignments",
695 runTime.constant(),
696 runTime,
697 IOobject::NO_READ,
698 IOobject::AUTO_WRITE
699 ),
700 filterFarPoints(mesh, alignments)
701 );
702
703 pointsIO.write();
704 sizesIO.write();
705 alignmentsIO.write();
706
707 Info<< nl;
708 runTime.printExecutionTime(Info);
709
710 Info<< "\nEnd\n" << endl;
711
712 return 0;
713}
714
715
716// ************************************************************************* //
reduce(hasMovingMesh, orOp< bool >())
Triangulation::Finite_vertices_iterator Finite_vertices_iterator
Definition: DelaunayMesh.H:73
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
Generic templated field type.
Definition: Field.H:82
A primitive field of type <T> with automated input and output.
Definition: IOField.H:58
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
void transfer(List< T > &list)
Definition: List.C:447
void setSize(const label n)
Alias for resize()
Definition: List.H:218
void resize(const label len)
Adjust allocated size of list.
Definition: ListI.H:139
Output to file stream, using an OSstream.
Definition: OFstream.H:57
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
label index() const noexcept
Return the hit index.
bool hit() const noexcept
Is there a hit?
A list of faces which address into the list of points.
Random number generator.
Definition: Random.H:60
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Definition: UListI.H:427
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
void set(T *p) noexcept
Deprecated(2018-02) Identical to reset().
Definition: autoPtr.H:288
Store a background polyMesh to use for the decomposition of space and queries for parallel conformalV...
const PtrList< extendedFeatureEdgeMesh > & features() const
Return the object holding the feature points and edges.
const searchableSurfaces & geometry() const
Return reference to the searchableSurfaces object containing all.
void findEdgeNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &edgeHit, label &featureHit) const
Find the nearest point on any feature edge.
void findFeaturePointNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &fpHit, label &featureHit) const
Find the nearest point on any feature edge.
const labelList & surfaces() const
Return the surface indices.
Field< bool > inside(const pointField &samplePts) const
Check if points are inside surfaces to conform to.
const labelListList & edgeNormals() const
Return the indices of the normals that are adjacent to the.
const labelListList & featurePointNormals() const
Return the indices of the normals that are adjacent to the.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
Container for searchableSurfaces. The collection is specified as a dictionary. For example,...
A class for managing temporary objects.
Definition: tmp.H:65
Representation of a 3D Cartesian coordinate system as a Vector of row vectors.
Definition: triad.H:67
void align(const vector &v)
Align this triad with the given vector v.
Definition: triad.C:240
bool set(const direction d) const
Is the vector in the direction d set.
Definition: triadI.H:70
void normalize()
Same as normalise.
Definition: triad.H:170
void orthogonalize()
Same as orthogonalise.
Definition: triad.H:167
const volScalarField & T
rDeltaTY field()
dynamicFvMesh & mesh
engineTime & runTime
const pointField & points
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
Definition: BitOps.H:78
const wordList surface
Standard surface field types (scalar, vector, tensor, etc)
Namespace for OpenFOAM.
pointFromPoint topoint(const Point &P)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
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
uint8_t direction
Definition: direction.H:56
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
Random rndGen
Definition: createFields.H:23