searchableSurfaceControl.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 
28 #include "cellSizeFunction.H"
29 #include "triSurfaceMesh.H"
30 #include "searchableBox.H"
31 #include "tetrahedron.H"
32 #include "vectorTools.H"
33 #include "quaternion.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 
40 defineTypeNameAndDebug(searchableSurfaceControl, 0);
42 (
43  cellSizeAndAlignmentControl,
44  searchableSurfaceControl,
45  dictionary
46 );
47 
48 }
49 
50 
51 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
52 
53 //Foam::tensor Foam::surfaceControl::requiredAlignment
54 //(
55 // const Foam::point& pt,
56 // const vectorField& ptNormals
57 //) const
58 //{
92 //
98 //
99 // vector np = ptNormals[0];
100 //
101 // const tensor Rp = rotationTensor(vector(0,0,1), np);
102 //
103 // vector na = Zero;
104 //
105 // scalar smallestAngle = GREAT;
106 //
107 // for (label pnI = 1; pnI < ptNormals.size(); ++pnI)
108 // {
109 // const vector& nextNormal = ptNormals[pnI];
110 //
111 // const scalar cosPhi = vectorTools::cosPhi(np, nextNormal);
112 //
113 // if (mag(cosPhi) < smallestAngle)
114 // {
115 // na = nextNormal;
116 // smallestAngle = cosPhi;
117 // }
118 // }
119 //
120 // // Secondary alignment
121 // vector ns = np ^ na;
122 //
123 // if (mag(ns) < SMALL)
124 // {
125 // WarningInFunction
126 // << "Parallel normals detected in spoke search." << nl
127 // << "point: " << pt << nl
128 // << "np : " << np << nl
129 // << "na : " << na << nl
130 // << "ns : " << ns << nl
131 // << endl;
132 //
133 // ns = np;
134 // }
135 //
136 // ns /= mag(ns);
137 //
138 // tensor Rs = rotationTensor((Rp & vector(0,1,0)), ns);
139 //
146 //
147 // return (Rs & Rp);
148 //}
149 
150 
151 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
152 
153 Foam::searchableSurfaceControl::searchableSurfaceControl
154 (
155  const Time& runTime,
156  const word& name,
157  const dictionary& controlFunctionDict,
158  const conformationSurfaces& geometryToConformTo,
159  const scalar& defaultCellSize
160 )
161 :
162  cellSizeAndAlignmentControl
163  (
164  runTime,
165  name,
166  controlFunctionDict,
167  geometryToConformTo,
168  defaultCellSize
169  ),
170  surfaceName_(controlFunctionDict.lookupOrDefault<word>("surface", name)),
171  searchableSurface_(geometryToConformTo.geometry()[surfaceName_]),
172  geometryToConformTo_(geometryToConformTo),
173  cellSizeFunctions_(1),
174  regionToCellSizeFunctions_(searchableSurface_.regions().size(), -1),
175  maxPriority_(-1)
176 {
177  Info<< indent << "Master settings:" << endl;
178  Info<< incrIndent;
179 
180  cellSizeFunctions_.set
181  (
182  0,
184  (
185  controlFunctionDict,
186  searchableSurface_,
188  labelList()
189  )
190  );
191 
192  Info<< decrIndent;
193 
194  PtrList<cellSizeFunction> regionCellSizeFunctions;
195 
196  DynamicList<label> defaultCellSizeRegions;
197 
198  label nRegionCellSizeFunctions = 0;
199 
200  // Loop over regions - if any entry is not specified they should
201  // inherit values from the parent surface.
202  if (controlFunctionDict.found("regions"))
203  {
204  const dictionary& regionsDict = controlFunctionDict.subDict("regions");
205  const wordList& regionNames = searchableSurface_.regions();
206 
207  label nRegions = regionsDict.size();
208 
209  regionCellSizeFunctions.setSize(nRegions);
210  defaultCellSizeRegions.setCapacity(nRegions);
211 
212  forAll(regionNames, regionI)
213  {
214  const word& regionName = regionNames[regionI];
215 
216  label regionID = geometryToConformTo_.geometry().findSurfaceRegionID
217  (
218  this->name(),
219  regionName
220  );
221 
222  if (regionsDict.found(regionName))
223  {
224  // Get the dictionary for region
225  const dictionary& regionDict = regionsDict.subDict(regionName);
226 
227  Info<< indent << "Region " << regionName
228  << " (ID = " << regionID << ")" << " settings:"
229  << endl;
230  Info<< incrIndent;
231 
232  regionCellSizeFunctions.set
233  (
234  nRegionCellSizeFunctions,
236  (
237  regionDict,
238  searchableSurface_,
240  labelList(1, regionID)
241  )
242  );
243  Info<< decrIndent;
244 
245  regionToCellSizeFunctions_[regionID] = nRegionCellSizeFunctions;
246 
247  nRegionCellSizeFunctions++;
248  }
249  else
250  {
251  // Add to default list
252  defaultCellSizeRegions.append(regionID);
253  }
254  }
255  }
256 
257  if (defaultCellSizeRegions.empty() && !regionCellSizeFunctions.empty())
258  {
259  cellSizeFunctions_.transfer(regionCellSizeFunctions);
260  }
261  else if (nRegionCellSizeFunctions > 0)
262  {
263  regionCellSizeFunctions.setSize(nRegionCellSizeFunctions + 1);
264 
265  regionCellSizeFunctions.set
266  (
267  nRegionCellSizeFunctions,
269  (
270  controlFunctionDict,
271  searchableSurface_,
273  labelList()
274  )
275  );
276 
277  const wordList& regionNames = searchableSurface_.regions();
278 
279  forAll(regionNames, regionI)
280  {
281  if (regionToCellSizeFunctions_[regionI] == -1)
282  {
283  regionToCellSizeFunctions_[regionI] = nRegionCellSizeFunctions;
284  }
285  }
286 
287  cellSizeFunctions_.transfer(regionCellSizeFunctions);
288  }
289  else
290  {
291  const wordList& regionNames = searchableSurface_.regions();
292 
293  forAll(regionNames, regionI)
294  {
295  if (regionToCellSizeFunctions_[regionI] == -1)
296  {
297  regionToCellSizeFunctions_[regionI] = 0;
298  }
299  }
300  }
301 
302 
303  forAll(cellSizeFunctions_, funcI)
304  {
305  const label funcPriority = cellSizeFunctions_[funcI].priority();
306 
307  if (funcPriority > maxPriority_)
308  {
309  maxPriority_ = funcPriority;
310  }
311  }
312 
313  // Sort controlFunctions_ by maxPriority
314  SortableList<label> functionPriorities(cellSizeFunctions_.size());
315 
316  forAll(cellSizeFunctions_, funcI)
317  {
318  functionPriorities[funcI] = cellSizeFunctions_[funcI].priority();
319  }
320 
321  functionPriorities.reverseSort();
322 
323  labelList invertedFunctionPriorities =
324  invert(functionPriorities.size(), functionPriorities.indices());
325 
326  cellSizeFunctions_.reorder(invertedFunctionPriorities);
327 
328  Info<< nl << indent << "There are " << cellSizeFunctions_.size()
329  << " region control functions" << endl;
330 }
331 
332 
333 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
334 
336 {}
337 
338 
339 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
340 
342 (
343  pointField& pts,
344  scalarField& sizes,
345  triadField& alignments
346 ) const
347 {
348  pts = searchableSurface_.points();
349  sizes.setSize(pts.size());
350  alignments.setSize(pts.size());
351 
352  const scalar nearFeatDistSqrCoeff = 1e-8;
353 
354  forAll(pts, pI)
355  {
356  // Is the point in the extendedFeatureEdgeMesh? If so get the
357  // point normal, otherwise get the surface normal from
358  // searchableSurface
359 
360  pointIndexHit info;
361  label infoFeature;
362  geometryToConformTo_.findFeaturePointNearest
363  (
364  pts[pI],
365  nearFeatDistSqrCoeff,
366  info,
367  infoFeature
368  );
369 
370  scalar limitedCellSize = GREAT;
371 
372  autoPtr<triad> pointAlignment;
373 
374  if (info.hit())
375  {
376  const extendedFeatureEdgeMesh& features =
377  geometryToConformTo_.features()[infoFeature];
378 
379  vectorField norms = features.featurePointNormals(info.index());
380 
381  // Create a triad from these norms.
382  pointAlignment.set(new triad());
383  forAll(norms, nI)
384  {
385  pointAlignment() += norms[nI];
386  }
387 
388  pointAlignment().normalize();
389  pointAlignment().orthogonalize();
390  }
391  else
392  {
393  geometryToConformTo_.findEdgeNearest
394  (
395  pts[pI],
396  nearFeatDistSqrCoeff,
397  info,
398  infoFeature
399  );
400 
401  if (info.hit())
402  {
403  const extendedFeatureEdgeMesh& features =
404  geometryToConformTo_.features()[infoFeature];
405 
406  vectorField norms = features.edgeNormals(info.index());
407 
408  // Create a triad from these norms.
409  pointAlignment.set(new triad());
410  forAll(norms, nI)
411  {
412  pointAlignment() += norms[nI];
413  }
414 
415  pointAlignment().normalize();
416  pointAlignment().orthogonalize();
417  }
418  else
419  {
420  pointField ptField(1, pts[pI]);
421  scalarField distField(1, nearFeatDistSqrCoeff);
422  List<pointIndexHit> infoList(1, pointIndexHit());
423 
424  searchableSurface_.findNearest(ptField, distField, infoList);
425 
426  vectorField normals(1);
427  searchableSurface_.getNormal(infoList, normals);
428 
429  if (mag(normals[0]) < SMALL)
430  {
431  normals[0] = vector(1, 1, 1);
432  }
433 
434  pointAlignment.set(new triad(normals[0]));
435 
436  if (infoList[0].hit())
437  {
438  // Limit cell size
439  const vector vN =
440  infoList[0].hitPoint()
441  - 2.0*normals[0]*defaultCellSize_;
442 
443  List<pointIndexHit> intersectionList;
444  searchableSurface_.findLineAny
445  (
446  ptField,
447  pointField(1, vN),
448  intersectionList
449  );
450  }
451 
452 // if (intersectionList[0].hit())
453 // {
454 // scalar dist =
455 // mag(intersectionList[0].hitPoint() - pts[pI]);
456 //
457 // limitedCellSize = dist/2.0;
458 // }
459  }
460  }
461 
462  label priority = -1;
463  if (!cellSize(pts[pI], sizes[pI], priority))
464  {
465  sizes[pI] = defaultCellSize_;
466 // FatalErrorInFunction
467 // << "Could not calculate cell size"
468 // << abort(FatalError);
469  }
470 
471  sizes[pI] = min(limitedCellSize, sizes[pI]);
472 
473  alignments[pI] = pointAlignment();
474  }
475 }
476 
477 
479 (
480  DynamicList<Foam::point>& pts,
481  DynamicList<scalar>& sizes
482 ) const
483 {
484  const tmp<pointField> tmpPoints = searchableSurface_.points();
485  const pointField& points = tmpPoints();
486 
487  const scalar nearFeatDistSqrCoeff = 1e-8;
488 
489  pointField ptField(1, vector::min);
490  scalarField distField(1, nearFeatDistSqrCoeff);
491  List<pointIndexHit> infoList(1, pointIndexHit());
492 
493  vectorField normals(1);
494  labelList region(1, label(-1));
495 
496  forAll(points, pI)
497  {
498  // Is the point in the extendedFeatureEdgeMesh? If so get the
499  // point normal, otherwise get the surface normal from
500  // searchableSurface
501  ptField[0] = points[pI];
502 
503  searchableSurface_.findNearest(ptField, distField, infoList);
504 
505  if (infoList[0].hit())
506  {
507  searchableSurface_.getNormal(infoList, normals);
508  searchableSurface_.getRegion(infoList, region);
509 
510  const cellSizeFunction& sizeFunc =
511  sizeFunctions()[regionToCellSizeFunctions_[region[0]]];
512 
513  pointField extraPts;
514  scalarField extraSizes;
515  sizeFunc.sizeLocations
516  (
517  infoList[0],
518  normals[0],
519  extraPts,
520  extraSizes
521  );
522 
523  pts.append(extraPts);
524  sizes.append(extraSizes);
525  }
526  }
527 }
528 
529 
531 (
532  const Foam::point& pt,
533  scalar& cellSize,
534  label& priority
535 ) const
536 {
537  bool anyFunctionFound = false;
538 
539  forAll(sizeFunctions(), funcI)
540  {
541  const cellSizeFunction& sizeFunc = sizeFunctions()[funcI];
542 
543  if (sizeFunc.priority() < priority)
544  {
545  continue;
546  }
547 
548  scalar sizeI = -1;
549 
550  if (sizeFunc.cellSize(pt, sizeI))
551  {
552  anyFunctionFound = true;
553 
554  if (sizeFunc.priority() == priority)
555  {
556  if (sizeI < cellSize)
557  {
558  cellSize = sizeI;
559  }
560  }
561  else
562  {
563  cellSize = sizeI;
564 
565  priority = sizeFunc.priority();
566  }
567 
568  if (debug)
569  {
570  Info<< " sizeI " << sizeI
571  <<" minSize " << cellSize << endl;
572  }
573  }
574  }
575 
576  return anyFunctionFound;
577 }
578 
579 
580 // ************************************************************************* //
virtual void cellSizeFunctionVertices(DynamicList< Foam::point > &pts, DynamicList< scalar > &sizes) const
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:223
const double e
Elementary charge.
Definition: doubleFloat.H:78
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
virtual void initialVertices(pointField &pts, scalarField &sizes, triadField &alignments) const
Macros for easy insertion into run-time selection tables.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
Definition: ListOps.C:37
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:97
List< label > labelList
A List of labels.
Definition: labelList.H:56
Field< triad > triadField
Specialisation of Field<T> for triad.
Definition: triadField.H:49
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
static const char nl
Definition: Ostream.H:262
defineTypeNameAndDebug(combustionModel, 0)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:237
bool cellSize(const Foam::point &pt, scalar &cellSize, label &priority) const
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
List< word > wordList
A List of words.
Definition: fileName.H:54
~searchableSurfaceControl()
Destructor.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
static autoPtr< cellSizeFunction > New(const dictionary &cellSizeFunctionDict, const searchableSurface &surface, const scalar &defaultCellSize, const labelList regionIndices)
Return a reference to the selected cellSizeFunction.
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:230
Namespace for OpenFOAM.