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