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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2012-2017 OpenFOAM Foundation
9 Copyright (C) 2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
31#include "cellSizeFunction.H"
32#include "triSurfaceMesh.H"
33#include "searchableBox.H"
34#include "tetPointRef.H"
35#include "vectorTools.H"
36#include "quaternion.H"
37
38// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39
40namespace Foam
41{
42
43defineTypeNameAndDebug(searchableSurfaceControl, 0);
45(
46 cellSizeAndAlignmentControl,
47 searchableSurfaceControl,
48 dictionary
49);
50
51}
52
53
54// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
55
56//Foam::tensor Foam::surfaceControl::requiredAlignment
57//(
58// const Foam::point& pt,
59// const vectorField& ptNormals
60//) const
61//{
95//
101//
102// vector np = ptNormals[0];
103//
104// const tensor Rp = rotationTensor(vector(0,0,1), np);
105//
106// vector na = Zero;
107//
108// scalar smallestAngle = GREAT;
109//
110// for (label pnI = 1; pnI < ptNormals.size(); ++pnI)
111// {
112// const vector& nextNormal = ptNormals[pnI];
113//
114// const scalar cosPhi = vectorTools::cosPhi(np, nextNormal);
115//
116// if (mag(cosPhi) < smallestAngle)
117// {
118// na = nextNormal;
119// smallestAngle = cosPhi;
120// }
121// }
122//
123// // Secondary alignment
124// vector ns = np ^ na;
125//
126// if (mag(ns) < SMALL)
127// {
128// WarningInFunction
129// << "Parallel normals detected in spoke search." << nl
130// << "point: " << pt << nl
131// << "np : " << np << nl
132// << "na : " << na << nl
133// << "ns : " << ns << nl
134// << endl;
135//
136// ns = np;
137// }
138//
139// ns /= mag(ns);
140//
141// tensor Rs = rotationTensor((Rp & vector(0,1,0)), ns);
142//
149//
150// return (Rs & Rp);
151//}
152
153
154// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
155
157(
158 const Time& runTime,
159 const word& name,
160 const dictionary& controlFunctionDict,
161 const conformationSurfaces& geometryToConformTo,
162 const scalar& defaultCellSize
163)
164:
165 cellSizeAndAlignmentControl
166 (
167 runTime,
168 name,
169 controlFunctionDict,
170 geometryToConformTo,
171 defaultCellSize
172 ),
173 surfaceName_(controlFunctionDict.getOrDefault<word>("surface", name)),
174 searchableSurface_(geometryToConformTo.geometry()[surfaceName_]),
175 geometryToConformTo_(geometryToConformTo),
176 cellSizeFunctions_(1),
177 regionToCellSizeFunctions_(searchableSurface_.regions().size(), -1),
178 maxPriority_(-1)
179{
180 Info<< indent << "Master settings:" << endl;
182
183 cellSizeFunctions_.set
184 (
185 0,
186 cellSizeFunction::New
187 (
188 controlFunctionDict,
189 searchableSurface_,
190 defaultCellSize_,
191 labelList()
192 )
193 );
194
196
197 PtrList<cellSizeFunction> regionCellSizeFunctions;
198
199 DynamicList<label> defaultCellSizeRegions;
200
201 label nRegionCellSizeFunctions = 0;
202
203 // Loop over regions - if any entry is not specified they should
204 // inherit values from the parent surface.
205 if (controlFunctionDict.found("regions"))
206 {
207 const dictionary& regionsDict = controlFunctionDict.subDict("regions");
208 const wordList& regionNames = searchableSurface_.regions();
209
210 label nRegions = regionsDict.size();
211
212 regionCellSizeFunctions.setSize(nRegions);
213 defaultCellSizeRegions.setCapacity(nRegions);
214
215 forAll(regionNames, regionI)
216 {
217 const word& regionName = regionNames[regionI];
218
219 label regionID = geometryToConformTo_.geometry().findSurfaceRegionID
220 (
221 this->name(),
223 );
224
225 if (regionsDict.found(regionName))
226 {
227 // Get the dictionary for region
228 const dictionary& regionDict = regionsDict.subDict(regionName);
229
230 Info<< indent << "Region " << regionName
231 << " (ID = " << regionID << ")" << " settings:"
232 << endl;
234
235 regionCellSizeFunctions.set
236 (
237 nRegionCellSizeFunctions,
238 cellSizeFunction::New
239 (
240 regionDict,
241 searchableSurface_,
242 defaultCellSize_,
243 labelList(1, regionID)
244 )
245 );
247
248 regionToCellSizeFunctions_[regionID] = nRegionCellSizeFunctions;
249
250 nRegionCellSizeFunctions++;
251 }
252 else
253 {
254 // Add to default list
255 defaultCellSizeRegions.append(regionID);
256 }
257 }
258 }
259
260 if (defaultCellSizeRegions.empty() && !regionCellSizeFunctions.empty())
261 {
262 cellSizeFunctions_.transfer(regionCellSizeFunctions);
263 }
264 else if (nRegionCellSizeFunctions > 0)
265 {
266 regionCellSizeFunctions.setSize(nRegionCellSizeFunctions + 1);
267
268 regionCellSizeFunctions.set
269 (
270 nRegionCellSizeFunctions,
271 cellSizeFunction::New
272 (
273 controlFunctionDict,
274 searchableSurface_,
275 defaultCellSize_,
276 labelList()
277 )
278 );
279
280 const wordList& regionNames = searchableSurface_.regions();
281
282 forAll(regionNames, regionI)
283 {
284 if (regionToCellSizeFunctions_[regionI] == -1)
285 {
286 regionToCellSizeFunctions_[regionI] = nRegionCellSizeFunctions;
287 }
288 }
289
290 cellSizeFunctions_.transfer(regionCellSizeFunctions);
291 }
292 else
293 {
294 const wordList& regionNames = searchableSurface_.regions();
295
296 forAll(regionNames, regionI)
297 {
298 if (regionToCellSizeFunctions_[regionI] == -1)
299 {
300 regionToCellSizeFunctions_[regionI] = 0;
301 }
302 }
303 }
304
305
306 forAll(cellSizeFunctions_, funcI)
307 {
308 const label funcPriority = cellSizeFunctions_[funcI].priority();
309
310 if (funcPriority > maxPriority_)
311 {
312 maxPriority_ = funcPriority;
313 }
314 }
315
316 // Sort controlFunctions_ by maxPriority
317 SortableList<label> functionPriorities(cellSizeFunctions_.size());
318
319 forAll(cellSizeFunctions_, funcI)
320 {
321 functionPriorities[funcI] = cellSizeFunctions_[funcI].priority();
322 }
323
324 functionPriorities.reverseSort();
325
326 labelList invertedFunctionPriorities =
327 invert(functionPriorities.size(), functionPriorities.indices());
328
329 cellSizeFunctions_.reorder(invertedFunctionPriorities);
330
331 Info<< nl << indent << "There are " << cellSizeFunctions_.size()
332 << " region control functions" << endl;
333}
334
335
336// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
337
339{}
340
341
342// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
343
345(
346 pointField& pts,
347 scalarField& sizes,
348 triadField& alignments
349) const
350{
351 pts = searchableSurface_.points();
352 sizes.setSize(pts.size());
353 alignments.setSize(pts.size());
354
355 const scalar nearFeatDistSqrCoeff = 1e-8;
356
357 forAll(pts, pI)
358 {
359 // Is the point in the extendedFeatureEdgeMesh? If so get the
360 // point normal, otherwise get the surface normal from
361 // searchableSurface
362
363 pointIndexHit info;
364 label infoFeature;
365 geometryToConformTo_.findFeaturePointNearest
366 (
367 pts[pI],
368 nearFeatDistSqrCoeff,
369 info,
370 infoFeature
371 );
372
373 scalar limitedCellSize = GREAT;
374
375 autoPtr<triad> pointAlignment;
376
377 if (info.hit())
378 {
379 const extendedFeatureEdgeMesh& features =
380 geometryToConformTo_.features()[infoFeature];
381
382 vectorField norms = features.featurePointNormals(info.index());
383
384 // Create a triad from these norms.
385 pointAlignment.reset(new triad());
386 forAll(norms, nI)
387 {
388 pointAlignment() += norms[nI];
389 }
390
391 pointAlignment().normalize();
392 pointAlignment().orthogonalize();
393 }
394 else
395 {
396 geometryToConformTo_.findEdgeNearest
397 (
398 pts[pI],
399 nearFeatDistSqrCoeff,
400 info,
401 infoFeature
402 );
403
404 if (info.hit())
405 {
406 const extendedFeatureEdgeMesh& features =
407 geometryToConformTo_.features()[infoFeature];
408
409 vectorField norms = features.edgeNormals(info.index());
410
411 // Create a triad from these norms.
412 pointAlignment.reset(new triad());
413 forAll(norms, nI)
414 {
415 pointAlignment() += norms[nI];
416 }
417
418 pointAlignment().normalize();
419 pointAlignment().orthogonalize();
420 }
421 else
422 {
423 pointField ptField(1, pts[pI]);
424 scalarField distField(1, nearFeatDistSqrCoeff);
425 List<pointIndexHit> infoList(1, pointIndexHit());
426
427 searchableSurface_.findNearest(ptField, distField, infoList);
428
429 vectorField normals(1);
430 searchableSurface_.getNormal(infoList, normals);
431
432 if (mag(normals[0]) < SMALL)
433 {
434 normals[0] = vector::one;
435 }
436
437 pointAlignment.reset(new triad(normals[0]));
438
439 if (infoList[0].hit())
440 {
441 // Limit cell size
442 const vector vN =
443 infoList[0].hitPoint()
444 - 2.0*normals[0]*defaultCellSize_;
445
446 List<pointIndexHit> intersectionList;
447 searchableSurface_.findLineAny
448 (
449 ptField,
450 pointField(1, vN),
451 intersectionList
452 );
453 }
454
455// if (intersectionList[0].hit())
456// {
457// scalar dist =
458// mag(intersectionList[0].hitPoint() - pts[pI]);
459//
460// limitedCellSize = dist/2.0;
461// }
462 }
463 }
464
465 label priority = -1;
466 if (!cellSize(pts[pI], sizes[pI], priority))
467 {
468 sizes[pI] = defaultCellSize_;
469// FatalErrorInFunction
470// << "Could not calculate cell size"
471// << abort(FatalError);
472 }
473
474 sizes[pI] = min(limitedCellSize, sizes[pI]);
475
476 alignments[pI] = pointAlignment();
477 }
478}
479
480
482(
483 DynamicList<Foam::point>& pts,
484 DynamicList<scalar>& sizes
485) const
486{
487 const tmp<pointField> tmpPoints = searchableSurface_.points();
488 const pointField& points = tmpPoints();
489
490 const scalar nearFeatDistSqrCoeff = 1e-8;
491
492 pointField ptField(1, vector::min);
493 scalarField distField(1, nearFeatDistSqrCoeff);
494 List<pointIndexHit> infoList(1, pointIndexHit());
495
496 vectorField normals(1);
497 labelList region(1, label(-1));
498
499 forAll(points, pI)
500 {
501 // Is the point in the extendedFeatureEdgeMesh? If so get the
502 // point normal, otherwise get the surface normal from
503 // searchableSurface
504 ptField[0] = points[pI];
505
506 searchableSurface_.findNearest(ptField, distField, infoList);
507
508 if (infoList[0].hit())
509 {
510 searchableSurface_.getNormal(infoList, normals);
511 searchableSurface_.getRegion(infoList, region);
512
513 const cellSizeFunction& sizeFunc =
514 sizeFunctions()[regionToCellSizeFunctions_[region[0]]];
515
516 pointField extraPts;
517 scalarField extraSizes;
518 sizeFunc.sizeLocations
519 (
520 infoList[0],
521 normals[0],
522 extraPts,
523 extraSizes
524 );
525
526 pts.append(extraPts);
527 sizes.append(extraSizes);
528 }
529 }
530}
531
532
534(
535 const Foam::point& pt,
536 scalar& cellSize,
537 label& priority
538) const
539{
540 bool anyFunctionFound = false;
541
542 forAll(sizeFunctions(), funcI)
543 {
544 const cellSizeFunction& sizeFunc = sizeFunctions()[funcI];
545
546 if (sizeFunc.priority() < priority)
547 {
548 continue;
549 }
550
551 scalar sizeI = -1;
552
553 if (sizeFunc.cellSize(pt, sizeI))
554 {
555 anyFunctionFound = true;
556
557 if (sizeFunc.priority() == priority)
558 {
559 if (sizeI < cellSize)
560 {
561 cellSize = sizeI;
562 }
563 }
564 else
565 {
566 cellSize = sizeI;
567
568 priority = sizeFunc.priority();
569 }
570
571 if (debug)
572 {
573 Info<< " sizeI " << sizeI
574 <<" minSize " << cellSize << endl;
575 }
576 }
577 }
578
579 return anyFunctionFound;
580}
581
582
583// ************************************************************************* //
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
static const complex min
complex (-VGREAT,-VGREAT)
Definition: complex.H:292
virtual void cellSizeFunctionVertices(DynamicList< Foam::point > &pts, DynamicList< scalar > &sizes) const
~searchableSurfaceControl()
Destructor.
bool cellSize(const Foam::point &pt, scalar &cellSize, label &priority) const
virtual void initialVertices(pointField &pts, scalarField &sizes, triadField &alignments) const
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
engineTime & runTime
Foam::word regionName(Foam::polyMesh::defaultRegion)
wordList regionNames
const pointField & points
Namespace for OpenFOAM.
List< word > wordList
A List of words.
Definition: fileName.H:63
List< label > labelList
A List of labels.
Definition: List.H:66
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Field< triad > triadField
Specialisation of Field<T> for triad.
Definition: triadFieldFwd.H:45
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:349
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:342
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Field< vector > vectorField
Specialisation of Field<T> for vector.
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
Definition: ListOps.C:36
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:356
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333