34 void Foam::refinementFeatures::read
36 const objectRegistry& io,
37 const PtrList<dictionary>& featDicts
42 const dictionary& dict = featDicts[featI];
44 fileName featFileName(dict.lookup(
"file"));
53 "extendedFeatureEdgeMesh",
60 const fileName fName(typeFilePath<extendedFeatureEdgeMesh>(extFeatObj));
69 Info<<
"Read extendedFeatureEdgeMesh " << extFeatObj.name()
71 eMeshPtr().writeStats(
Info);
74 set(featI,
new extendedFeatureEdgeMesh(extFeatObj, eMeshPtr()));
91 const fileName fName(typeFilePath<featureEdgeMesh>(featObj));
98 ) <<
"Could not open " << featObj.objectPath()
105 const edgeMesh& eMesh = eMeshPtr();
107 Info<<
"Read edgeMesh " << featObj.name() <<
nl 109 eMesh.writeStats(
Info);
117 labelList oldToNew(eMesh.points().size(), -1);
118 DynamicField<point> newPoints(eMesh.points().size());
119 forAll(pointEdges, pointi)
121 if (pointEdges[pointi].
size() > 2)
123 oldToNew[pointi] = newPoints.size();
124 newPoints.append(eMesh.points()[pointi]);
129 label nFeatures = newPoints.size();
132 if (oldToNew[pointi] == -1)
134 oldToNew[pointi] = newPoints.size();
135 newPoints.append(eMesh.points()[pointi]);
140 const edgeList& edges = eMesh.edges();
144 const edge&
e = edges[edgeI];
145 newEdges[edgeI] = edge
156 extendedEdgeMesh eeMesh
168 List<extendedEdgeMesh::sideVolumeType>(0),
182 set(featI,
new extendedFeatureEdgeMesh(featObj, eeMesh));
185 const extendedEdgeMesh& eMesh =
operator[](featI);
187 if (dict.found(
"levels"))
189 List<Tuple2<scalar, label>> distLevels(dict[
"levels"]);
194 <<
" : levels should be at least size 1" << endl
195 <<
"levels : " << dict[
"levels"]
199 distances_[featI].
setSize(distLevels.size());
200 levels_[featI].
setSize(distLevels.size());
204 distances_[featI][j] = distLevels[j].
first();
205 levels_[featI][j] = distLevels[j].second();
212 (distances_[featI][j] <= distances_[featI][j-1])
213 || (levels_[featI][j] > levels_[featI][j-1])
217 <<
" : Refinement should be specified in order" 218 <<
" of increasing distance" 219 <<
" (and decreasing refinement level)." << endl
220 <<
"Distance:" << distances_[featI][j]
221 <<
" refinementLevel:" << levels_[featI][j]
234 Info<<
"Refinement level according to distance to " 235 << featFileName <<
" (" << eMesh.points().size() <<
" points, " 236 << eMesh.edges().size() <<
" edges)." <<
endl;
239 Info<<
" level " << levels_[featI][j]
240 <<
" for all cells within " << distances_[featI][j]
241 <<
" metre." <<
endl;
247 void Foam::refinementFeatures::buildTrees(
const label featI)
249 const extendedEdgeMesh& eMesh =
operator[](featI);
251 const edgeList& edges = eMesh.edges();
254 treeBoundBox bb(points);
262 bb.min() -=
point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
263 bb.max() +=
point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
268 new indexedOctree<treeDataEdge>
290 new indexedOctree<treeDataPoint>
292 treeDataPoint(points, featurePoints),
303 void Foam::refinementFeatures::findHigherLevel
310 const labelList& levels = levels_[featI];
321 label candidateI = 0;
327 if (levels[levelI] > maxLevel[pointi])
329 candidates[candidateI] = pt[pointi];
330 candidateMap[candidateI] = pointi;
331 candidateDistSqr[candidateI] =
sqr(distances[levelI]);
337 candidates.setSize(candidateI);
338 candidateMap.setSize(candidateI);
339 candidateDistSqr.setSize(candidateI);
342 const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
344 List<pointIndexHit>
nearInfo(candidates.size());
345 forAll(candidates, candidateI)
347 nearInfo[candidateI] = tree.findNearest
349 candidates[candidateI],
350 candidateDistSqr[candidateI]
355 forAll(nearInfo, candidateI)
357 if (nearInfo[candidateI].hit())
363 mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
366 label pointi = candidateMap[candidateI];
369 maxLevel[pointi] = levels[minDistI+1];
380 if (!regionEdgeTreesPtr_.valid())
382 regionEdgeTreesPtr_.reset
403 bb.
min() -=
point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
404 bb.
max() +=
point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
426 return regionEdgeTreesPtr_();
439 distances_(featDicts.
size()),
440 levels_(featDicts.
size()),
441 edgeTrees_(featDicts.
size()),
442 pointTrees_(featDicts.
size())
547 if (tree.
shapes().size() > 0)
551 const point& sample = samples[sampleI];
554 if (nearInfo[sampleI].hit())
556 distSqr =
magSqr(nearInfo[sampleI].hitPoint()-sample);
560 distSqr = nearestDistSqr[sampleI];
567 nearFeature[sampleI] = featI;
576 const edge&
e = td.edges()[nearInfo[sampleI].index()];
577 nearNormal[sampleI] = e.
vec(td.points());
578 nearNormal[sampleI] /=
mag(nearNormal[sampleI])+VSMALL;
606 forAll(regionTrees, featI)
612 const point& sample = samples[sampleI];
615 if (nearInfo[sampleI].hit())
617 distSqr =
magSqr(nearInfo[sampleI].hitPoint()-sample);
621 distSqr = nearestDistSqr[sampleI];
625 pointIndexHit info = regionTree.findNearest(sample, distSqr);
631 nearFeature[sampleI] = featI;
639 const edge&
e = td.
edges()[nearInfo[sampleI].index()];
640 nearNormal[sampleI] = e.vec(td.
points());
641 nearNormal[sampleI] /=
mag(nearNormal[sampleI])+VSMALL;
714 forAll(pointTrees_, featI)
718 if (tree.
shapes().pointLabels().size() > 0)
722 const point& sample = samples[sampleI];
725 if (nearFeature[sampleI] != -1)
727 distSqr =
magSqr(nearInfo[sampleI].hitPoint()-sample);
731 distSqr = nearestDistSqr[sampleI];
738 nearFeature[sampleI] = featI;
752 void Foam::refinementFeatures::findHigherLevel
764 findHigherLevel(pt, featI, maxLevel);
771 scalar overallMax = -GREAT;
774 overallMax =
max(overallMax,
max(distances_[featI]));
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.
const pointField & points() const
bool set(const label) const
Is element set.
errorManipArg< error, int > exit(error &err, const int errNo=1)
refinementFeatures(const objectRegistry &io, const PtrList< dictionary > &featDicts)
Construct from description.
scalar maxDistance() const
Highest distance of all features.
const double e
Elementary charge.
void findNearestPoint(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo) const
Find nearest feature point. Sets:
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)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Description of feature edges and points.
const T & operator[](const label) const
Return element const reference.
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.
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
const labelList & regionEdges() const
Return the feature edges which are on the boundary between.
T & first()
Return the first element of the list.
const Point & hitPoint() const
Return hit point.
vectorField pointField
pointField is a vectorField.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const edgeList & edges() const
Return edges.
cachedRandom rndGen(label(0), -1)
bool hit() const
Is there a hit.
const Type & shapes() const
Reference to shape.
List< label > labelList
A List of labels.
label readLabel(Istream &is)
Simple random number generator.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
vector vec(const pointField &) const
Return the vector (end - start)
const pointField & points() const
Return points.
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
treeBoundBox extend(Random &, const scalar s) const
Return slightly wider bounding box.
label size() const
Return the number of elements in the UPtrList.
void setSize(const label)
Reset size of List.
const point & max() const
Maximum describing the bounding box.
void findNearestRegionEdge(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo, vectorField &nearNormal) const
Find nearest point on nearest region edge. Sets:
const edgeList & edges() const
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...
static bool canRead(const fileName &, const bool verbose=false)
Can we read this file format?
const PtrList< indexedOctree< treeDataEdge > > & regionEdgeTrees() const
void findNearestEdge(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo, vectorField &nearNormal) const
Find nearest point on nearest feature edge. Sets:
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 > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
const point & min() const
Minimum describing the bounding box.
Registry of regIOobjects.
Ostream & incrIndent(Ostream &os)
Increment the indent level.
label index() const
Return index.
Tuple2< scalar, label > nearInfo
Private class for finding nearest.