smoothAlignmentSolver.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-2016 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
29
30// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31
32template<class Triangulation, class Type>
33Foam::tmp<Foam::Field<Type>> Foam::smoothAlignmentSolver::filterFarPoints
34(
35 const Triangulation& mesh,
36 const Field<Type>& field
37)
38{
39 tmp<Field<Type>> tNewField(new Field<Type>(field.size()));
40 Field<Type>& newField = tNewField.ref();
41
42 label added = 0;
43 label count = 0;
44
45 for
46 (
48 mesh.finite_vertices_begin();
49 vit != mesh.finite_vertices_end();
50 ++vit
51 )
52 {
53 if (vit->real())
54 {
55 newField[added++] = field[count];
56 }
57
58 count++;
59 }
60
61 newField.resize(added);
62
63 return tNewField;
64}
65
66
67template<class Triangulation>
68Foam::autoPtr<Foam::mapDistribute> Foam::smoothAlignmentSolver::buildReferredMap
69(
70 const Triangulation& mesh,
71 labelList& indices
72)
73{
74 globalIndex globalIndexing(mesh.vertexCount());
75
76 DynamicList<label> dynIndices(mesh.vertexCount()/10);
77
78 for
79 (
81 mesh.finite_vertices_begin();
82 vit != mesh.finite_vertices_end();
83 ++vit
84 )
85 {
86 if (vit->referred())
87 {
88 dynIndices.append
89 (
90 globalIndexing.toGlobal(vit->procIndex(), vit->index())
91 );
92 }
93 }
94
95 indices.transfer(dynIndices);
96
97 List<Map<label>> compactMap;
99 (
100 globalIndexing,
101 indices,
102 compactMap
103 );
104}
105
106
107template<class Triangulation>
108Foam::autoPtr<Foam::mapDistribute> Foam::smoothAlignmentSolver::buildMap
109(
110 const Triangulation& mesh,
111 labelListList& pointPoints
112)
113{
114 pointPoints.setSize(mesh.vertexCount());
115
116 globalIndex globalIndexing(mesh.vertexCount());
117
118 for
119 (
121 mesh.finite_vertices_begin();
122 vit != mesh.finite_vertices_end();
123 ++vit
124 )
125 {
126 if (!vit->real())
127 {
128 continue;
129 }
130
131 std::list<typename Triangulation::Vertex_handle> adjVerts;
132 mesh.finite_adjacent_vertices(vit, std::back_inserter(adjVerts));
133
134 DynamicList<label> indices(adjVerts.size());
135
136 for
137 (
138 typename std::list<typename Triangulation::Vertex_handle>::
139 const_iterator adjVertI = adjVerts.begin();
140 adjVertI != adjVerts.end();
141 ++adjVertI
142 )
143 {
144 typename Triangulation::Vertex_handle vh = *adjVertI;
145
146 if (!vh->farPoint())
147 {
148 indices.append
149 (
150 globalIndexing.toGlobal(vh->procIndex(), vh->index())
151 );
152 }
153 }
154
155 pointPoints[vit->index()].transfer(indices);
156 }
157
158 List<Map<label>> compactMap;
160 (
161 globalIndexing,
162 pointPoints,
163 compactMap
164 );
165}
166
167
168template<class Triangulation>
169Foam::tmp<Foam::triadField> Foam::smoothAlignmentSolver::buildAlignmentField
170(
171 const Triangulation& mesh
172)
173{
174 tmp<triadField> tAlignments
175 (
176 new triadField(mesh.vertexCount(), triad::unset)
177 );
178 triadField& alignments = tAlignments.ref();
179
180 for
181 (
183 mesh.finite_vertices_begin();
184 vit != mesh.finite_vertices_end();
185 ++vit
186 )
187 {
188 if (!vit->real())
189 {
190 continue;
191 }
192
193 alignments[vit->index()] = vit->alignment();
194 }
195
196 return tAlignments;
197}
198
199
200template<class Triangulation>
201Foam::tmp<Foam::pointField> Foam::smoothAlignmentSolver::buildPointField
202(
203 const Triangulation& mesh
204)
205{
206 tmp<pointField> tPoints
207 (
208 new pointField(mesh.vertexCount(), point(GREAT, GREAT, GREAT))
209 );
210 pointField& points = tPoints.ref();
211
212 for
213 (
215 mesh.finite_vertices_begin();
216 vit != mesh.finite_vertices_end();
217 ++vit
218 )
219 {
220 if (!vit->real())
221 {
222 continue;
223 }
224
225 points[vit->index()] = topoint(vit->point());
226 }
227
228 return tPoints;
229}
230
231
232void Foam::smoothAlignmentSolver::applyBoundaryConditions
233(
234 const triad& fixedAlignment,
235 triad& t
236) const
237{
238 label nFixed = 0;
239
240 forAll(fixedAlignment, dirI)
241 {
242 if (fixedAlignment.set(dirI))
243 {
244 nFixed++;
245 }
246 }
247
248 if (nFixed == 1)
249 {
250 forAll(fixedAlignment, dirI)
251 {
252 if (fixedAlignment.set(dirI))
253 {
254 t.align(fixedAlignment[dirI]);
255 }
256 }
257 }
258 else if (nFixed == 2)
259 {
260 forAll(fixedAlignment, dirI)
261 {
262 if (fixedAlignment.set(dirI))
263 {
264 t[dirI] = fixedAlignment[dirI];
265 }
266 else
267 {
268 t[dirI] = triad::unset[dirI];
269 }
270 }
271
272 t.orthogonalize();
273 }
274 else if (nFixed == 3)
275 {
276 forAll(fixedAlignment, dirI)
277 {
278 if (fixedAlignment.set(dirI))
279 {
280 t[dirI] = fixedAlignment[dirI];
281 }
282 }
283 }
284}
285
286
287// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
288
290:
291 mesh_(mesh)
292{}
293
294
295// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
296
298{}
299
300
301// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
302
304(
305 const label maxSmoothingIterations
306)
307{
308 scalar minResidual = 0;
309
310 labelListList pointPoints;
311 autoPtr<mapDistribute> meshDistributor = buildMap
312 (
313 mesh_,
314 pointPoints
315 );
316
317 triadField alignments(buildAlignmentField(mesh_));
318 pointField points(buildPointField(mesh_));
319
320 // Setup the sizes and alignments on each point
321 triadField fixedAlignments(mesh_.vertexCount(), triad::unset);
322
323 for
324 (
325 CellSizeDelaunay::Finite_vertices_iterator vit =
326 mesh_.finite_vertices_begin();
327 vit != mesh_.finite_vertices_end();
328 ++vit
329 )
330 {
331 if (vit->real())
332 {
333 fixedAlignments[vit->index()] = vit->alignment();
334 }
335 }
336
337 Info<< nl << "Smoothing alignments" << endl;
338
339
340 for (label iter = 0; iter < maxSmoothingIterations; iter++)
341 {
342 Info<< "Iteration " << iter;
343
344 meshDistributor().distribute(points);
345 meshDistributor().distribute(fixedAlignments);
346 meshDistributor().distribute(alignments);
347
348 scalar residual = 0;
349
350 triadField triadAv(alignments.size(), triad::unset);
351
352 forAll(pointPoints, pI)
353 {
354 const labelList& pPoints = pointPoints[pI];
355
356 if (pPoints.empty())
357 {
358 continue;
359 }
360
361 triad& newTriad = triadAv[pI];
362
363 forAll(pPoints, adjPointi)
364 {
365 const label adjPointIndex = pPoints[adjPointi];
366
367 scalar dist = mag(points[pI] - points[adjPointIndex]);
368
369 triad tmpTriad = alignments[adjPointIndex];
370
371 for (direction dir = 0; dir < 3; dir++)
372 {
373 if (tmpTriad.set(dir))
374 {
375 tmpTriad[dir] *= 1.0/(dist + SMALL);
376 }
377 }
378
379 newTriad += tmpTriad;
380 }
381 }
382
383 // Update the alignment field
384 forAll(alignments, pI)
385 {
386 const triad& oldTriad = alignments[pI];
387 triad& newTriad = triadAv[pI];
388
389 newTriad.normalize();
390 newTriad.orthogonalize();
391
392 // Enforce the boundary conditions
393 const triad& fixedAlignment = fixedAlignments[pI];
394
395 applyBoundaryConditions
396 (
397 fixedAlignment,
398 newTriad
399 );
400
401 newTriad = newTriad.sortxyz();
402
403 // Residual Calculation
404 for (direction dir = 0; dir < 3; ++dir)
405 {
406 if
407 (
408 newTriad.set(dir)
409 && oldTriad.set(dir)
410 && !fixedAlignment.set(dir)
411 )
412 {
413 residual += diff(oldTriad, newTriad);
414 }
415 }
416
417 alignments[pI] = newTriad;
418 }
419
420 reduce(residual, sumOp<scalar>());
421
422 Info<< ", Residual = "
423 << residual
424 /returnReduce(points.size(), sumOp<label>())
425 << endl;
426
427 if (iter > 0 && residual <= minResidual)
428 {
429 break;
430 }
431 }
432
433 meshDistributor().distribute(alignments);
434
435 for
436 (
437 CellSizeDelaunay::Finite_vertices_iterator vit =
438 mesh_.finite_vertices_begin();
439 vit != mesh_.finite_vertices_end();
440 ++vit
441 )
442 {
443 if (vit->real())
444 {
445 vit->alignment() = alignments[vit->index()];
446 }
447 }
448
449 labelList referredPoints;
450 autoPtr<mapDistribute> referredDistributor = buildReferredMap
451 (
452 mesh_,
453 referredPoints
454 );
455
456 alignments.setSize(mesh_.vertexCount());
457 referredDistributor().distribute(alignments);
458
459 label referredI = 0;
460 for
461 (
462 CellSizeDelaunay::Finite_vertices_iterator vit =
463 mesh_.finite_vertices_begin();
464 vit != mesh_.finite_vertices_end();
465 ++vit
466 )
467 {
468 if (vit->referred())
469 {
470 vit->alignment() = alignments[referredPoints[referredI++]];
471 }
472 }
473}
474
475
476// ************************************************************************* //
Triangulation::Finite_vertices_iterator Finite_vertices_iterator
Definition: DelaunayMesh.H:73
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
CellSizeDelaunay::Vertex_handle Vertex_handle
void smoothAlignments(const label maxSmoothingIterations)
Smooth the alignments on the mesh.
~smoothAlignmentSolver()
Destructor.
A class for managing temporary objects.
Definition: tmp.H:65
static const triad unset
Definition: triad.H:97
void orthogonalize()
Same as orthogonalise.
Definition: triad.H:167
rDeltaTY field()
dynamicFvMesh & mesh
const pointField & points
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
Definition: BitOps.H:78
List< label > labelList
A List of labels.
Definition: List.H:66
pointFromPoint topoint(const Point &P)
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
vector point
Point is a vector.
Definition: point.H:43
Field< triad > triadField
Specialisation of Field<T> for triad.
Definition: triadFieldFwd.H:45
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
Definition: triad.C:378
uint8_t direction
Definition: direction.H:56
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333