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