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

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.17