45 Foam::scalar Foam::controlMeshRefinement::calcFirstDerivative
48 const scalar& cellSizeA,
50 const scalar& cellSizeB
53 return (cellSizeA - cellSizeB)/
mag(a -
b);
71 bool Foam::controlMeshRefinement::detectEdge
77 const scalar secondDerivTolSqr
96 pointFound.hitPoint(midPoint);
102 scalar cellSizeA = sizeControls_.cellSize(a);
103 scalar cellSizeB = sizeControls_.cellSize(
b);
110 scalar cellSizeMid = sizeControls_.cellSize(midPoint);
114 const scalar cellSizeMid1 = sizeControls_.cellSize(midPoint1);
119 scalar secondDerivative1 =
132 const scalar cellSizeMid2 = sizeControls_.cellSize(midPoint2);
137 scalar secondDerivative2 =
152 magSqr(secondDerivative1) < secondDerivTolSqr
153 &&
magSqr(secondDerivative2) < secondDerivTolSqr
160 if (
magSqr(secondDerivative1) >
magSqr(secondDerivative2))
163 midPoint = midPoint1;
168 midPoint = midPoint2;
181 const scalar tolSqr =
sqr(1
e-3);
182 const scalar secondDerivTolSqr =
sqr(1
e-3);
202 Foam::controlMeshRefinement::controlMeshRefinement
204 cellShapeControl& shapeController
207 shapeController_(shapeController),
208 mesh_(shapeController.shapeControlMesh()),
209 sizeControls_(shapeController.sizeAndAlignment()),
210 geometryToConformTo_(sizeControls_.geometryToConformTo())
224 const autoPtr<backgroundMeshDecomposition>& decomposition
227 if (shapeController_.shapeControlMesh().vertexCount() > 0)
230 Info<<
"Cell size and alignment mesh already populated." <<
endl;
234 autoPtr<boundBox> overallBoundBox;
256 Map<label> priorityMap;
258 const PtrList<cellSizeAndAlignmentControl>& controlFunctions =
259 sizeControls_.controlFunctions();
261 forAll(controlFunctions, fI)
263 const cellSizeAndAlignmentControl& controlFunction =
264 controlFunctions[fI];
266 const Switch forceInsertion =
267 controlFunction.forceInitialPointInsertion();
269 Info<<
"Inserting points from " << controlFunction.name()
270 <<
" (" << controlFunction.type() <<
")" <<
endl;
271 Info<<
" Force insertion is " << forceInsertion.c_str() <<
endl;
277 controlFunction.initialVertices(
pts, sizes, alignments);
284 for (label vI = 0; vI <
pts.
size(); ++vI)
288 label maxPriority = -1;
289 scalar size = sizeControls_.cellSize(
pts[vI], maxPriority);
291 if (maxPriority > controlFunction.maxPriority())
296 shapeController_.minimumCellSize()
312 shapeController_.minimumCellSize()
316 vertices[vI].alignment() = alignments[vI];
319 Info<<
" Clipped minimum size" <<
endl;
335 keep = decomposition().positionOnThisProcessor(pt);
338 if (keep && geometryToConformTo_.wellOutside(pt, SMALL))
345 keepVertex.unset(vI);
351 const label preInsertedSize = mesh_.number_of_vertices();
357 bool insertPoint =
false;
363 mesh_.dimension() < 3
373 const scalar interpolatedCellSize = shapeController_.cellSize(pt);
374 const triad interpolatedAlignment =
375 shapeController_.cellAlignment(pt);
376 const scalar calculatedCellSize =
vertices[vI].targetCellSize();
377 const triad calculatedAlignment =
vertices[vI].alignment();
381 Info<<
"Point = " << pt <<
nl 382 <<
" Size(interp) = " << interpolatedCellSize <<
nl 383 <<
" Size(calc) = " << calculatedCellSize <<
nl 384 <<
" Align(interp) = " << interpolatedAlignment <<
nl 385 <<
" Align(calc) = " << calculatedAlignment <<
nl 389 const scalar sizeDiff =
390 mag(interpolatedCellSize - calculatedCellSize);
391 const scalar alignmentDiff =
392 diff(interpolatedAlignment, calculatedAlignment);
396 Info<<
" size difference = " << sizeDiff <<
nl 397 <<
", alignment difference = " << alignmentDiff <<
endl;
403 sizeDiff/interpolatedCellSize > 0.1
404 || alignmentDiff > 0.15
410 if (forceInsertion || insertPoint)
412 const label oldSize = mesh_.vertexCount();
422 if (oldSize == mesh_.vertexCount() - 1)
426 insertedVert->index(),
427 controlFunction.maxPriority()
438 label(mesh_.number_of_vertices()) - preInsertedSize,
447 forAll(controlFunctions, fI)
449 const cellSizeAndAlignmentControl& controlFunction =
450 controlFunctions[fI];
452 const Switch forceInsertion =
453 controlFunction.forceInitialPointInsertion();
455 Info<<
"Inserting points from " << controlFunction.name()
456 <<
" (" << controlFunction.type() <<
")" <<
endl;
457 Info<<
" Force insertion is " << forceInsertion.c_str() <<
endl;
459 DynamicList<Foam::point> extraPts;
460 DynamicList<scalar> extraSizes;
462 controlFunction.cellSizeFunctionVertices(extraPts, extraSizes);
467 for (label vI = 0; vI < extraPts.size(); ++vI)
471 label maxPriority = -1;
472 scalar size = sizeControls_.cellSize(extraPts[vI], maxPriority);
474 if (maxPriority > controlFunction.maxPriority())
479 shapeController_.minimumCellSize()
482 else if (maxPriority == controlFunction.maxPriority())
486 min(extraSizes[vI], size),
487 shapeController_.minimumCellSize()
495 shapeController_.minimumCellSize()
510 keep = decomposition().positionOnThisProcessor(pt);
513 if (keep && geometryToConformTo_.wellOutside(pt, SMALL))
520 keepVertex.unset(vI);
526 const label preInsertedSize = mesh_.number_of_vertices();
530 bool insertPoint =
false;
536 mesh_.dimension() < 3
546 const scalar interpolatedCellSize = shapeController_.cellSize(pt);
547 const scalar calculatedCellSize =
vertices[vI].targetCellSize();
551 Info<<
"Point = " << pt <<
nl 552 <<
" Size(interp) = " << interpolatedCellSize <<
nl 553 <<
" Size(calc) = " << calculatedCellSize <<
nl 557 const scalar sizeDiff =
558 mag(interpolatedCellSize - calculatedCellSize);
562 Info<<
" size difference = " << sizeDiff <<
endl;
566 if (sizeDiff/interpolatedCellSize > 0.1)
571 if (forceInsertion || insertPoint)
629 Info<<
" Inserted extra points " 632 label(mesh_.number_of_vertices()) - preInsertedSize,
691 const autoPtr<backgroundMeshDecomposition>& decomposition
694 Info<<
"Iterate over " 695 <<
returnReduce(label(mesh_.number_of_finite_edges()), sumOp<label>())
696 <<
" cell size mesh edges" <<
endl;
698 DynamicList<Vb> verts(mesh_.number_of_vertices());
704 CellSizeDelaunay::Finite_edges_iterator eit =
705 mesh_.finite_edges_begin();
706 eit != mesh_.finite_edges_end();
710 if (
count % 10000 == 0)
712 Info<<
count <<
" edges, inserted " << verts.size()
713 <<
" Time = " << mesh_.time().elapsedCpuTime()
718 CellSizeDelaunay::Cell_handle
c = eit->first;
719 CellSizeDelaunay::Vertex_handle vA =
c->vertex(eit->second);
720 CellSizeDelaunay::Vertex_handle vB =
c->vertex(eit->third);
724 mesh_.is_infinite(vA)
725 || mesh_.is_infinite(vB)
726 || (vA->referred() && vB->referred())
727 || (vA->referred() && (vA->procIndex() > vB->procIndex()))
728 || (vB->referred() && (vB->procIndex() > vA->procIndex()))
739 const pointHit hitPt = findDiscontinuities(l);
745 if (!geometryToConformTo_.inside(pt))
752 if (!decomposition().positionOnThisProcessor(pt))
767 verts.last().targetCellSize() = sizeControls_.cellSize(pt);
772 mesh_.insertPoints(verts,
false);
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
void size(const label n)
Older name for setAddressableSize.
pointFromPoint topoint(const Point &P)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
static bool & parRun() noexcept
Test if this a parallel run.
void inplaceSubset(const BoolListType &select, ListType &input, const bool invert=false)
Inplace extract elements of the input list when select is true.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
#define forAll(list, i)
Loop across all elements in list.
void initialMeshPopulation(const autoPtr< backgroundMeshDecomposition > &decomposition)
pointField vertices(const blockVertexList &bvl)
CGAL::indexedVertex< K > Vb
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
label refineMesh(const autoPtr< backgroundMeshDecomposition > &decomposition)
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
line< point, const point & > linePointRef
A line using referred points.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
void clear()
Clear the list, i.e. set size to zero.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
PointFrompoint toPoint(const Foam::point &p)
Describes the interaction of a object and a (templated) point. It carries the info of a successful hi...
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
~controlMeshRefinement()
Destructor.
Field< triad > triadField
Specialisation of Field<T> for triad.
CellSizeDelaunay::Vertex_handle Vertex_handle
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
vector point
Point is a vector.
static constexpr direction size() noexcept
The number of elements in the VectorSpace = Ncmpts.
const dimensionedScalar c
Speed of light in a vacuum.
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
PointHit< point > pointHit
A PointHit with a 3D point.