insertSurfaceNearestPointPairs.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
30// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31
32bool Foam::CV2D::dualCellSurfaceIntersection
33(
35) const
36{
37 Triangulation::Edge_circulator ecStart = incident_edges(vit);
38 Triangulation::Edge_circulator ec = ecStart;
39
40 do
41 {
42 if (!is_infinite(ec))
43 {
44 Foam::point e0 = toPoint3D(circumcenter(ec->first));
45
46 // If edge end is outside bounding box then edge cuts boundary
47 if (!qSurf_.globalBounds().contains(e0))
48 {
49 return true;
50 }
51
52 Foam::point e1 =
53 toPoint3D(circumcenter(ec->first->neighbor(ec->second)));
54
55 // If other edge end is outside bounding box then edge cuts boundary
56 if (!qSurf_.globalBounds().contains(e1))
57 {
58 return true;
59 }
60
61 if (magSqr(e1 - e0) > meshControls().minEdgeLen2())
62 {
63 if (qSurf_.findSurfaceAnyIntersection(e0, e1))
64 {
65 return true;
66 }
67 }
68 }
69
70 } while (++ec != ecStart);
71
72 return false;
73}
74
75
76void Foam::CV2D::insertPointPairs
77(
78 const DynamicList<point2D>& nearSurfacePoints,
79 const DynamicList<point2D>& surfacePoints,
80 const DynamicList<label>& surfaceTris,
81 const DynamicList<label>& surfaceHits,
82 const fileName fName
83)
84{
85 if (meshControls().mirrorPoints())
86 {
87 forAll(surfacePoints, ppi)
88 {
89 insertMirrorPoint
90 (
91 nearSurfacePoints[ppi],
92 surfacePoints[ppi]
93 );
94 }
95 }
96 else
97 {
98 forAll(surfacePoints, ppi)
99 {
100 pointIndexHit pHit
101 (
102 true,
103 toPoint3D(surfacePoints[ppi]),
104 surfaceTris[ppi]
105 );
106
107 vectorField norm(1);
108 qSurf_.geometry()[surfaceHits[ppi]].getNormal
109 (
110 List<pointIndexHit>(1, pHit),
111 norm
112 );
113
114 insertPointPair
115 (
116 meshControls().ppDist(),
117 surfacePoints[ppi],
118 toPoint2D(norm[0])
119 );
120 }
121 }
122
123 Info<< surfacePoints.size() << " point-pairs inserted" << endl;
124
125 if (meshControls().objOutput())
126 {
127 OFstream str(fName);
128 label vertI = 0;
129
130 forAll(surfacePoints, ppi)
131 {
132 meshTools::writeOBJ(str, toPoint3D(surfacePoints[ppi]));
133 vertI++;
134 }
135
136 Info<< "insertPointPairs: Written " << surfacePoints.size()
137 << " inserted point-pair locations to file "
138 << str.name() << endl;
139 }
140}
141
142
143void Foam::CV2D::insertSurfaceNearestPointPairs()
144{
145 Info<< "insertSurfaceNearestPointPairs: ";
146
147 label nSurfacePointsEst =
148 min
149 (
150 label(number_of_vertices()),
151 label(10*sqrt(scalar(number_of_vertices())))
152 );
153
154 DynamicList<point2D> nearSurfacePoints(nSurfacePointsEst);
155 DynamicList<point2D> surfacePoints(nSurfacePointsEst);
156 DynamicList<label> surfaceTris(nSurfacePointsEst);
157 DynamicList<label> surfaceHits(nSurfacePointsEst);
158
159 // Local references to surface mesh addressing
160// const pointField& localPoints = qSurf_.localPoints();
161// const labelListList& edgeFaces = qSurf_.edgeFaces();
162// const vectorField& faceNormals = qSurf_.faceNormals();
163// const labelListList& faceEdges = qSurf_.faceEdges();
164
165 for
166 (
167 Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
168 vit != finite_vertices_end();
169 vit++
170 )
171 {
172 if (vit->internalPoint())
173 {
174 point2DFromPoint vert(toPoint2D(vit->point()));
175
176 pointIndexHit pHit;
177 label hitSurface = -1;
178
179 qSurf_.findSurfaceNearest
180 (
181 toPoint3D(vert),
182 4*meshControls().minCellSize2(),
183 pHit,
184 hitSurface
185 );
186
187 if (pHit.hit())
188 {
189 vit->setNearBoundary();
190
191 // Reference to the nearest triangle
192// const labelledTri& f = qSurf_[hitSurface];
193
194// // Find where point is on triangle.
195// // Note tolerance needed is relative one
196// // (used in comparing normalized [0..1] triangle coordinates).
197// label nearType, nearLabel;
198// triPointRef
199// (
200// localPoints[f[0]],
201// localPoints[f[1]],
202// localPoints[f[2]]
203// ).classify(pHit.hitPoint(), nearType, nearLabel);
204
205// // If point is on a edge check if it is an internal feature
206
207// bool internalFeatureEdge = false;
208
209// if (nearType == triPointRef::EDGE)
210// {
211// label edgeI = faceEdges[hitSurface][nearLabel];
212// const labelList& eFaces = edgeFaces[edgeI];
213
214// if
215// (
216// eFaces.size() == 2
217// && (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
218// < -0.2
219// )
220// {
221// internalFeatureEdge = true;
222// }
223// }
224
225 if (dualCellSurfaceIntersection(vit)) //&& !internalFeatureEdge)
226 {
227 nearSurfacePoints.append(vert);
228 surfacePoints.append(toPoint2D(pHit.hitPoint()));
229 surfaceTris.append(pHit.index());
230 surfaceHits.append(hitSurface);
231 }
232 }
233 }
234 }
235
236 insertPointPairs
237 (
238 nearSurfacePoints,
239 surfacePoints,
240 surfaceTris,
241 surfaceHits,
242 "surfaceNearestIntersections.obj"
243 );
244}
245
246
247// ************************************************************************* //
const cv2DControls & meshControls() const
Definition: CV2DI.H:119
Foam::point toPoint3D(const point2D &) const
Definition: CV2DI.H:142
Triangulation::Finite_vertices_iterator Finite_vertices_iterator
Definition: DelaunayMesh.H:73
bool findSurfaceAnyIntersection(const point &start, const point &end) const
const treeBoundBox & globalBounds() const
Return the global bounds.
bool contains(const vector &dir, const point &) const
Contains point (inside or on edge) and moving in direction.
Definition: treeBoundBox.C:320
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensionedScalar sqrt(const dimensionedScalar &ds)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Field< vector > vectorField
Specialisation of Field<T> for vector.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333