cellShapeControl.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-2017 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
28#include "cellShapeControl.H"
29#include "pointField.H"
30#include "scalarField.H"
31#include "triadField.H"
34#include "cellSizeFunction.H"
35#include "indexedVertexOps.H"
36
37// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38
39namespace Foam
40{
41 defineTypeNameAndDebug(cellShapeControl, 0);
42}
43
44
45// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
46
48(
49 const Time& runTime,
50 const cvControls& foamyHexMeshControls,
51 const searchableSurfaces& allGeometry,
52 const conformationSurfaces& geometryToConformTo
53)
54:
55 dictionary
56 (
57 foamyHexMeshControls.foamyHexMeshDict().subDict("motionControl")
58 ),
59 geometryToConformTo_(geometryToConformTo),
60 defaultCellSize_(foamyHexMeshControls.defaultCellSize()),
61 minimumCellSize_(foamyHexMeshControls.minimumCellSize()),
62 shapeControlMesh_(runTime),
63 aspectRatio_(*this),
64 sizeAndAlignment_
65 (
66 runTime,
67 subDict("shapeControlFunctions"),
68 geometryToConformTo_,
69 defaultCellSize_
70 )
71{}
72
73
74// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
75
77(
78 const pointField& pts
79) const
80{
81 scalarField cellSizes(pts.size());
82
83 forAll(pts, i)
84 {
85 cellSizes[i] = cellSize(pts[i]);
86 }
87
88 return cellSizes;
89}
90
91
92Foam::scalar Foam::cellShapeControl::cellSize(const point& pt) const
93{
94 barycentric bary;
96
97 shapeControlMesh_.barycentricCoords(pt, bary, ch);
98
99 scalar size = 0;
100
101 if (shapeControlMesh_.dimension() < 3)
102 {
103 size = sizeAndAlignment_.cellSize(pt);
104 }
105 else if (shapeControlMesh_.is_infinite(ch))
106 {
107// if (nFarPoints != 0)
108// {
109// for (label pI = 0; pI < 4; ++pI)
110// {
111// if (!ch->vertex(pI)->farPoint())
112// {
113// size = ch->vertex(pI)->targetCellSize();
114// return size;
115// }
116// }
117// }
118
119// cellShapeControlMesh::Vertex_handle nearV =
120// shapeControlMesh_.nearest_vertex_in_cell
121// (
122// toPoint<cellShapeControlMesh::Point>(pt),
123// ch
124// );
125//
126// size = nearV->targetCellSize();
127
128 // Find nearest surface. This can be quite slow if there are a lot of
129 // surfaces
130 size = sizeAndAlignment_.cellSize(pt);
131 }
132 else
133 {
134 label nFarPoints = 0;
135 for (label pI = 0; pI < 4; ++pI)
136 {
137 if (ch->vertex(pI)->farPoint())
138 {
139 nFarPoints++;
140 }
141 }
142
143 if (nFarPoints != 0)
144 {
145 for (label pI = 0; pI < 4; ++pI)
146 {
147 if (!CGAL::indexedVertexOps::uninitialised(ch->vertex(pI)))
148 {
149 size = ch->vertex(pI)->targetCellSize();
150 return size;
151 }
152 }
153 }
154 else
155 {
156 forAll(bary, pI)
157 {
158 size += bary[pI]*ch->vertex(pI)->targetCellSize();
159 }
160 }
161 }
162
163 return size;
164}
165
166
168{
169 barycentric bary;
171
172 shapeControlMesh_.barycentricCoords(pt, bary, ch);
173
174 tensor alignment = Zero;
175
176 if (shapeControlMesh_.dimension() < 3 || shapeControlMesh_.is_infinite(ch))
177 {
178 alignment = tensor::I;
179 }
180 else
181 {
182 label nFarPoints = 0;
183 for (label pI = 0; pI < 4; ++pI)
184 {
185 if (ch->vertex(pI)->farPoint())
186 {
187 nFarPoints++;
188 }
189 }
190
191// if (nFarPoints != 0)
192// {
193// for (label pI = 0; pI < 4; ++pI)
194// {
195// if (!ch->vertex(pI)->farPoint())
196// {
197// alignment = ch->vertex(pI)->alignment();
198// }
199// }
200// }
201// else
202 {
203 triad tri;
204
205 for (label pI = 0; pI < 4; ++pI)
206 {
207 if (bary[pI] > SMALL)
208 {
209 tri += triad(bary[pI]*ch->vertex(pI)->alignment());
210 }
211 }
212
213 tri.normalize();
214 tri.orthogonalize();
215 tri = tri.sortxyz();
216
217 alignment = tri;
218 }
219
220// cellShapeControlMesh::Vertex_handle nearV =
221// shapeControlMesh_.nearest_vertex_in_cell
222// (
223// toPoint<cellShapeControlMesh::Point>(pt),
224// ch
225// );
226//
227// alignment = nearV->alignment();
228 }
229
230 return alignment;
231}
232
233
235(
236 const point& pt,
237 scalar& size,
238 tensor& alignment
239) const
240{
241 barycentric bary;
243
244 shapeControlMesh_.barycentricCoords(pt, bary, ch);
245
246 alignment = Zero;
247 size = 0;
248
249 if (shapeControlMesh_.dimension() < 3 || shapeControlMesh_.is_infinite(ch))
250 {
251 // Find nearest surface
252 size = sizeAndAlignment_.cellSize(pt);
253 alignment = tensor::I;
254 }
255 else
256 {
257 label nFarPoints = 0;
258 for (label pI = 0; pI < 4; ++pI)
259 {
260 if (ch->vertex(pI)->farPoint())
261 {
262 nFarPoints++;
263 }
264 }
265
266 if (nFarPoints != 0)
267 {
268 for (label pI = 0; pI < 4; ++pI)
269 {
270 if (!CGAL::indexedVertexOps::uninitialised(ch->vertex(pI)))
271 {
272 size = ch->vertex(pI)->targetCellSize();
273 alignment = ch->vertex(pI)->alignment();
274 }
275 }
276 }
277 else
278 {
279 triad tri;
280
281 for (label pI = 0; pI < 4; ++pI)
282 {
283 size += bary[pI]*ch->vertex(pI)->targetCellSize();
284
285 if (bary[pI] > SMALL)
286 {
287 tri += triad(bary[pI]*ch->vertex(pI)->alignment());
288 }
289 }
290
291 tri.normalize();
292 tri.orthogonalize();
293 tri = tri.sortxyz();
294
295 alignment = tri;
296
297// cellShapeControlMesh::Vertex_handle nearV =
298// shapeControlMesh_.nearest_vertex
299// (
300// toPoint<cellShapeControlMesh::Point>(pt)
301// );
302//
303// alignment = nearV->alignment();
304 }
305 }
306
307 for (label dir = 0; dir < 3; dir++)
308 {
309 triad v = alignment;
310
311 if (!v.set(dir) || size == 0)
312 {
313 // Force orthogonalization of triad.
314
315 scalar dotProd = GREAT;
316 if (dir == 0)
317 {
318 dotProd = v[1] & v[2];
319
320 v[dir] = v[1] ^ v[2];
321 }
322 if (dir == 1)
323 {
324 dotProd = v[0] & v[2];
325
326 v[dir] = v[0] ^ v[2];
327 }
328 if (dir == 2)
329 {
330 dotProd = v[0] & v[1];
331
332 v[dir] = v[0] ^ v[1];
333 }
334
335 v.normalize();
336 v.orthogonalize();
337
338 Pout<< "Dot prod = " << dotProd << endl;
339 Pout<< "Alignment = " << v << endl;
340
341 alignment = v;
342
343// FatalErrorInFunction
344// << "Point has bad alignment! "
345// << pt << " " << size << " " << alignment << nl
346// << "Bary Coords = " << bary << nl
347// << ch->vertex(0)->info() << nl
348// << ch->vertex(1)->info() << nl
349// << ch->vertex(2)->info() << nl
350// << ch->vertex(3)->info()
351// << abort(FatalError);
352 }
353 }
354}
355
356
357// ************************************************************************* //
static const Tensor I
Definition: Tensor.H:81
CellSizeDelaunay::Cell_handle Cell_handle
tensor cellAlignment(const point &pt) const
Return the cell alignment at the given location.
void cellSizeAndAlignment(const point &pt, scalar &size, tensor &alignment) const
scalar cellSize(const point &pt) const
Return the cell size at the given location.
Tensor of scalars, i.e. Tensor<scalar>.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
engineTime & runTime
bool uninitialised(const VertexType &v)
Namespace for OpenFOAM.
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
Definition: barycentric.H:51
vector point
Point is a vector.
Definition: point.H:43
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333