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