insertFeaturePoints.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-2015 OpenFOAM Foundation
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*----------------------------------------------------------------------------*/
27 
28 #include "CV2D.H"
29 #include "plane.H"
30 #include "unitConversion.H"
31 
32 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
33 
34 bool Foam::CV2D::on2DLine(const point2D& p, const linePointRef& line)
35 {
36  const point2D& a = toPoint2D(line.start());
37  const point2D& b = toPoint2D(line.end());
38 
39  if
40  (
41  p.x() < min(a.x(), b.x())
42  || p.x() > max(a.x(), b.x())
43  || p.y() < min(a.y(), b.y())
44  || p.y() > max(a.y(), b.y())
45  )
46  {
47  return false;
48  }
49 
50  return true;
51 }
52 
53 
54 void Foam::CV2D::insertFeaturePoints()
55 {
56  featurePoints_.clear();
57  label nVert = number_of_vertices();
58 
59  const PtrList<extendedFeatureEdgeMesh>& feMeshes
60  (
61  qSurf_.features()
62  );
63 
64  if (feMeshes.empty())
65  {
66  WarningInFunction
67  << "Extended Feature Edge Mesh is empty so no feature points will "
68  << "be found." << nl
69  << " Use: featureMethod extendedFeatureEdgeMesh;" << nl
70  << endl;
71  }
72 
73  forAll(feMeshes, i)
74  {
75  const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
76  const edgeList& edges = feMesh.edges();
77  const pointField& points = feMesh.points();
78 
79  if (debug)
80  {
81  label nConvex = feMesh.concaveStart() - feMesh.convexStart();
82  label nConcave = feMesh.mixedStart() - feMesh.concaveStart();
83  label nMixed = feMesh.nonFeatureStart() - feMesh.mixedStart();
84  label nExternal = feMesh.internalStart() - feMesh.externalStart();
85  label nInternal = feMesh.flatStart() - feMesh.internalStart();
86  label nFlat = feMesh.openStart() - feMesh.flatStart();
87  label nOpen = feMesh.multipleStart() - feMesh.openStart();
88  label nMultiple = edges.size() - feMesh.multipleStart();
89 
90  Info<< "Inserting Feature Points:" << nl
91  << " Convex points: " << nConvex << nl
92  << " Concave points: " << nConcave << nl
93  << " Mixed points: " << nMixed << nl
94  << " External edges: " << nExternal << nl
95  << " Internal edges: " << nInternal << nl
96  << " Flat edges: " << nFlat << nl
97  << " Open edges: " << nOpen << nl
98  << " Multiple edges: " << nMultiple << endl;
99  }
100 
101  // Args: (base point, normal)
102  // TODO: allow user to input this
103  plane zPlane(vector(0, 0, z_), vector(0, 0, 1));
104 
105  if (debug)
106  {
107  Info<< " plane: " << zPlane << " " << z_ << endl;
108  }
109 
110  forAll(edges, edgeI)
111  {
112  const edge& e = feMesh.edges()[edgeI];
113 
114  const Foam::point& ep0 = points[e.start()];
115  const Foam::point& ep1 = points[e.end()];
116 
117  const linePointRef line(ep0, ep1);
118 
119  scalar intersect = zPlane.lineIntersect(line);
120 
121  point2D featPoint = toPoint2D(intersect * (ep1 - ep0) + ep0);
122 
123  if (on2DLine(featPoint, line))
124  {
125  vector2DField fpn = toPoint2D(feMesh.edgeNormals(edgeI));
126 
127  vector2D cornerNormal = sum(fpn);
128  cornerNormal.normalise();
129 
130  if (debug)
131  {
132  Info<< nl << " line: " << line << nl
133  << " vec: " << line.vec() << nl
134  << " featurePoint: " << featPoint << nl
135  << " line length: " << line.mag() << nl
136  << " intersect: " << intersect << endl;
137  }
138 
139  if
140  (
141  feMesh.getEdgeStatus(edgeI)
142  == extendedFeatureEdgeMesh::EXTERNAL
143  )
144  {
145  // Convex Point
146  Foam::point2D internalPt =
147  featPoint - meshControls().ppDist()*cornerNormal;
148 
149  if (debug)
150  {
151  Info<< "PREC: " << internalPt << nl
152  << " : " << featPoint << nl
153  << " : " << meshControls().ppDist() << nl
154  << " : " << cornerNormal << endl;
155  }
156 
157  featurePoints_.push_back
158  (
159  Vb
160  (
161  toPoint(internalPt),
162  nVert,
163  nVert + 1
164  )
165  );
166  label masterPtIndex = nVert++;
167 
168  forAll(fpn, nI)
169  {
170  const vector n3D(fpn[nI][0], fpn[nI][1], 0.0);
171 
172  plane planeN = plane(toPoint3D(featPoint), n3D);
173 
174  Foam::point2D externalPt =
175  internalPt
176  + (
177  2.0
178  * planeN.distance(toPoint3D(internalPt))
179  * fpn[nI]
180  );
181 
182  featurePoints_.push_back
183  (
184  Vb
185  (
186  toPoint(externalPt),
187  nVert++,
188  masterPtIndex
189  )
190  );
191 
192  if (debug)
193  {
194  Info<< " side point: " << externalPt << endl;
195  }
196  }
197 
198  if (debug)
199  {
200  Info<< "Convex Point: " << featPoint << nl
201  << " corner norm: " << cornerNormal << nl
202  << " reference: " << internalPt << endl;
203  }
204  }
205  else if
206  (
207  feMesh.getEdgeStatus(edgeI)
208  == extendedFeatureEdgeMesh::INTERNAL
209  )
210  {
211  // Concave Point
212  Foam::point2D externalPt =
213  featPoint + meshControls().ppDist()*cornerNormal;
214 
215  Foam::point2D refPt =
216  featPoint - meshControls().ppDist()*cornerNormal;
217 
218  label slavePointIndex = 0;
219 
220  scalar totalAngle =
221  radToDeg
222  (
223  constant::mathematical::pi
224  + acos(mag(fpn[0] & fpn[1]))
225  );
226 
227  // Number of quadrants the angle should be split into
228  int nQuads =
229  int(totalAngle/meshControls().maxQuadAngle()) + 1;
230 
231  // The number of additional master points needed to
232  //obtain the required number of quadrants.
233  int nAddPoints = min(max(nQuads - 2, 0), 2);
234 
235  // index of reflMaster
236  label reflectedMaster = nVert + 2 + nAddPoints;
237 
238  if (debug)
239  {
240  Info<< "Concave Point: " << featPoint << nl
241  << " corner norm: " << cornerNormal << nl
242  << " external: " << externalPt << nl
243  << " reference: " << refPt << nl
244  << " angle: " << totalAngle << nl
245  << " nQuads: " << nQuads << nl
246  << " nAddPoints: " << nAddPoints << endl;
247  }
248 
249  forAll(fpn, nI)
250  {
251  const vector n3D(fpn[nI][0], fpn[nI][1], 0.0);
252 
253  plane planeN = plane(toPoint3D(featPoint), n3D);
254 
255  Foam::point2D internalPt =
256  externalPt
257  - (
258  2.0
259  * planeN.distance(toPoint3D(externalPt))
260  * fpn[nI]
261  );
262 
263  featurePoints_.push_back
264  (
265  Vb
266  (
267  toPoint(internalPt),
268  nVert,
269  reflectedMaster
270  )
271  );
272  slavePointIndex = nVert++;
273 
274  if (debug)
275  {
276  Info<< "Internal Point: " << internalPt << endl;
277  }
278  }
279 
280  if (nAddPoints == 1)
281  {
282  // One additional point is the reflection of the slave
283  // point, i.e., the original reference point
284  featurePoints_.push_back
285  (
286  Vb
287  (
288  toPoint(refPt),
289  nVert++,
290  reflectedMaster
291  )
292  );
293 
294  if (debug)
295  {
296  Info<< "ref Point: " << refPt << endl;
297  }
298  }
299  else if (nAddPoints == 2)
300  {
301  point2D reflectedAa =
302  refPt - ((featPoint - externalPt) & fpn[1])*fpn[1];
303 
304  featurePoints_.push_back
305  (
306  Vb
307  (
308  toPoint(reflectedAa),
309  nVert++,
310  reflectedMaster
311  )
312  );
313 
314  point2D reflectedBb =
315  refPt - ((featPoint - externalPt) & fpn[0])*fpn[0];
316 
317  featurePoints_.push_back
318  (
319  Vb
320  (
321  toPoint(reflectedBb),
322  nVert++,
323  reflectedMaster
324  )
325  );
326 
327  if (debug)
328  {
329  Info<< "refA Point: " << reflectedAa << nl
330  << "refb Point: " << reflectedBb << endl;
331  }
332  }
333 
334  featurePoints_.push_back
335  (
336  Vb
337  (
338  toPoint(externalPt),
339  nVert++,
340  slavePointIndex
341  )
342  );
343  }
344  else
345  {
346  WarningInFunction
347  << "Feature Edge " << edges[edgeI] << nl
348  << " points(" << points[edges[edgeI].start()]
349  << ", " << points[edges[edgeI].end()] << ")" << nl
350  << " is not labelled as either concave or convex, it"
351  << " is labelled as (#2 = flat): "
352  << feMesh.getEdgeStatus(edgeI) << endl;
353  }
354  }
355  else
356  {
357  WarningInFunction
358  << "Point " << featPoint << " is not on the line "
359  << line << endl;
360  }
361  }
362  }
363 
364  // Insert the feature points.
365  reinsertFeaturePoints();
366 
367  if (meshControls().objOutput())
368  {
369  writePoints("feat_allPoints.obj", false);
370  writeFaces("feat_allFaces.obj", false);
371  writeFaces("feat_faces.obj", true);
372  writeTriangles("feat_triangles.obj", true);
373  }
374 }
375 
376 
377 void Foam::CV2D::reinsertFeaturePoints()
378 {
379  for
380  (
381  std::list<Vb>::iterator vit=featurePoints_.begin();
382  vit != featurePoints_.end();
383  ++vit
384  )
385  {
386  insertPoint
387  (
388  toPoint2D(vit->point()),
389  vit->index(),
390  vit->type()
391  );
392  }
393 }
394 
395 
396 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::edgeList
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
CGAL::indexedVertex
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Definition: indexedVertex.H:54
unitConversion.H
Unit conversion functions.
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::Vector2D< scalar >
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::CV2D::toPoint2D
const point2D & toPoint2D(const Foam::point &) const
Definition: CV2DI.H:125
Foam::extendedEdgeMesh::INTERNAL
"Concave" edge
Definition: extendedEdgeMesh.H:107
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
plane.H
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::radToDeg
constexpr scalar radToDeg(const scalar rad) noexcept
Conversion from radians to degrees.
Definition: unitConversion.H:54
Foam::vector2DField
Field< vector2D > vector2DField
Forward declarations of the specialisation of Field<T> for vector2D.
Definition: vector2DFieldFwd.H:49
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::toPoint
PointFrompoint toPoint(const Foam::point &p)
Definition: pointConversion.H:82
Foam::point2D
vector2D point2D
Point2D is a vector.
Definition: point2D.H:43
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::linePointRef
line< point, const point & > linePointRef
A line using referred points.
Definition: linePointRef.H:47
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::vector2D
Vector2D< scalar > vector2D
A 2D vector of scalars obtained from the generic Vector2D.
Definition: vector2D.H:51
Foam::Vector
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
Definition: Vector.H:62
Foam::acos
dimensionedScalar acos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:268
Foam::extendedEdgeMesh::EXTERNAL
"Convex" edge
Definition: extendedEdgeMesh.H:106
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
CV2D.H