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