44Foam::scalar Foam::controlMeshRefinement::calcFirstDerivative
47 const scalar& cellSizeA,
49 const scalar& cellSizeB
52 return (cellSizeA - cellSizeB)/
mag(a -
b);
70bool Foam::controlMeshRefinement::detectEdge
76 const scalar secondDerivTolSqr
84 label nIterations = 0;
95 pointFound.setPoint(midPoint);
103 scalar cellSizeA = sizeControls_.cellSize(a);
104 scalar cellSizeB = sizeControls_.cellSize(
b);
111 scalar cellSizeMid = sizeControls_.cellSize(midPoint);
115 const scalar cellSizeMid1 = sizeControls_.cellSize(midPoint1);
120 scalar secondDerivative1 =
133 const scalar cellSizeMid2 = sizeControls_.cellSize(midPoint2);
138 scalar secondDerivative2 =
153 magSqr(secondDerivative1) < secondDerivTolSqr
154 &&
magSqr(secondDerivative2) < secondDerivTolSqr
161 if (
magSqr(secondDerivative1) >
magSqr(secondDerivative2))
164 midPoint = midPoint1;
169 midPoint = midPoint2;
182 const scalar tolSqr =
sqr(1
e-3);
183 const scalar secondDerivTolSqr =
sqr(1
e-3);
205 cellShapeControl& shapeController
208 shapeController_(shapeController),
209 mesh_(shapeController.shapeControlMesh()),
210 sizeControls_(shapeController.sizeAndAlignment()),
211 geometryToConformTo_(sizeControls_.geometryToConformTo())
225 const autoPtr<backgroundMeshDecomposition>& decomposition
228 if (shapeController_.shapeControlMesh().vertexCount() > 0)
231 Info<<
"Cell size and alignment mesh already populated." <<
endl;
235 autoPtr<boundBox> overallBoundBox;
257 Map<label> priorityMap;
259 const PtrList<cellSizeAndAlignmentControl>& controlFunctions =
260 sizeControls_.controlFunctions();
262 forAll(controlFunctions, fI)
264 const cellSizeAndAlignmentControl& controlFunction =
265 controlFunctions[fI];
267 const Switch& forceInsertion =
268 controlFunction.forceInitialPointInsertion();
270 Info<<
"Inserting points from " << controlFunction.name()
271 <<
" (" << controlFunction.type() <<
")" <<
endl;
272 Info<<
" Force insertion is " << forceInsertion.c_str() <<
endl;
278 controlFunction.initialVertices(pts, sizes, alignments);
280 Info<<
" Got initial vertices list of size " << pts.size() <<
endl;
285 for (label vI = 0; vI < pts.size(); ++vI)
289 label maxPriority = -1;
290 scalar size = sizeControls_.cellSize(pts[vI], maxPriority);
292 if (maxPriority > controlFunction.maxPriority())
297 shapeController_.minimumCellSize()
313 shapeController_.minimumCellSize()
317 vertices[vI].alignment() = alignments[vI];
320 Info<<
" Clipped minimum size" <<
endl;
336 keep = decomposition().positionOnThisProcessor(pt);
339 if (keep && geometryToConformTo_.wellOutside(pt, SMALL))
346 keepVertex.unset(vI);
352 const label preInsertedSize = mesh_.number_of_vertices();
358 bool insertPoint =
false;
364 mesh_.dimension() < 3
374 const scalar interpolatedCellSize = shapeController_.cellSize(pt);
375 const triad interpolatedAlignment =
376 shapeController_.cellAlignment(pt);
377 const scalar calculatedCellSize =
vertices[vI].targetCellSize();
378 const triad calculatedAlignment =
vertices[vI].alignment();
382 Info<<
"Point = " << pt <<
nl
383 <<
" Size(interp) = " << interpolatedCellSize <<
nl
384 <<
" Size(calc) = " << calculatedCellSize <<
nl
385 <<
" Align(interp) = " << interpolatedAlignment <<
nl
386 <<
" Align(calc) = " << calculatedAlignment <<
nl
390 const scalar sizeDiff =
391 mag(interpolatedCellSize - calculatedCellSize);
392 const scalar alignmentDiff =
393 diff(interpolatedAlignment, calculatedAlignment);
397 Info<<
" size difference = " << sizeDiff <<
nl
398 <<
", alignment difference = " << alignmentDiff <<
endl;
404 sizeDiff/interpolatedCellSize > 0.1
405 || alignmentDiff > 0.15
411 if (forceInsertion || insertPoint)
413 const label oldSize = mesh_.vertexCount();
423 if (oldSize == mesh_.vertexCount() - 1)
427 insertedVert->index(),
428 controlFunction.maxPriority()
439 label(mesh_.number_of_vertices()) - preInsertedSize,
448 forAll(controlFunctions, fI)
450 const cellSizeAndAlignmentControl& controlFunction =
451 controlFunctions[fI];
453 const Switch& forceInsertion =
454 controlFunction.forceInitialPointInsertion();
456 Info<<
"Inserting points from " << controlFunction.name()
457 <<
" (" << controlFunction.type() <<
")" <<
endl;
458 Info<<
" Force insertion is " << forceInsertion.c_str() <<
endl;
460 DynamicList<Foam::point> extraPts;
461 DynamicList<scalar> extraSizes;
463 controlFunction.cellSizeFunctionVertices(extraPts, extraSizes);
468 for (label vI = 0; vI < extraPts.size(); ++vI)
472 label maxPriority = -1;
473 scalar size = sizeControls_.cellSize(extraPts[vI], maxPriority);
475 if (maxPriority > controlFunction.maxPriority())
480 shapeController_.minimumCellSize()
483 else if (maxPriority == controlFunction.maxPriority())
487 min(extraSizes[vI], size),
488 shapeController_.minimumCellSize()
496 shapeController_.minimumCellSize()
511 keep = decomposition().positionOnThisProcessor(pt);
514 if (keep && geometryToConformTo_.wellOutside(pt, SMALL))
521 keepVertex.unset(vI);
527 const label preInsertedSize = mesh_.number_of_vertices();
531 bool insertPoint =
false;
537 mesh_.dimension() < 3
547 const scalar interpolatedCellSize = shapeController_.cellSize(pt);
548 const scalar calculatedCellSize =
vertices[vI].targetCellSize();
552 Info<<
"Point = " << pt <<
nl
553 <<
" Size(interp) = " << interpolatedCellSize <<
nl
554 <<
" Size(calc) = " << calculatedCellSize <<
nl
558 const scalar sizeDiff =
559 mag(interpolatedCellSize - calculatedCellSize);
563 Info<<
" size difference = " << sizeDiff <<
endl;
567 if (sizeDiff/interpolatedCellSize > 0.1)
572 if (forceInsertion || insertPoint)
630 Info<<
" Inserted extra points "
633 label(mesh_.number_of_vertices()) - preInsertedSize,
692 const autoPtr<backgroundMeshDecomposition>& decomposition
695 Info<<
"Iterate over "
696 <<
returnReduce(label(mesh_.number_of_finite_edges()), sumOp<label>())
697 <<
" cell size mesh edges" <<
endl;
699 DynamicList<Vb> verts(mesh_.number_of_vertices());
705 CellSizeDelaunay::Finite_edges_iterator eit =
706 mesh_.finite_edges_begin();
707 eit != mesh_.finite_edges_end();
711 if (count % 10000 == 0)
713 Info<<
count <<
" edges, inserted " << verts.size()
714 <<
" Time = " << mesh_.time().elapsedCpuTime()
719 CellSizeDelaunay::Cell_handle
c = eit->first;
720 CellSizeDelaunay::Vertex_handle vA =
c->vertex(eit->second);
721 CellSizeDelaunay::Vertex_handle vB =
c->vertex(eit->third);
725 mesh_.is_infinite(vA)
726 || mesh_.is_infinite(vB)
727 || (vA->referred() && vB->referred())
728 || (vA->referred() && (vA->procIndex() > vB->procIndex()))
729 || (vB->referred() && (vB->procIndex() > vA->procIndex()))
740 const pointHit hitPt = findDiscontinuities(l);
746 if (!geometryToConformTo_.inside(pt))
753 if (!decomposition().positionOnThisProcessor(pt))
768 verts.last().targetCellSize() = sizeControls_.cellSize(pt);
773 mesh_.insertPoints(verts,
false);
CGAL::indexedVertex< K > Vb
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Describes the interaction of a face and a point. It carries the info of a successful hit and (if succ...
void size(const label n)
Older name for setAddressableSize.
static bool & parRun() noexcept
Test if this a parallel run.
static constexpr direction size() noexcept
The number of elements in the VectorSpace = Ncmpts.
CellSizeDelaunay::Vertex_handle Vertex_handle
void initialMeshPopulation(const autoPtr< backgroundMeshDecomposition > &decomposition)
~controlMeshRefinement()
Destructor.
label refineMesh(const autoPtr< backgroundMeshDecomposition > &decomposition)
static const complex max
complex (VGREAT,VGREAT)
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
const dimensionedScalar c
Speed of light in a vacuum.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
void inplaceSubset(const BoolListType &select, ListType &input, const bool invert=false)
Inplace extract elements of the input list when select is true.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointFrompoint toPoint(const Foam::point &p)
line< point, const point & > linePointRef
A line using referred points.
pointFromPoint topoint(const Point &P)
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
Field< triad > triadField
Specialisation of Field<T> for triad.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
pointField vertices(const blockVertexList &bvl)
PointHit< point > pointHit
A PointIndexHit for 3D points.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.