insertBoundaryConformPointPairs.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 void Foam::CV2D::insertPointPair
31 (
32  Triangulation::Finite_vertices_iterator& vit,
33  const point2D& p,
34  const label trii,
35  const label hitSurface
36 )
37 {
38  if
39  (
40  !meshControls().mirrorPoints()
41  || !insertMirrorPoint(toPoint2D(vit->point()), p)
42  )
43  {
44  pointIndexHit pHit
45  (
46  true,
47  toPoint3D(p),
48  trii
49  );
50 
51  vectorField norm(1);
52  qSurf_.geometry()[hitSurface].getNormal
53  (
54  List<pointIndexHit>(1, pHit),
55  norm
56  );
57 
58  insertPointPair
59  (
60  meshControls().ppDist(),
61  p,
62  toPoint2D(norm[0])
63  );
64  }
65 
66  vit = Triangulation::Finite_vertices_iterator
67  (
68  CGAL::Filter_iterator
69  <
70  Triangulation::All_vertices_iterator,
71  Triangulation::Infinite_tester
72  >(finite_vertices_end(), vit.predicate(), vit.base())
73  );
74 }
75 
76 
77 bool Foam::CV2D::insertPointPairAtIntersection
78 (
79  Triangulation::Finite_vertices_iterator& vit,
80  const point2D& defVert,
81  const point2D vertices[],
82  const scalar maxProtSize2
83 )
84 {
85  bool found = false;
86  point2D interPoint;
87  label interTri = -1;
88  label interHitSurface = -1;
89  scalar interDist2 = 0;
90 
91  Face_circulator fcStart = incident_faces(vit);
92  Face_circulator fc = fcStart;
93  label vi = 0;
94 
95  do
96  {
97  if (!is_infinite(fc))
98  {
99  pointIndexHit pHit;
100  label hitSurface = -1;
101 
103  (
104  toPoint3D(defVert),
105  toPoint3D(vertices[vi]),
106  pHit,
107  hitSurface
108  );
109 
110  if (pHit.hit())
111  {
112  scalar dist2 =
113  magSqr(toPoint2D(pHit.hitPoint()) - vertices[vi]);
114 
115  // Check the point is further away than the furthest so far
116  if (dist2 > interDist2)
117  {
118  scalar mps2 = maxProtSize2;
119 
120  // If this is a boundary triangle reset the tolerance
121  // to avoid finding a hit point very close to a boundary
122  // vertex
123  if (boundaryTriangle(fc))
124  {
125  mps2 = meshControls().maxNotchLen2();
126  }
127 
128  if (dist2 > mps2)
129  {
130  found = true;
131  interPoint = toPoint2D(pHit.hitPoint());
132  interTri = pHit.index();
133  interDist2 = dist2;
134  interHitSurface = hitSurface;
135  }
136  }
137  }
138 
139  vi++;
140  }
141  } while (++fc != fcStart);
142 
143  if (found)
144  {
145  insertPointPair(vit, interPoint, interTri, interHitSurface);
146  return true;
147  }
148  else
149  {
150  return false;
151  }
152 }
153 
154 
155 Foam::label Foam::CV2D::insertBoundaryConformPointPairs
156 (
157  const fileName& fName
158 )
159 {
160  label nIntersections = 0;
161 
162  for
163  (
164  Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
165  vit != finite_vertices_end();
166  vit++
167  )
168  {
169  // Consider only those points part of point-pairs or near boundary
170  if (!vit->nearOrOnBoundary())
171  {
172  continue;
173  }
174 
175  // Counter-clockwise circulator
176  Face_circulator fcStart = incident_faces(vit);
177  Face_circulator fc = fcStart;
178 
179  bool infinite = false;
180  bool changed = false;
181 
182  do
183  {
184  if (is_infinite(fc))
185  {
186  infinite = true;
187  break;
188  }
189  else if (fc->faceIndex() < Fb::UNCHANGED)
190  {
191  changed = true;
192  break;
193  }
194  } while (++fc != fcStart);
195 
196  // If the dual-cell is connected to the infinite point or none of the
197  // faces whose circumcentres it uses have changed ignore
198  if (infinite || !changed) continue;
199 
200  fc = fcStart;
201  label nVerts = 0;
202 
203  do
204  {
205  vertices[nVerts++] = toPoint2D(circumcenter(fc));
206 
207  if (nVerts == maxNvert)
208  {
209  break;
210  }
211  } while (++fc != fcStart);
212 
213  // Check if dual-cell has a large number of faces in which case
214  // assumed to be in the far-field and reject
215  if (nVerts == maxNvert) continue;
216 
217  // Set n+1 vertex to the first vertex for easy circulating
218  vertices[nVerts] = vertices[0];
219 
220  // Convert triangle vertex to OpenFOAM point
221  point2DFromPoint defVert = toPoint2D(vit->point());
222 
223  scalar maxProtSize2 = meshControls().maxNotchLen2();
224 
225  if (vit->internalOrBoundaryPoint())
226  {
227  // Calculate metrics of the dual-cell
228  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
229 
230  // The perimeter of the dual-cell
231  scalar perimeter = 0;
232 
233  // Twice the area of the dual-cell
234  scalar areaT2 = 0;
235 
236  for (int vi=0; vi<nVerts; vi++)
237  {
238  vector2D edge(vertices[vi+1] - vertices[vi]);
239  perimeter += mag(edge);
240  vector2D otherEdge = defVert - vertices[vi];
241  areaT2 += mag(edge.x()*otherEdge.y() - edge.y()*otherEdge.x());
242  }
243 
244  // If the dual-cell is very small reject refinement
245  if (areaT2 < meshControls().minEdgeLen2()) continue;
246 
247  // Estimate the cell width
248  scalar cellWidth = areaT2/perimeter;
249 
250 
251  // Check dimensions of dual-cell
252  /*
253  // Quick rejection of dual-cell refinement based on it's perimeter
254  if (perimeter < 2*meshControls().minCellSize()) continue;
255 
256  // Also check the area of the cell and reject refinement
257  // if it is less than that allowed
258  if (areaT2 < meshControls().minCellSize2()) continue;
259 
260  // Estimate the cell width and reject refinement if it is less than
261  // that allowed
262  if (cellWidth < 0.5*meshControls().minEdgeLen()) continue;
263  */
264 
265  if
266  (
267  perimeter > 2*meshControls().minCellSize()
268  && areaT2 > meshControls().minCellSize2()
269  && cellWidth > 0.5*meshControls().minEdgeLen()
270  )
271  {
272  maxProtSize2 = 0.25*meshControls().maxNotchLen2();
273  }
274  }
275 
276  if (insertPointPairAtIntersection(vit, defVert, vertices, maxProtSize2))
277  {
278  nIntersections++;
279  }
280  }
281 
282  return nIntersections;
283 }
284 
285 
286 void Foam::CV2D::markNearBoundaryPoints()
287 {
288  label count = 0;
289  for
290  (
291  Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
292  vit != finite_vertices_end();
293  vit++
294  )
295  {
296  if (vit->internalPoint())
297  {
298  point vert(toPoint3D(vit->point()));
299 
300  pointIndexHit pHit;
301  label hitSurface = -1;
302 
303  qSurf_.findSurfaceNearest
304  (
305  vert,
306  4*meshControls().minCellSize2(),
307  pHit,
308  hitSurface
309  );
310 
311  if (pHit.hit())
312  {
313  vit->setNearBoundary();
314  ++count;
315  }
316  }
317  }
318 
319  Info<< count << " points marked as being near a boundary" << endl;
320 }
321 
322 
323 // ************************************************************************* //
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
vector2D point2D
Point2D is a vector.
Definition: point2D.H:41
Vector2D< scalar > vector2D
vector2D obtained from generic Vector2D
Definition: vector2D.H:49
scalar maxNotchLen2() const
Return the maxNotchLen squared.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
point toPoint3D(const point2D &) const
Definition: CV2DI.H:141
const cv2DControls & meshControls() const
Definition: CV2DI.H:118
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
const point2D & toPoint2D(const point &) const
Definition: CV2DI.H:124
void findSurfaceNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &surfHit, label &hitSurface) const
Find the nearest point to the sample and return it to the.
const point2D & point2DFromPoint
Definition: CV2D.H:354
dimensioned< scalar > magSqr(const dimensioned< Type > &)
scalar minCellSize2() const
Return the square of the minimum cell size.
Definition: cv2DControlsI.H:32
scalar minEdgeLen2() const
Return the minEdgeLen squared.
scalar minCellSize() const
Return the minimum cell size.
Definition: cv2DControlsI.H:26
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
scalar minEdgeLen() const
Return the minEdgeLen.
void findSurfaceNearestIntersection(const point &start, const point &end, pointIndexHit &surfHit, label &hitSurface) const
Finding the nearestIntersection of the surface to start.
Field< vector > vectorField
Specialisation of Field<T> for vector.
bool boundaryTriangle(const CV2D::Face_handle fc)
Definition: CV2DI.H:205