42 Foam::scalar Foam::controlMeshRefinement::calcFirstDerivative
45 const scalar& cellSizeA,
47 const scalar& cellSizeB
50 return (cellSizeA - cellSizeB)/
mag(a - b);
68 bool Foam::controlMeshRefinement::detectEdge
74 const scalar secondDerivTolSqr
82 label nIterations = 0;
93 pointFound.setPoint(midPoint);
101 scalar cellSizeA = sizeControls_.cellSize(a);
102 scalar cellSizeB = sizeControls_.cellSize(b);
109 scalar cellSizeMid = sizeControls_.cellSize(midPoint);
113 const scalar cellSizeMid1 = sizeControls_.cellSize(midPoint1);
118 scalar secondDerivative1 =
131 const scalar cellSizeMid2 = sizeControls_.cellSize(midPoint2);
136 scalar secondDerivative2 =
151 magSqr(secondDerivative1) < secondDerivTolSqr
152 &&
magSqr(secondDerivative2) < secondDerivTolSqr
159 if (
magSqr(secondDerivative1) >
magSqr(secondDerivative2))
162 midPoint = midPoint1;
167 midPoint = midPoint2;
180 const scalar tolSqr =
sqr(1
e-3);
181 const scalar secondDerivTolSqr =
sqr(1
e-3);
201 Foam::controlMeshRefinement::controlMeshRefinement
203 cellShapeControl& shapeController
206 shapeController_(shapeController),
207 mesh_(shapeController.shapeControlMesh()),
208 sizeControls_(shapeController.sizeAndAlignment()),
209 geometryToConformTo_(sizeControls_.geometryToConformTo())
223 const autoPtr<backgroundMeshDecomposition>& decomposition
226 if (shapeController_.shapeControlMesh().vertexCount() > 0)
229 Info<<
"Cell size and alignment mesh already populated." <<
endl;
233 autoPtr<boundBox> overallBoundBox;
255 Map<label> priorityMap;
257 const PtrList<cellSizeAndAlignmentControl>& controlFunctions =
258 sizeControls_.controlFunctions();
260 forAll(controlFunctions, fI)
263 controlFunctions[fI];
265 const Switch& forceInsertion =
266 controlFunction.forceInitialPointInsertion();
268 Info<<
"Inserting points from " << controlFunction.name()
269 <<
" (" << controlFunction.type() <<
")" <<
endl;
270 Info<<
" Force insertion is " << forceInsertion.asText() <<
endl;
276 controlFunction.initialVertices(pts, sizes, alignments);
278 Info<<
" Got initial vertices list of size " << pts.size() <<
endl;
280 List<Vb> vertices(pts.size());
283 for (
label vI = 0; vI < pts.size(); ++vI)
287 label maxPriority = -1;
288 scalar size = sizeControls_.cellSize(pts[vI], maxPriority);
290 if (maxPriority > controlFunction.maxPriority())
292 vertices[vI].targetCellSize() =
max 295 shapeController_.minimumCellSize()
308 vertices[vI].targetCellSize() =
max 311 shapeController_.minimumCellSize()
315 vertices[vI].alignment() = alignments[vI];
318 Info<<
" Clipped minimum size" <<
endl;
324 PackedBoolList keepVertex(vertices.size(),
true);
334 keep = decomposition().positionOnThisProcessor(pt);
337 if (keep && geometryToConformTo_.
wellOutside(pt, SMALL))
344 keepVertex[vI] =
false;
350 const label preInsertedSize = mesh_.number_of_vertices();
356 bool insertPoint =
false;
362 mesh_.dimension() < 3
365 mesh_.locate(vertices[vI].point())
372 const scalar interpolatedCellSize = shapeController_.cellSize(pt);
373 const triad interpolatedAlignment =
374 shapeController_.cellAlignment(pt);
375 const scalar calculatedCellSize = vertices[vI].targetCellSize();
376 const triad calculatedAlignment = vertices[vI].alignment();
380 Info<<
"Point = " << pt <<
nl 381 <<
" Size(interp) = " << interpolatedCellSize <<
nl 382 <<
" Size(calc) = " << calculatedCellSize <<
nl 383 <<
" Align(interp) = " << interpolatedAlignment <<
nl 384 <<
" Align(calc) = " << calculatedAlignment <<
nl 388 const scalar sizeDiff =
389 mag(interpolatedCellSize - calculatedCellSize);
390 const scalar alignmentDiff =
391 diff(interpolatedAlignment, calculatedAlignment);
395 Info<<
" size difference = " << sizeDiff <<
nl 396 <<
", alignment difference = " << alignmentDiff <<
endl;
402 sizeDiff/interpolatedCellSize > 0.1
403 || alignmentDiff > 0.15
409 if (forceInsertion || insertPoint)
411 const label oldSize = mesh_.vertexCount();
417 vertices[vI].alignment(),
421 if (oldSize == mesh_.vertexCount() - 1)
425 insertedVert->index(),
426 controlFunction.maxPriority()
437 label(mesh_.number_of_vertices()) - preInsertedSize,
446 forAll(controlFunctions, fI)
449 controlFunctions[fI];
451 const Switch& forceInsertion =
452 controlFunction.forceInitialPointInsertion();
454 Info<<
"Inserting points from " << controlFunction.name()
455 <<
" (" << controlFunction.type() <<
")" <<
endl;
456 Info<<
" Force insertion is " << forceInsertion.asText() <<
endl;
458 DynamicList<Foam::point> extraPts;
459 DynamicList<scalar> extraSizes;
461 controlFunction.cellSizeFunctionVertices(extraPts, extraSizes);
463 List<Vb> vertices(extraPts.size());
466 for (
label vI = 0; vI < extraPts.size(); ++vI)
470 label maxPriority = -1;
471 scalar size = sizeControls_.cellSize(extraPts[vI], maxPriority);
473 if (maxPriority > controlFunction.maxPriority())
475 vertices[vI].targetCellSize() =
max 478 shapeController_.minimumCellSize()
481 else if (maxPriority == controlFunction.maxPriority())
483 vertices[vI].targetCellSize() =
max 485 min(extraSizes[vI], size),
486 shapeController_.minimumCellSize()
491 vertices[vI].targetCellSize() =
max 494 shapeController_.minimumCellSize()
499 PackedBoolList keepVertex(vertices.size(),
true);
509 keep = decomposition().positionOnThisProcessor(pt);
512 if (keep && geometryToConformTo_.
wellOutside(pt, SMALL))
519 keepVertex[vI] =
false;
525 const label preInsertedSize = mesh_.number_of_vertices();
529 bool insertPoint =
false;
535 mesh_.dimension() < 3
538 mesh_.locate(vertices[vI].point())
545 const scalar interpolatedCellSize = shapeController_.cellSize(pt);
546 const scalar calculatedCellSize = vertices[vI].targetCellSize();
550 Info<<
"Point = " << pt <<
nl 551 <<
" Size(interp) = " << interpolatedCellSize <<
nl 552 <<
" Size(calc) = " << calculatedCellSize <<
nl 556 const scalar sizeDiff =
557 mag(interpolatedCellSize - calculatedCellSize);
561 Info<<
" size difference = " << sizeDiff <<
endl;
565 if (sizeDiff/interpolatedCellSize > 0.1)
570 if (forceInsertion || insertPoint)
610 vertices[vI].alignment(),
628 Info<<
" Inserted extra points " 631 label(mesh_.number_of_vertices()) - preInsertedSize,
690 const autoPtr<backgroundMeshDecomposition>& decomposition
693 Info<<
"Iterate over " 695 <<
" cell size mesh edges" << endl;
697 DynamicList<Vb> verts(mesh_.number_of_vertices());
703 CellSizeDelaunay::Finite_edges_iterator eit =
704 mesh_.finite_edges_begin();
705 eit != mesh_.finite_edges_end();
709 if (count % 10000 == 0)
711 Info<< count <<
" edges, inserted " << verts.size()
712 <<
" Time = " << mesh_.time().elapsedCpuTime()
717 CellSizeDelaunay::Cell_handle c = eit->first;
718 CellSizeDelaunay::Vertex_handle vA = c->vertex(eit->second);
719 CellSizeDelaunay::Vertex_handle vB = c->vertex(eit->third);
723 mesh_.is_infinite(vA)
724 || mesh_.is_infinite(vB)
725 || (vA->referred() && vB->referred())
726 || (vA->referred() && (vA->procIndex() > vB->procIndex()))
727 || (vB->referred() && (vB->procIndex() > vA->procIndex()))
738 const pointHit hitPt = findDiscontinuities(l);
744 if (!geometryToConformTo_.
inside(pt))
751 if (!decomposition().positionOnThisProcessor(pt))
766 verts.last().targetCellSize() = sizeControls_.cellSize(pt);
771 mesh_.insertPoints(verts,
false);
void inplaceSubset(const UList< T > &select, const T &value, ListType &)
Inplace extract elements of List when select is a certain value.
static const Vector< scalar > max
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
pointFromPoint topoint(const Point &P)
const double e
Elementary charge.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
cellSizeAndAlignmentControl(const Time &runTime, const word &name, const dictionary &controlFunctionDict, const conformationSurfaces &geometryToConformTo, const scalar &defaultCellSize)
Construct from dictionary and references to conformalVoronoiMesh and.
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
void initialMeshPopulation(const autoPtr< backgroundMeshDecomposition > &decomposition)
CGAL::indexedVertex< K > Vb
label refineMesh(const autoPtr< backgroundMeshDecomposition > &decomposition)
vectorField pointField
pointField is a vectorField.
line< point, const point & > linePointRef
Line using referred points.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
PointFrompoint toPoint(const Foam::point &p)
~controlMeshRefinement()
Destructor.
Field< triad > triadField
Specialisation of Field<T> for triad.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
CellSizeDelaunay::Vertex_handle Vertex_handle
static direction size()
Return the number of elements in the VectorSpace = Ncmpts.
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
static bool & parRun()
Is this a parallel run?
vector point
Point is a vector.
dimensioned< scalar > mag(const dimensioned< Type > &)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
PointHit< point > pointHit