featurePointConformer.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) 2013-2016 OpenFOAM Foundation
9 Copyright (C) 2019-2022 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
30#include "cvControls.H"
33#include "cellShapeControl.H"
34#include "DelaunayMeshTools.H"
35#include "Circulator.H"
37#include "autoPtr.H"
38#include "mapDistribute.H"
39
40// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41
42namespace Foam
43{
44defineTypeNameAndDebug(featurePointConformer, 0);
45}
46
47const Foam::scalar Foam::featurePointConformer::tolParallel = 1e-3;
48
49
50// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
51
52Foam::vector Foam::featurePointConformer::sharedFaceNormal
53(
54 const extendedFeatureEdgeMesh& feMesh,
55 const label edgeI,
56 const label nextEdgeI
57) const
58{
59 const labelList& edgeInormals = feMesh.edgeNormals()[edgeI];
60 const labelList& nextEdgeInormals = feMesh.edgeNormals()[nextEdgeI];
61
62 const vector& A1 = feMesh.normals()[edgeInormals[0]];
63 const vector& A2 = feMesh.normals()[edgeInormals[1]];
64
65 const vector& B1 = feMesh.normals()[nextEdgeInormals[0]];
66 const vector& B2 = feMesh.normals()[nextEdgeInormals[1]];
67
68// Info<< " A1 = " << A1 << endl;
69// Info<< " A2 = " << A2 << endl;
70// Info<< " B1 = " << B1 << endl;
71// Info<< " B2 = " << B2 << endl;
72
73 const scalar A1B1 = mag((A1 & B1) - 1.0);
74 const scalar A1B2 = mag((A1 & B2) - 1.0);
75 const scalar A2B1 = mag((A2 & B1) - 1.0);
76 const scalar A2B2 = mag((A2 & B2) - 1.0);
77
78// Info<< " A1B1 = " << A1B1 << endl;
79// Info<< " A1B2 = " << A1B2 << endl;
80// Info<< " A2B1 = " << A2B1 << endl;
81// Info<< " A2B2 = " << A2B2 << endl;
82
83 if (A1B1 < A1B2 && A1B1 < A2B1 && A1B1 < A2B2)
84 {
85 return 0.5*(A1 + B1);
86 }
87 else if (A1B2 < A1B1 && A1B2 < A2B1 && A1B2 < A2B2)
88 {
89 return 0.5*(A1 + B2);
90 }
91 else if (A2B1 < A1B1 && A2B1 < A1B2 && A2B1 < A2B2)
92 {
93 return 0.5*(A2 + B1);
94 }
95 else
96 {
97 return 0.5*(A2 + B2);
98 }
99}
100
101
102Foam::label Foam::featurePointConformer::getSign
103(
105) const
106{
108 {
109 return -1;
110 }
111 else if (eStatus == extendedFeatureEdgeMesh::INTERNAL)
112 {
113 return 1;
114 }
115
116 return 0;
117}
118
119
120void Foam::featurePointConformer::addMasterAndSlavePoints
121(
122 const DynamicList<Foam::point>& masterPoints,
123 const DynamicList<Foam::indexedVertexEnum::vertexType>& masterPointsTypes,
124 const Map<DynamicList<autoPtr<plane>>>& masterPointReflections,
125 DynamicList<Vb>& pts,
126 const label ptI
127) const
128{
129 typedef DynamicList<autoPtr<plane>> planeDynList;
130 typedef Foam::indexedVertexEnum::vertexType vertexType;
131
132 forAll(masterPoints, pI)
133 {
134 // Append master to the list of points
135
136 const Foam::point& masterPt = masterPoints[pI];
137 const vertexType masterType = masterPointsTypes[pI];
138
139 pts.append
140 (
141 Vb
142 (
143 masterPt,
144 foamyHexMesh_.vertexCount() + pts.size(),
145 masterType,
147 )
148 );
149
150 const label masterIndex = pts.last().index();
151
152 //meshTools::writeOBJ(strMasters, masterPt);
153
154 const planeDynList& masterPointPlanes = masterPointReflections[pI];
155
156 forAll(masterPointPlanes, planeI)
157 {
158 // Reflect master points in the planes and insert the slave points
159
160 const plane& reflPlane = masterPointPlanes[planeI]();
161
162 const Foam::point slavePt = reflPlane.mirror(masterPt);
163
164 const vertexType slaveType =
165 (
166 masterType == Vb::vtInternalFeaturePoint
169 );
170
171 pts.append
172 (
173 Vb
174 (
175 slavePt,
176 foamyHexMesh_.vertexCount() + pts.size(),
177 slaveType,
179 )
180 );
181
182 ftPtPairs_.addPointPair
183 (
184 masterIndex,
185 pts.last().index()
186 );
187
188 //meshTools::writeOBJ(strSlaves, slavePt);
189 }
190 }
191}
192
193
194void Foam::featurePointConformer::createMasterAndSlavePoints
195(
196 const extendedFeatureEdgeMesh& feMesh,
197 const label ptI,
198 DynamicList<Vb>& pts
199) const
200{
201 typedef DynamicList<autoPtr<plane>> planeDynList;
202 typedef indexedVertexEnum::vertexType vertexType;
203 typedef extendedFeatureEdgeMesh::edgeStatus edgeStatus;
204
205 const Foam::point& featPt = feMesh.points()[ptI];
206
207 if
208 (
209 (
211 && !foamyHexMesh_.decomposition().positionOnThisProcessor(featPt)
212 )
213 || geometryToConformTo_.outside(featPt)
214 )
215 {
216 return;
217 }
218
219 const scalar ppDist = foamyHexMesh_.pointPairDistance(featPt);
220
221 // Maintain a list of master points and the planes to relect them in
222 DynamicList<Foam::point> masterPoints;
223 DynamicList<vertexType> masterPointsTypes;
224 Map<planeDynList> masterPointReflections;
225
226 const labelList& featPtEdges = feMesh.featurePointEdges()[ptI];
227
228 pointFeatureEdgesTypes pointEdgeTypes(feMesh, ptI);
229
230 const List<extendedFeatureEdgeMesh::edgeStatus> allEdStat =
231 pointEdgeTypes.calcPointFeatureEdgesTypes();
232
233// Info<< nl << featPt << " " << pointEdgeTypes;
234
235 ConstCirculator<labelList> circ(featPtEdges);
236
237 // Loop around the edges of the feature point
238 if (circ.size()) do
239 {
240// const edgeStatus eStatusPrev = feMesh.getEdgeStatus(circ.prev());
241 const edgeStatus eStatusCurr = feMesh.getEdgeStatus(circ.curr());
242// const edgeStatus eStatusNext = feMesh.getEdgeStatus(circ.next());
243
244// Info<< " Prev = "
245// << extendedFeatureEdgeMesh::edgeStatusNames_[eStatusPrev]
246// << " Curr = "
247// << extendedFeatureEdgeMesh::edgeStatusNames_[eStatusCurr]
250// << endl;
251
252 // Get the direction in which to move the point in relation to the
253 // feature point
254 label sign = getSign(eStatusCurr);
255
256 const vector n = sharedFaceNormal(feMesh, circ.curr(), circ.next());
257
258 const vector pointMotionDirection = sign*0.5*ppDist*n;
259
260// Info<< " Shared face normal = " << n << endl;
261// Info<< " Direction to move point = " << pointMotionDirection
262// << endl;
263
264 if (masterPoints.empty())
265 {
266 // Initialise with the first master point
267 Foam::point pt = featPt + pointMotionDirection;
268
269 planeDynList firstPlane;
270 firstPlane.append(autoPtr<plane>::New(featPt, n));
271
272 masterPoints.append(pt);
273
274 masterPointsTypes.append
275 (
276 sign == 1
279 );
280
281 masterPointReflections.insert
282 (
283 masterPoints.size() - 1,
284 std::move(firstPlane)
285 );
286 }
287// else if
288// (
289// eStatusPrev == extendedFeatureEdgeMesh::INTERNAL
290// && eStatusCurr == extendedFeatureEdgeMesh::EXTERNAL
291// )
292// {
293// // Insert a new master point.
294// Foam::point pt = featPt + pointMotionDirection;
295//
296// planeDynList firstPlane;
297// firstPlane.append(autoPtr<plane>::New(featPt, n));
298//
299// masterPoints.append(pt);
300//
301// masterPointsTypes.append
302// (
303// sign == 1
304// ? Vb::vtExternalFeaturePoint // true
305// : Vb::vtInternalFeaturePoint // false
306// );
307//
308// masterPointReflections.insert
309// (
310// masterPoints.size() - 1,
311// firstPlane
312// );
313// }
314// else if
315// (
316// eStatusPrev == extendedFeatureEdgeMesh::EXTERNAL
317// && eStatusCurr == extendedFeatureEdgeMesh::INTERNAL
318// )
319// {
320//
321// }
322 else
323 {
324 // Just add this face contribution to the latest master point
325
326 masterPoints.last() += pointMotionDirection;
327
328 masterPointReflections[masterPoints.size() - 1].append
329 (
330 autoPtr<plane>::New(featPt, n)
331 );
332 }
333
334 } while (circ.circulate(CirculatorBase::CLOCKWISE));
335
336 addMasterAndSlavePoints
337 (
338 masterPoints,
339 masterPointsTypes,
340 masterPointReflections,
341 pts,
342 ptI
343 );
344}
345
346
347void Foam::featurePointConformer::createMixedFeaturePoints
348(
349 DynamicList<Vb>& pts
350) const
351{
352 if (foamyHexMeshControls_.mixedFeaturePointPPDistanceCoeff() < 0)
353 {
354 Info<< nl
355 << "Skipping specialised handling for mixed feature points" << endl;
356
357 return;
358 }
359
360 const PtrList<extendedFeatureEdgeMesh>& feMeshes
361 (
362 geometryToConformTo_.features()
363 );
364
365 forAll(feMeshes, i)
366 {
367 const extendedFeatureEdgeMesh& feMesh = feMeshes[i];
368 const labelListList& pointsEdges = feMesh.pointEdges();
369 const pointField& points = feMesh.points();
370
371 for
372 (
373 label ptI = feMesh.mixedStart();
374 ptI < feMesh.nonFeatureStart();
375 ptI++
376 )
377 {
378 const Foam::point& featPt = points[ptI];
379
380 if
381 (
383 && !foamyHexMesh_.decomposition().positionOnThisProcessor(featPt))
384 {
385 continue;
386 }
387
388 const labelList& pEds = pointsEdges[ptI];
389
390 pointFeatureEdgesTypes pFEdgeTypes(feMesh, ptI);
391
392 const List<extendedFeatureEdgeMesh::edgeStatus> allEdStat =
393 pFEdgeTypes.calcPointFeatureEdgesTypes();
394
395 bool specialisedSuccess = false;
396
397 if (foamyHexMeshControls_.specialiseFeaturePoints())
398 {
399 specialisedSuccess =
400 createSpecialisedFeaturePoint
401 (
402 feMesh, pEds, pFEdgeTypes, allEdStat, ptI, pts
403 );
404 }
405
406 if (!specialisedSuccess && foamyHexMeshControls_.edgeAiming())
407 {
408 // Specialisations available for some mixed feature points. For
409 // non-specialised feature points, inserting mixed internal and
410 // external edge groups at feature point.
411
412 // Skipping unsupported mixed feature point types
413// bool skipEdge = false;
414//
415// forAll(pEds, e)
416// {
417// const label edgeI = pEds[e];
418//
419// const extendedFeatureEdgeMesh::edgeStatus edStatus
420// = feMesh.getEdgeStatus(edgeI);
421//
422// if
423// (
424// edStatus == extendedFeatureEdgeMesh::OPEN
425// || edStatus == extendedFeatureEdgeMesh::MULTIPLE
426// )
427// {
428// Info<< "Edge type " << edStatus
429// << " found for mixed feature point " << ptI
430// << ". Not supported."
431// << endl;
432//
433// skipEdge = true;
434// }
435// }
436//
437// if (skipEdge)
438// {
439// Info<< "Skipping point " << ptI << nl << endl;
440//
441// continue;
442// }
443
444// createFeaturePoints(feMesh, ptI, pts, types);
445
446 const Foam::point& pt = points[ptI];
447
448 const scalar edgeGroupDistance =
449 foamyHexMesh_.mixedFeaturePointDistance(pt);
450
451 forAll(pEds, e)
452 {
453 const label edgeI = pEds[e];
454
455 const Foam::point edgePt =
456 pt + edgeGroupDistance*feMesh.edgeDirection(edgeI, ptI);
457
458 const pointIndexHit edgeHit(true, edgePt, edgeI);
459
460 foamyHexMesh_.createEdgePointGroup(feMesh, edgeHit, pts);
461 }
462 }
463 }
464 }
465}
466
467
468void Foam::featurePointConformer::createFeaturePoints(DynamicList<Vb>& pts)
469{
470 const PtrList<extendedFeatureEdgeMesh>& feMeshes
471 (
472 geometryToConformTo_.features()
473 );
474
475 forAll(feMeshes, i)
476 {
477 const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
478
479 for
480 (
481 label ptI = feMesh.convexStart();
482 ptI < feMesh.mixedStart();
483// ptI < feMesh.nonFeatureStart();
484 ptI++
485 )
486 {
487 createMasterAndSlavePoints(feMesh, ptI, pts);
488 }
489
490 if (foamyHexMeshControls_.guardFeaturePoints())
491 {
492 for
493 (
494 //label ptI = feMesh.convexStart();
495 label ptI = feMesh.mixedStart();
496 ptI < feMesh.nonFeatureStart();
497 ptI++
498 )
499 {
500 pts.append
501 (
502 Vb
503 (
504 feMesh.points()[ptI],
506 )
507 );
508 }
509 }
510 }
511}
512
513
514// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
515
517(
518 const conformalVoronoiMesh& foamyHexMesh
519)
520:
521 foamyHexMesh_(foamyHexMesh),
522 foamyHexMeshControls_(foamyHexMesh.foamyHexMeshControls()),
523 geometryToConformTo_(foamyHexMesh.geometryToConformTo()),
524 featurePointVertices_(),
525 ftPtPairs_(foamyHexMesh)
526{
527 Info<< nl
528 << "Conforming to feature points" << endl;
529
531 << indent << "Circulating edges is: "
532 << foamyHexMeshControls_.circulateEdges().c_str() << nl
533 << indent << "Guarding feature points is: "
534 << foamyHexMeshControls_.guardFeaturePoints().c_str() << nl
535 << indent << "Snapping to feature points is: "
536 << foamyHexMeshControls_.snapFeaturePoints().c_str() << nl
537 << indent << "Specialising feature points is: "
538 << foamyHexMeshControls_.specialiseFeaturePoints().c_str()
539 << decrIndent
540 << endl;
541
542 DynamicList<Vb> pts;
543
544 createFeaturePoints(pts);
545
546 createMixedFeaturePoints(pts);
547
548 // Points added using the createEdgePointGroup function will be labelled as
549 // internal/external feature edges. Relabel them as feature points,
550 // otherwise they are inserted as both feature points and surface points.
551 forAll(pts, pI)
552 {
553 Vb& pt = pts[pI];
554
555 //if (pt.featureEdgePoint())
556 {
557 if (pt.internalBoundaryPoint())
558 {
560 }
561 else if (pt.externalBoundaryPoint())
562 {
564 }
565 }
566 }
567
568 if (foamyHexMeshControls_.objOutput())
569 {
570 DelaunayMeshTools::writeOBJ("featureVertices.obj", pts);
571 }
572
573 featurePointVertices_.transfer(pts);
574}
575
576
577// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
578
580{}
581
582
583// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
584
586(
587 const backgroundMeshDecomposition& decomposition
588)
589{
590 // Distribute points to their appropriate processor
591 decomposition.distributePoints(featurePointVertices_);
592
593 // Update the processor indices of the points to the new processor number
594 forAll(featurePointVertices_, vI)
595 {
596 featurePointVertices_[vI].procIndex() = Pstream::myProcNo();
597 }
598
599 // Also update the feature point pairs
600}
601
602
604(
605 const Map<label>& oldToNewIndices
606)
607{
608 forAll(featurePointVertices_, vI)
609 {
610 const label currentIndex = featurePointVertices_[vI].index();
611
612 const auto newIndexIter = oldToNewIndices.cfind(currentIndex);
613
614 if (newIndexIter.found())
615 {
616 featurePointVertices_[vI].index() = *newIndexIter;
617 }
618 }
619
620 ftPtPairs_.reIndex(oldToNewIndices);
621}
622
623
624// ************************************************************************* //
label n
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Definition: indexedVertex.H:88
vertexType & type()
bool externalBoundaryPoint() const
bool internalBoundaryPoint() const
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
The Delaunay vertices required to conform to a feature point can be determined upon initialisation be...
void reIndexPointPairs(const Map< label > &oldToNewIndices)
Reindex the feature point pairs using the map.
~featurePointConformer()
Destructor.
void distribute(const backgroundMeshDecomposition &decomposition)
Distribute the feature point vertices according to the.
int myProcNo() const noexcept
Return processor number.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const pointField & points
Namespace for OpenFOAM.
dimensionedScalar sign(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
Definition: List.H:66
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:349
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:342
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:356
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