36#include "indexedVertex.H"
60template<
class Triangulation,
class Type>
63 const Triangulation&
mesh,
76 mesh.finite_vertices_begin();
77 vit !=
mesh.finite_vertices_end();
108 typename T::Finite_vertices_iterator vit =
mesh.finite_vertices_begin();
109 vit !=
mesh.finite_vertices_end();
118 std::list<typename T::Vertex_handle> adjVerts;
119 mesh.finite_adjacent_vertices(vit, std::back_inserter(adjVerts));
125 typename std::list<typename T::Vertex_handle>::const_iterator
126 adjVertI = adjVerts.begin();
127 adjVertI != adjVerts.end();
131 typename T::Vertex_handle vh = *adjVertI;
137 globalIndexing.toGlobal(vh->procIndex(), vh->index())
142 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();
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);
267int 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.getOrDefault(
"singleRegionName",
true)
322 foamyHexMeshDict.subDict(
"surfaceConformation")
326 if (Pstream::parRun())
335 foamyHexMeshDict.subDict(
"backgroundMeshDecomposition")
351 const label surfI = geometryToConformTo.
surfaces()[sI];
354 geometryToConformTo.
geometry()[surfI];
356 Info<<
nl <<
"Inserting points from surface " <<
surface.name()
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();
433 surface.findNearest(ptField, distField, infoList);
436 surface.getNormal(infoList, normals);
438 pointAlignment.
set(
new triad(normals[0]));
442 if (Pstream::parRun())
446 CellSizeDelaunay::Vertex_handle vh =
mesh.insert
457 CellSizeDelaunay::Vertex_handle vh =
mesh.insert
474 maxRefinementIterations,
479 if (Pstream::parRun())
492 mesh.printInfo(Info);
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];
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];
595 newTriad[dirI] = triad::unset[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];
656 meshTools::writeOBJ(str,
points[pI], tri[dirI] +
points[pI]);
687 filterFarPoints(
mesh, sizes)
700 filterFarPoints(
mesh, alignments)
705 alignmentsIO.write();
708 runTime.printExecutionTime(Info);
reduce(hasMovingMesh, orOp< bool >())
Triangulation::Finite_vertices_iterator Finite_vertices_iterator
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
void append(const T &val)
Copy append an element to the end of this list.
Generic templated field type.
A primitive field of type <T> with automated input and output.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Defines the attributes of an object for which implicit objectRegistry management is supported,...
void transfer(List< T > &list)
void setSize(const label n)
Alias for resize()
void resize(const label len)
Adjust allocated size of list.
Output to file stream, using an OSstream.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
label index() const noexcept
Return the hit index.
bool hit() const noexcept
Is there a hit?
A list of faces which address into the list of points.
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
void size(const label n)
Older name for setAddressableSize.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
void set(T *p) noexcept
Deprecated(2018-02) Identical to reset().
Store a background polyMesh to use for the decomposition of space and queries for parallel conformalV...
const labelListList & edgeNormals() const
Return the indices of the normals that are adjacent to the.
const labelListList & featurePointNormals() const
Return the indices of the normals that are adjacent to the.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
Container for searchableSurfaces. The collection is specified as a dictionary. For example,...
A class for managing temporary objects.
Representation of a 3D Cartesian coordinate system as a Vector of row vectors.
void align(const vector &v)
Align this triad with the given vector v.
bool set(const direction d) const
Is the vector in the direction d set.
void normalize()
Same as normalise.
void orthogonalize()
Same as orthogonalise.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
const wordList surface
Standard surface field types (scalar, vector, tensor, etc)
pointFromPoint topoint(const Point &P)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.