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().normalize();
398 pointAlignment().orthogonalize();
405 nearFeatDistSqrCoeff,
413 geometryToConformTo.
features()[infoFeature];
421 pointAlignment() += norms[nI];
424 pointAlignment().normalize();
425 pointAlignment().orthogonalize();
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)
646 OFstream str(runTime.path()/
"alignments.obj");
650 const triad& tri = alignments[pI];
673 filterFarPoints(mesh, points)
687 filterFarPoints(mesh, sizes)
700 filterFarPoints(mesh, alignments)
705 alignmentsIO.write();
707 Info<<
nl <<
"ExecutionTime = " << runTime.elapsedCpuTime() <<
" s" 708 <<
" ClockTime = " << runTime.elapsedClockTime() <<
" s"
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.
pointFromPoint topoint(const Point &P)
const labelListList & edgeNormals() const
Return the indices of the normals that are adjacent to the.
const double e
Elementary charge.
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.
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 orthogonalize()
Orthogonalize this triad so that it is ortho-normal.
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 printInfo(Ostream &os) const
Write mesh statistics to stream.
void normalize()
Normalize each set axis vector to have a unit magnitude.
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.
cachedRandom rndGen(label(0), -1)
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)
Simple random number generator.
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?
Class containing processor-to-processor mapping information.
label vertexCount() const
Return the vertex count (the next unique vertex index)
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)
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...
A class for managing temporary objects.
virtual bool write(const bool valid=true) const
Write using setting from DB.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
PrimitivePatch< face, List, const pointField > bMesh
Holder of faceList and points. (v.s. e.g. primitivePatch which references points) ...
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.