33 #include "indexedVertex.H" 57 template<
class Triangulation,
class Type>
60 const Triangulation& mesh,
72 typename Triangulation::Finite_vertices_iterator vit =
73 mesh.finite_vertices_begin();
74 vit != mesh.finite_vertices_end();
80 newField[added++] = field[
count];
99 pointPoints.
setSize(mesh.vertexCount());
105 typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
106 vit != mesh.finite_vertices_end();
115 std::list<typename T::Vertex_handle> adjVerts;
116 mesh.finite_adjacent_vertices(vit, std::back_inserter(adjVerts));
122 typename std::list<typename T::Vertex_handle>::const_iterator
123 adjVertI = adjVerts.begin();
124 adjVertI != adjVerts.end();
128 typename T::Vertex_handle vh = *adjVertI;
134 globalIndexing.toGlobal(vh->procIndex(), vh->index())
139 pointPoints[vit->index()].
transfer(indices);
167 typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
168 vit != mesh.finite_vertices_end();
177 alignments[vit->index()] = vit->alignment();
195 typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
196 vit != mesh.finite_vertices_end();
205 points[vit->index()] =
topoint(vit->point());
216 const label maxRefinementIterations,
217 const scalar defaultCellSize
220 for (
label iter = 0; iter < maxRefinementIterations; ++iter)
226 CellSizeDelaunay::Finite_cells_iterator cit =
227 mesh.finite_cells_begin();
228 cit != mesh.finite_cells_end();
232 const point newPoint =
237 cit->vertex(0)->point(),
238 cit->vertex(1)->point(),
239 cit->vertex(2)->point(),
240 cit->vertex(3)->point()
244 if (geometryToConformTo.
inside(newPoint))
246 ptsToInsert.
append(newPoint);
267 int main(
int argc,
char *argv[])
274 label maxRefinementIterations = 2;
275 label maxSmoothingIterations = 200;
276 scalar minResidual = 0;
277 scalar defaultCellSize = 0.001;
278 scalar nearFeatDistSqrCoeff = 1
e-8;
306 "cvSearchableSurfaces",
313 foamyHexMeshDict.subDict(
"geometry"),
314 foamyHexMeshDict.lookupOrDefault(
"singleRegionName",
true)
322 foamyHexMeshDict.subDict(
"surfaceConformation")
335 foamyHexMeshDict.subDict(
"backgroundMeshDecomposition")
354 geometryToConformTo.
geometry()[surfI];
356 Info<<
nl <<
"Inserting points from surface " << surface.
name()
357 <<
" (" << surface.type() <<
")" <<
endl;
362 Info<<
" Number of points = " << points.size() <<
endl;
375 nearFeatDistSqrCoeff,
386 geometryToConformTo.
features()[infoFeature];
394 pointAlignment() += norms[nI];
397 pointAlignment().normalise();
398 pointAlignment().orthogonalise();
405 nearFeatDistSqrCoeff,
413 geometryToConformTo.
features()[infoFeature];
421 pointAlignment() += norms[nI];
424 pointAlignment().normalise();
425 pointAlignment().orthogonalise();
438 pointAlignment.
set(
new triad(normals[0]));
444 if (
bMesh().positionOnThisProcessor(points[pI]))
446 CellSizeDelaunay::Vertex_handle vh = mesh.
insert 457 CellSizeDelaunay::Vertex_handle vh = mesh.
insert 474 maxRefinementIterations,
489 triadField alignments(buildAlignmentField(mesh));
500 CellSizeDelaunay::Finite_vertices_iterator vit =
501 mesh.finite_vertices_begin();
502 vit != mesh.finite_vertices_end();
506 if (vit->nearBoundary())
508 fixedAlignments[vit->index()] = vit->alignment();
514 for (
label iter = 0; iter < maxSmoothingIterations; iter++)
516 Info<<
"Iteration " << iter;
518 meshDistributor().distribute(points);
519 meshDistributor().distribute(alignments);
527 const labelList& pPoints = pointPoints[pI];
534 const triad& oldTriad = alignments[pI];
535 triad& newTriad = triadAv[pI];
538 const triad& fixedAlignment = fixedAlignments[pI];
540 forAll(pPoints, adjPointi)
542 const label adjPointIndex = pPoints[adjPointi];
544 scalar dist =
mag(points[pI] - points[adjPointIndex]);
548 triad tmpTriad = alignments[adjPointIndex];
552 if (tmpTriad.
set(dir))
554 tmpTriad[dir] *= (1.0/dist);
558 newTriad += tmpTriad;
567 forAll(fixedAlignment, dirI)
569 if (fixedAlignment.
set(dirI))
577 forAll(fixedAlignment, dirI)
579 if (fixedAlignment.
set(dirI))
581 newTriad.
align(fixedAlignment[dirI]);
585 else if (nFixed == 2)
587 forAll(fixedAlignment, dirI)
589 if (fixedAlignment.
set(dirI))
591 newTriad[dirI] = fixedAlignment[dirI];
601 else if (nFixed == 3)
603 forAll(fixedAlignment, dirI)
605 if (fixedAlignment.
set(dirI))
607 newTriad[dirI] = fixedAlignment[dirI];
621 scalar dotProd = (oldTriad[dir] & newTriad[dir]);
622 scalar
diff =
mag(dotProd) - 1.0;
624 residual +=
mag(diff);
631 alignments[pI] = triadAv[pI].sortxyz();
636 Info<<
", Residual = " << residual <<
endl;
638 if (residual <= minResidual)
650 const triad& tri = alignments[pI];
673 filterFarPoints(mesh, points)
687 filterFarPoints(mesh, sizes)
700 filterFarPoints(mesh, alignments)
705 alignmentsIO.write();
const labelListList & featurePointNormals() const
Return the indices of the normals that are adjacent to the.
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.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
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 word & name() const
Return name.
bool set(const label) const
Is element set.
pointFromPoint topoint(const Point &P)
const labelListList & edgeNormals() const
Return the indices of the normals that are adjacent to the.
void align(const vector &v)
Align this triad with the given vector v.
T & ref() const
Return non-const reference or generate a fatal error.
void size(const label)
Override size to be inconsistent with allocated storage.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Ostream & endl(Ostream &os)
Add newline and flush stream.
label count(const ListType &l, typename ListType::const_reference x)
Count the number of occurrences of a value in a list.
Vertex_handle insert(const Foam::point &pt, const scalar &size, const triad &alignment, const Foam::indexedVertexEnum::vertexType type=Vb::vtInternal)
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
PointIndexHit< point > pointIndexHit
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
void resize(const label)
Alias for setSize(const label)
bool set(const direction d) const
Is the vector in the direction d set.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
void orthogonalise()
Orthogonalise this triad so that it is ortho-normal.
void printInfo(Ostream &os) const
Write mesh statistics to stream.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
vectorField pointField
pointField is a vectorField.
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const =0
Pre-declare SubField and related Field type.
static const word & geometryDir()
Return the geometry directory name.
const word & constant() const
Return constant name.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
bool hit() const
Is there a hit.
Container for searchableSurfaces.
void distribute(const backgroundMeshDecomposition &decomposition)
const word & system() const
Return system name.
Field< triad > triadField
Specialisation of Field<T> for triad.
virtual tmp< pointField > points() const =0
Get the points that define the surface.
void set(T *)
Set pointer to that given.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Representation of a 3D Cartesian coordinate system as a Vector of vectors.
void setSize(const label)
Reset size of List.
static bool & parRun()
Is this a parallel run?
time_t elapsedClockTime() const
Returns wall-clock time from clock instantiation.
Class containing processor-to-processor mapping information.
label vertexCount() const
Return the vertex count (the next unique vertex index)
double elapsedCpuTime() const
Return CPU time (in seconds) from the start.
vector point
Point is a vector.
void normalise()
Normalise each set axis vector to have a unit magnitude.
const polyBoundaryMesh & bMesh
dimensioned< scalar > mag(const dimensioned< Type > &)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Store a background polyMesh to use for the decomposition of space and queries for parallel conformalV...
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
const doubleScalar e
Elementary charge.
virtual bool write(const bool write=true) const
Write using setting from DB.
A class for managing temporary objects.
fileName path() const
Return path.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const =0
From a set of points and indices get the normal.
A primitive field of type <T> with automated input and output.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
label index() const
Return index.