DelaunayMesh.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) 2012-2015 OpenFOAM Foundation
9 Copyright (C) 2021 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 "DelaunayMesh.H"
30#include "polyMesh.H"
31#include "labelPairHashes.H"
32#include "PrintTable.H"
33#include "pointIOField.H"
34#include "scalarIOField.H"
35#include "labelIOField.H"
36#include "pointConversion.H"
37#include <algorithm>
38#include <random>
39
40// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
41
42template<class Triangulation>
44:
45 Triangulation(),
46 vertexCount_(0),
47 cellCount_(0),
48 runTime_(runTime)
49{}
50
51
52template<class Triangulation>
54(
55 const Time& runTime,
56 const word& meshName
57)
58:
59 Triangulation(),
60 vertexCount_(0),
61 cellCount_(0),
62 runTime_(runTime)
63{
64 Info<< "Reading " << meshName << " from " << runTime.timeName() << endl;
65
66 pointIOField pts
67 (
68 IOobject
69 (
70 "points",
71 runTime.timeName(),
72 meshName/polyMesh::meshSubDir,
73 runTime,
74 IOobject::READ_IF_PRESENT,
75 IOobject::NO_WRITE
76 )
77 );
78
79 if (pts.typeHeaderOk<pointIOField>(true))
80 {
81 labelIOField types
82 (
83 IOobject
84 (
85 "types",
86 runTime.timeName(),
87 meshName,
88 runTime,
89 IOobject::MUST_READ,
90 IOobject::NO_WRITE
91 )
92 );
93
94// Do not read in indices
95// labelIOField indices
96// (
97// IOobject
98// (
99// "indices",
100// runTime.timeName(),
101// meshName,
102// runTime,
103// IOobject::MUST_READ,
104// IOobject::NO_WRITE
105// )
106// );
107
108 labelIOField processorIndices
109 (
110 IOobject
111 (
112 "processorIndices",
113 runTime.timeName(),
114 meshName,
115 runTime,
116 IOobject::MUST_READ,
117 IOobject::NO_WRITE
118 )
119 );
120
121 List<Vb> pointsToInsert(pts.size());
122
123 forAll(pointsToInsert, pI)
124 {
125 pointsToInsert[pI] =
126 Vb
127 (
128 toPoint(pts[pI]),
129 pI,
130 static_cast<indexedVertexEnum::vertexType>(types[pI]),
131 processorIndices[pI]
132 );
133 }
134
135 rangeInsertWithInfo
136 (
137 pointsToInsert.begin(),
138 pointsToInsert.end(),
139 false,
140 false
141 );
142
143 vertexCount_ = Triangulation::number_of_vertices();
144 }
145}
146
147
148// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
149
150template<class Triangulation>
152{}
153
154
155// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
156
157
158// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
159
160
161// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
162
163template<class Triangulation>
165{
166 Info<< "Clearing triangulation" << endl;
167
168 DynamicList<Vb> vertices;
169
170 for
171 (
172 Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
173 vit != Triangulation::finite_vertices_end();
174 ++vit
175 )
176 {
177 if (vit->fixed())
178 {
180 (
181 Vb
182 (
183 vit->point(),
184 vit->index(),
185 vit->type(),
186 vit->procIndex()
187 )
188 );
189
190 vertices.last().fixed() = vit->fixed();
191 }
192 }
193
194 this->clear();
195
196 resetVertexCount();
197 resetCellCount();
198
199 insertPoints(vertices, false);
200
201 Info<< "Inserted " << vertexCount() << " fixed points" << endl;
202}
203
204
205template<class Triangulation>
207(
208 const List<Vb>& vertices,
209 const bool reIndex
210)
211{
212 return rangeInsertWithInfo
213 (
214 vertices.begin(),
215 vertices.end(),
216 false,
217 reIndex
218 );
219}
220
221
222template<class Triangulation>
224operator()
225(
226 const Point_3& p,
227 const Point_3& q
228) const
229{
230 return typename Gt::Less_x_3()(*(p.first), *(q.first));
231}
232
233template<class Triangulation>
235operator()
236(
237 const Point_3& p,
238 const Point_3& q
239) const
240{
241 return typename Gt::Less_y_3()(*(p.first), *(q.first));
242}
243
244template<class Triangulation>
246operator()
247(
248 const Point_3& p,
249 const Point_3& q
250) const
251{
252 return typename Gt::Less_z_3()(*(p.first), *(q.first));
253}
254
255template<class Triangulation>
258const
259{
260 return Less_x_3();
261}
262
263template<class Triangulation>
266const
267{
268 return Less_y_3();
269}
270
271template<class Triangulation>
274const
275{
276 return Less_z_3();
277}
278
279
280template<class Triangulation>
281template<class PointIterator>
283(
284 PointIterator begin,
285 PointIterator end,
286 bool printErrors,
287 bool reIndex
288)
289{
290 typedef DynamicList
291 <
292 std::pair
293 <
294 const typename Triangulation::Point*,
295 label
296 >
297 > vectorPairPointIndex;
298
299 vectorPairPointIndex points;
300
301 label count = 0;
302 for (PointIterator it = begin; it != end; ++it)
303 {
304 points.append
305 (
306 std::make_pair(&(it->point()), count++)
307 );
308 }
309
310 std::shuffle(points.begin(), points.end(), std::default_random_engine());
311
312 spatial_sort
313 (
314 points.begin(),
315 points.end(),
316 Traits_for_spatial_sort()
317 );
318
319 Vertex_handle hint;
320
321 Map<label> oldToNewIndex(points.size());
322
323 for
324 (
325 typename vectorPairPointIndex::const_iterator p = points.begin();
326 p != points.end();
327 ++p
328 )
329 {
330 const size_t checkInsertion = Triangulation::number_of_vertices();
331
332 hint = this->insert(*(p->first), hint);
333
334 const Vb& vert = *(begin + p->second);
335
336 if (checkInsertion != Triangulation::number_of_vertices() - 1)
337 {
338 if (printErrors)
339 {
340 Vertex_handle nearV =
341 Triangulation::nearest_vertex(*(p->first));
342
343 Pout<< "Failed insertion : " << vert.info()
344 << " nearest : " << nearV->info();
345 }
346 }
347 else
348 {
349 const label oldIndex = vert.index();
350 hint->index() = getNewVertexIndex();
351
352 if (reIndex)
353 {
354 oldToNewIndex.insert(oldIndex, hint->index());
355 }
356
357 hint->type() = vert.type();
358 hint->procIndex() = vert.procIndex();
359 hint->targetCellSize() = vert.targetCellSize();
360 hint->alignment() = vert.alignment();
361 }
362 }
363
364 return oldToNewIndex;
365}
366
367
368// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
369
370#include "DelaunayMeshIO.C"
371
372// ************************************************************************* //
CGAL::indexedVertex< K > Vb
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Definition: indexedVertex.H:88
Foam::scalar & targetCellSize()
vertexType & type()
Foam::label & index()
Foam::InfoProxy< indexedVertex< Gt, Vb > > info() const
Info proxy, to print information to a stream.
Foam::tensor & alignment()
The vertex and cell classes must have an index defined.
Definition: DelaunayMesh.H:63
Map< label > insertPoints(const List< Vb > &vertices, const bool reIndex)
Insert the list of vertices (calls rangeInsertWithInfo)
Map< label > rangeInsertWithInfo(PointIterator begin, PointIterator end, bool printErrors=false, bool reIndex=true)
Function inserting points into a triangulation and setting the.
~DelaunayMesh()
Destructor.
void reset()
Clear the entire triangulation.
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
iterator begin() noexcept
Return an iterator to begin traversing the UList.
Definition: UListI.H:329
iterator end() noexcept
Return an iterator to end traversing the UList.
Definition: UListI.H:350
T & last()
Return the last element of the list.
Definition: UListI.H:216
CellSizeDelaunay::Point Point
friend Ostream & operator(Ostream &, const faMatrix< Type > &)
volScalarField & p
patchWriters clear()
engineTime & runTime
const pointField & points
A HashTable to objects of type <T> with a labelPair key. The hashing is based on labelPair (FixedList...
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
Definition: BitOps.H:78
PointFrompoint toPoint(const Foam::point &p)
messageStream Info
Information stream (stdout output on master, null elsewhere)
IOField< label > labelIOField
labelField with IO.
Definition: labelIOField.H:44
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
vectorIOField pointIOField
pointIOField is a vectorIOField.
Definition: pointIOField.H:44
pointField vertices(const blockVertexList &bvl)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:158
constexpr auto begin(C &c) -> decltype(c.begin())
Return iterator to the beginning of the container c.
Definition: stdFoam.H:134
srcOptions insert("case", fileName(rootDirSource/caseDirSource))
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333