insertBoundaryConformPointPairs.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 Copyright (C) 2019 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*----------------------------------------------------------------------------*/
28
29#include "CV2D.H"
30
31// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32
33void Foam::CV2D::insertPointPair
34(
36 const point2D& p,
37 const label trii,
38 const label hitSurface
39)
40{
41 if
42 (
43 !meshControls().mirrorPoints()
44 || !insertMirrorPoint(toPoint2D(vit->point()), p)
45 )
46 {
47 pointIndexHit pHit
48 (
49 true,
50 toPoint3D(p),
51 trii
52 );
53
54 vectorField norm(1);
55 qSurf_.geometry()[hitSurface].getNormal
56 (
57 List<pointIndexHit>(1, pHit),
58 norm
59 );
60
61 insertPointPair
62 (
63 meshControls().ppDist(),
64 p,
65 toPoint2D(norm[0])
66 );
67 }
68
70 (
71 CGAL::Filter_iterator
72 <
73 Triangulation::All_vertices_iterator,
74 Triangulation::Infinite_tester
75 >(finite_vertices_end(), vit.predicate(), vit.base())
76 );
77}
78
79
80bool Foam::CV2D::insertPointPairAtIntersection
81(
83 const point2D& defVert,
84 const point2D vertices[],
85 const scalar maxProtSize2
86)
87{
88 bool found = false;
89 point2D interPoint;
90 label interTri = -1;
91 label interHitSurface = -1;
92 scalar interDist2 = 0;
93
94 Face_circulator fcStart = incident_faces(vit);
95 Face_circulator fc = fcStart;
96 label vi = 0;
97
98 do
99 {
100 if (!is_infinite(fc))
101 {
102 pointIndexHit pHit;
103 label hitSurface = -1;
104
105 qSurf_.findSurfaceNearestIntersection
106 (
107 toPoint3D(defVert),
108 toPoint3D(vertices[vi]),
109 pHit,
110 hitSurface
111 );
112
113 if (pHit.hit())
114 {
115 scalar dist2 =
116 magSqr(toPoint2D(pHit.hitPoint()) - vertices[vi]);
117
118 // Check the point is further away than the furthest so far
119 if (dist2 > interDist2)
120 {
121 scalar mps2 = maxProtSize2;
122
123 // If this is a boundary triangle reset the tolerance
124 // to avoid finding a hit point very close to a boundary
125 // vertex
126 if (boundaryTriangle(fc))
127 {
128 mps2 = meshControls().maxNotchLen2();
129 }
130
131 if (dist2 > mps2)
132 {
133 found = true;
134 interPoint = toPoint2D(pHit.hitPoint());
135 interTri = pHit.index();
136 interDist2 = dist2;
137 interHitSurface = hitSurface;
138 }
139 }
140 }
141
142 vi++;
143 }
144 } while (++fc != fcStart);
145
146 if (found)
147 {
148 insertPointPair(vit, interPoint, interTri, interHitSurface);
149 return true;
150 }
151
152 return false;
153}
154
155
156Foam::label Foam::CV2D::insertBoundaryConformPointPairs
157(
158 const fileName& fName
159)
160{
161 label nIntersections = 0;
162
163 for
164 (
165 Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
166 vit != finite_vertices_end();
167 vit++
168 )
169 {
170 // Consider only those points part of point-pairs or near boundary
171 if (!vit->nearOrOnBoundary())
172 {
173 continue;
174 }
175
176 // Counter-clockwise circulator
177 Face_circulator fcStart = incident_faces(vit);
178 Face_circulator fc = fcStart;
179
180 bool infinite = false;
181 bool changed = false;
182
183 do
184 {
185 if (is_infinite(fc))
186 {
187 infinite = true;
188 break;
189 }
190 else if (fc->faceIndex() < Fb::UNCHANGED)
191 {
192 changed = true;
193 break;
194 }
195 } while (++fc != fcStart);
196
197 // If the dual-cell is connected to the infinite point or none of the
198 // faces whose circumcentres it uses have changed ignore
199 if (infinite || !changed) continue;
200
201 fc = fcStart;
202 label nVerts = 0;
203
204 do
205 {
206 vertices[nVerts++] = toPoint2D(circumcenter(fc));
207
208 if (nVerts == maxNvert)
209 {
210 break;
211 }
212 } while (++fc != fcStart);
213
214 // Check if dual-cell has a large number of faces in which case
215 // assumed to be in the far-field and reject
216 if (nVerts == maxNvert) continue;
217
218 // Set n+1 vertex to the first vertex for easy circulating
219 vertices[nVerts] = vertices[0];
220
221 // Convert triangle vertex to OpenFOAM point
222 point2DFromPoint defVert = toPoint2D(vit->point());
223
224 scalar maxProtSize2 = meshControls().maxNotchLen2();
225
226 if (vit->internalOrBoundaryPoint())
227 {
228 // Calculate metrics of the dual-cell
229 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
230
231 // The perimeter of the dual-cell
232 scalar perimeter = 0;
233
234 // Twice the area of the dual-cell
235 scalar areaT2 = 0;
236
237 for (int vi=0; vi<nVerts; vi++)
238 {
239 vector2D edge(vertices[vi+1] - vertices[vi]);
240 perimeter += mag(edge);
241 vector2D otherEdge = defVert - vertices[vi];
242 areaT2 += mag(edge.x()*otherEdge.y() - edge.y()*otherEdge.x());
243 }
244
245 // If the dual-cell is very small reject refinement
246 if (areaT2 < meshControls().minEdgeLen2()) continue;
247
248 // Estimate the cell width
249 scalar cellWidth = areaT2/perimeter;
250
251
252 // Check dimensions of dual-cell
253 /*
254 // Quick rejection of dual-cell refinement based on it's perimeter
255 if (perimeter < 2*meshControls().minCellSize()) continue;
256
257 // Also check the area of the cell and reject refinement
258 // if it is less than that allowed
259 if (areaT2 < meshControls().minCellSize2()) continue;
260
261 // Estimate the cell width and reject refinement if it is less than
262 // that allowed
263 if (cellWidth < 0.5*meshControls().minEdgeLen()) continue;
264 */
265
266 if
267 (
268 perimeter > 2*meshControls().minCellSize()
269 && areaT2 > meshControls().minCellSize2()
270 && cellWidth > 0.5*meshControls().minEdgeLen()
271 )
272 {
273 maxProtSize2 = 0.25*meshControls().maxNotchLen2();
274 }
275 }
276
277 if (insertPointPairAtIntersection(vit, defVert, vertices, maxProtSize2))
278 {
279 nIntersections++;
280 }
281 }
282
283 return nIntersections;
284}
285
286
287void Foam::CV2D::markNearBoundaryPoints()
288{
289 label count = 0;
290 for
291 (
292 Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
293 vit != finite_vertices_end();
294 vit++
295 )
296 {
297 if (vit->internalPoint())
298 {
299 Foam::point vert(toPoint3D(vit->point()));
300
301 pointIndexHit pHit;
302 label hitSurface = -1;
303
304 qSurf_.findSurfaceNearest
305 (
306 vert,
307 4*meshControls().minCellSize2(),
308 pHit,
309 hitSurface
310 );
311
312 if (pHit.hit())
313 {
314 vit->setNearBoundary();
315 ++count;
316 }
317 }
318 }
319
320 Info<< count << " points marked as being near a boundary" << endl;
321}
322
323
324// ************************************************************************* //
bool found
const cv2DControls & meshControls() const
Definition: CV2DI.H:119
Foam::point toPoint3D(const point2D &) const
Definition: CV2DI.H:142
const point2D & toPoint2D(const Foam::point &) const
Definition: CV2DI.H:125
Triangulation::Finite_vertices_iterator Finite_vertices_iterator
Definition: DelaunayMesh.H:73
const searchableSurfaces & geometry() const
Return reference to the searchableSurfaces object containing all.
volScalarField & p
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
Definition: BitOps.H:78
label otherEdge(const primitiveMesh &mesh, const labelList &edgeLabels, const label thisEdgeI, const label thisVertI)
Return label of other edge (out of candidates edgeLabels)
Definition: meshTools.C:521
Vector2D< scalar > vector2D
A 2D vector of scalars obtained from the generic Vector2D.
Definition: vector2D.H:51
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
bool boundaryTriangle(const CV2D::Face_handle fc)
Definition: CV2DI.H:206
messageStream Info
Information stream (stdout output on master, null elsewhere)
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)
pointField vertices(const blockVertexList &bvl)
Field< vector > vectorField
Specialisation of Field<T> for vector.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)