33 void Foam::refinementFeatures::read
35 const objectRegistry& io,
36 const PtrList<dictionary>& featDicts
41 const dictionary& dict = featDicts[featI];
43 fileName featFileName(dict.lookup(
"file"));
52 "extendedFeatureEdgeMesh",
59 const fileName fName(extFeatObj.filePath());
68 Info<<
"Read extendedFeatureEdgeMesh " << extFeatObj.name()
70 eMeshPtr().writeStats(
Info);
73 set(featI,
new extendedFeatureEdgeMesh(extFeatObj, eMeshPtr()));
90 const fileName fName(featObj.filePath());
97 ) <<
"Could not open " << featObj.objectPath()
104 const edgeMesh& eMesh = eMeshPtr();
106 Info<<
"Read edgeMesh " << featObj.name() <<
nl 108 eMesh.writeStats(
Info);
116 labelList oldToNew(eMesh.points().size(), -1);
117 DynamicField<point> newPoints(eMesh.points().size());
118 forAll(pointEdges, pointi)
120 if (pointEdges[pointi].
size() > 2)
122 oldToNew[pointi] = newPoints.size();
123 newPoints.append(eMesh.points()[pointi]);
128 label nFeatures = newPoints.size();
131 if (oldToNew[pointi] == -1)
133 oldToNew[pointi] = newPoints.size();
134 newPoints.append(eMesh.points()[pointi]);
139 const edgeList& edges = eMesh.edges();
143 const edge&
e = edges[edgeI];
144 newEdges[edgeI] = edge
155 extendedEdgeMesh eeMesh
167 List<extendedEdgeMesh::sideVolumeType>(0),
181 set(featI,
new extendedFeatureEdgeMesh(featObj, eeMesh));
184 const extendedEdgeMesh& eMesh =
operator[](featI);
186 if (dict.found(
"levels"))
188 List<Tuple2<scalar, label>> distLevels(dict[
"levels"]);
193 <<
" : levels should be at least size 1" << endl
194 <<
"levels : " << dict[
"levels"]
198 distances_[featI].
setSize(distLevels.size());
199 levels_[featI].
setSize(distLevels.size());
203 distances_[featI][j] = distLevels[j].
first();
204 levels_[featI][j] = distLevels[j].second();
211 (distances_[featI][j] <= distances_[featI][j-1])
212 || (levels_[featI][j] > levels_[featI][j-1])
216 <<
" : Refinement should be specified in order" 217 <<
" of increasing distance" 218 <<
" (and decreasing refinement level)." << endl
219 <<
"Distance:" << distances_[featI][j]
220 <<
" refinementLevel:" << levels_[featI][j]
233 Info<<
"Refinement level according to distance to " 234 << featFileName <<
" (" << eMesh.points().size() <<
" points, " 235 << eMesh.edges().size() <<
" edges)." <<
endl;
238 Info<<
" level " << levels_[featI][j]
239 <<
" for all cells within " << distances_[featI][j]
240 <<
" metre." <<
endl;
246 void Foam::refinementFeatures::buildTrees(
const label featI)
248 const extendedEdgeMesh& eMesh =
operator[](featI);
250 const edgeList& edges = eMesh.edges();
253 treeBoundBox bb(points);
261 bb.min() -=
point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
262 bb.max() +=
point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
267 new indexedOctree<treeDataEdge>
289 new indexedOctree<treeDataPoint>
291 treeDataPoint(points, featurePoints),
304 if (!regionEdgeTreesPtr_.valid())
306 regionEdgeTreesPtr_.reset
327 bb.
min() -=
point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
328 bb.
max() +=
point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
350 return regionEdgeTreesPtr_();
355 void Foam::refinementFeatures::findHigherLevel
373 label candidateI = 0;
379 if (levels[levelI] > maxLevel[pointi])
381 candidates[candidateI] = pt[pointi];
382 candidateMap[candidateI] = pointi;
383 candidateDistSqr[candidateI] =
sqr(distances[levelI]);
389 candidates.setSize(candidateI);
390 candidateMap.setSize(candidateI);
391 candidateDistSqr.setSize(candidateI);
397 forAll(candidates, candidateI)
399 nearInfo[candidateI] = tree.findNearest
401 candidates[candidateI],
402 candidateDistSqr[candidateI]
407 forAll(nearInfo, candidateI)
409 if (nearInfo[candidateI].hit())
415 mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
418 label pointi = candidateMap[candidateI];
421 maxLevel[pointi] = levels[minDistI+1];
436 distances_(featDicts.
size()),
437 levels_(featDicts.
size()),
438 edgeTrees_(featDicts.
size()),
439 pointTrees_(featDicts.
size())
544 if (tree.
shapes().size() > 0)
548 const point& sample = samples[sampleI];
551 if (nearInfo[sampleI].hit())
553 distSqr =
magSqr(nearInfo[sampleI].hitPoint()-sample);
557 distSqr = nearestDistSqr[sampleI];
564 nearFeature[sampleI] = featI;
573 const edge&
e = td.edges()[nearInfo[sampleI].index()];
574 nearNormal[sampleI] = e.
vec(td.points());
575 nearNormal[sampleI] /=
mag(nearNormal[sampleI])+VSMALL;
603 forAll(regionTrees, featI)
609 const point& sample = samples[sampleI];
612 if (nearInfo[sampleI].hit())
614 distSqr =
magSqr(nearInfo[sampleI].hitPoint()-sample);
618 distSqr = nearestDistSqr[sampleI];
622 pointIndexHit info = regionTree.findNearest(sample, distSqr);
628 nearFeature[sampleI] = featI;
636 const edge&
e = td.
edges()[nearInfo[sampleI].index()];
637 nearNormal[sampleI] = e.vec(td.
points());
638 nearNormal[sampleI] /=
mag(nearNormal[sampleI])+VSMALL;
711 forAll(pointTrees_, featI)
715 if (tree.
shapes().pointLabels().size() > 0)
719 const point& sample = samples[sampleI];
722 if (nearFeature[sampleI] != -1)
724 distSqr =
magSqr(nearInfo[sampleI].hitPoint()-sample);
728 distSqr = nearestDistSqr[sampleI];
735 nearFeature[sampleI] = featI;
749 void Foam::refinementFeatures::findHigherLevel
761 findHigherLevel(pt, featI, maxLevel);
768 scalar overallMax = -GREAT;
771 overallMax =
max(overallMax,
max(distances_[featI]));
cachedRandom rndGen(label(0),-1)
List< labelList > labelListList
A List of labelList.
static autoPtr< extendedEdgeMesh > New(const fileName &, const word &ext)
Select constructed from filename (explicit extension)
#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.
errorManipArg< error, int > exit(error &err, const int errNo=1)
refinementFeatures(const objectRegistry &io, const PtrList< dictionary > &featDicts)
Construct from description.
const double e
Elementary charge.
const point & min() const
Minimum describing the bounding box.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
static autoPtr< edgeMesh > New(const fileName &, const word &ext)
Select constructed from filename (explicit extension)
bool hit() const
Is there a hit.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Description of feature edges and points.
treeBoundBox extend(Random &, const scalar s) const
Return slightly wider bounding box.
Ostream & endl(Ostream &os)
Add newline and flush stream.
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Holds data for octree to work on an edges subset.
const point & max() const
Maximum describing the bounding box.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
PointIndexHit< point > pointIndexHit
T & first()
Return the first element of the list.
const edgeList & edges() const
const pointField & points() const
Return points.
vectorField pointField
pointField is a vectorField.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
scalar maxDistance() const
Highest distance of all features.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void findNearestPoint(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo) const
Find nearest feature point. Sets:
bool set(const label) const
Is element set.
List< label > labelList
A List of labels.
label readLabel(Istream &is)
Simple random number generator.
const labelList & regionEdges() const
Return the feature edges which are on the boundary between.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const List< scalarField > & distances() const
Per featureEdgeMesh the list of ranges.
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
const T & operator[](const label) const
Return element const reference.
void setSize(const label)
Reset size of List.
vector point
Point is a vector.
Non-pointer based hierarchical recursive searching.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
const Type & shapes() const
Reference to shape.
label index() const
Return index.
static bool canRead(const fileName &, const bool verbose=false)
Can we read this file format?
const edgeList & edges() const
Return edges.
const Point & hitPoint() const
Return hit point.
vector vec(const pointField &) const
Return the vector (end - start)
const pointField & points() const
Standard boundBox + extra functionality for use in octree.
label findLower(const ListType &, typename ListType::const_reference, const label stary, const BinaryOp &bop)
Find last element < given value in sorted list and return index,.
dimensioned< scalar > mag(const dimensioned< Type > &)
const PtrList< indexedOctree< treeDataEdge > > & regionEdgeTrees() const
Field< vector > vectorField
Specialisation of Field<T> for vector.
Registry of regIOobjects.
const labelListList & levels() const
Per featureEdgeMesh the list of level.
Ostream & incrIndent(Ostream &os)
Increment the indent level.
void findNearestEdge(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo, vectorField &nearNormal) const
Find nearest point on nearest feature edge. Sets:
void findNearestRegionEdge(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo, vectorField &nearNormal) const
Find nearest point on nearest region edge. Sets:
Tuple2< scalar, label > nearInfo
Private class for finding nearest.
label size() const
Return the number of elements in the UPtrList.