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-------------------------------------------------------------------------------
10License
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
34bool 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
54void 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 {
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)
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)
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 =
222 (
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 {
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 {
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
377void 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// ************************************************************************* //
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Definition: indexedVertex.H:88
const point2D & toPoint2D(const Foam::point &) const
Definition: CV2DI.H:125
iterator end() noexcept
Return an iterator to end traversing the UList.
Definition: UListI.H:350
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
volScalarField & p
const pointField & points
#define WarningInFunction
Report a warning using Foam::Warning.
constexpr scalar pi(M_PI)
Vector2D< scalar > vector2D
A 2D vector of scalars obtained from the generic Vector2D.
Definition: vector2D.H:51
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
PointFrompoint toPoint(const Foam::point &p)
line< point, const point & > linePointRef
A line using referred points.
Definition: linePointRef.H:47
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Field< vector2D > vector2DField
Forward declarations of the specialisation of Field<T> for vector2D.
vector2D point2D
Point2D is a vector.
Definition: point2D.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
constexpr scalar radToDeg() noexcept
Multiplication factor for radians to degrees conversion.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
Vector< scalar > vector
Definition: vector.H:61
dimensionedScalar acos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
volScalarField & b
Definition: createFields.H:27
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Unit conversion functions.