222 CellSizeDelaunay::Finite_vertices_iterator vit =
223 finite_vertices_begin();
224 vit != finite_vertices_end();
228 std::list<Vertex_handle> verts;
229 adjacent_vertices(vit, std::back_inserter(verts));
231 bool removePt =
true;
234 std::list<Vertex_handle>::iterator aVit = verts.begin();
239 Vertex_handle avh = *aVit;
242 mag(avh->targetCellSize() - vit->targetCellSize())
243 /
max(vit->targetCellSize(), 1
e-6);
264 tmp<pointField> tcellCentres(
new pointField(number_of_finite_cells()));
270 CellSizeDelaunay::Finite_cells_iterator c = finite_cells_begin();
271 c != finite_cells_end();
275 if (c->hasFarPoint())
285 CGAL::centroid<baseK>
287 c->vertex(0)->point(),
288 c->vertex(1)->point(),
289 c->vertex(2)->point(),
290 c->vertex(3)->point()
294 cellCentres[count++] = centre;
297 cellCentres.resize(count);
307 "refinementTriangulation_" 314 Info<<
"Write refinementTriangulation" <<
endl;
318 CellSizeDelaunay::Finite_edges_iterator
e = finite_edges_begin();
319 e != finite_edges_end();
323 Cell_handle c =
e->first;
324 Vertex_handle vA = c->vertex(
e->second);
325 Vertex_handle vB = c->vertex(
e->third);
328 if (vA->farPoint() || vB->farPoint())
334 if (vA->referred() && vB->referred())
347 Info<<
" Triangulation is valid" <<
endl;
352 <<
"Triangulation is not valid" 369 if (this->vertexCount())
383 if (mesh.nPoints() == this->vertexCount())
405 mesh.time().timeName(),
416 sizes.size() == this->vertexCount()
417 && alignments.size() == this->vertexCount()
422 Finite_vertices_iterator vit = finite_vertices_begin();
423 vit != finite_vertices_end();
427 vit->targetCellSize() = sizes[vit->index()];
428 vit->alignment() = alignments[vit->index()];
434 <<
"Cell size point field is not the same size as the " 462 if (dimension() > 2 && !is_infinite(ch))
468 topoint(ch->vertex(0)->point()),
469 topoint(ch->vertex(1)->point()),
470 topoint(ch->vertex(2)->point()),
471 topoint(ch->vertex(3)->point())
474 bary = tet.pointToBarycentric(pt);
481 DynamicList<Foam::point> pts(number_of_vertices());
485 Finite_vertices_iterator vit = finite_vertices_begin();
486 vit != finite_vertices_end();
492 pts.append(
topoint(vit->point()));
504 const backgroundMeshDecomposition& decomposition
507 DynamicList<Foam::point>
points(number_of_vertices());
508 DynamicList<scalar> sizes(number_of_vertices());
509 DynamicList<tensor> alignments(number_of_vertices());
511 DynamicList<Vb> farPts(8);
515 Finite_vertices_iterator vit = finite_vertices_begin();
516 vit != finite_vertices_end();
523 sizes.append(vit->targetCellSize());
524 alignments.append(vit->alignment());
526 else if (vit->farPoint())
539 farPts.last().targetCellSize() = vit->targetCellSize();
540 farPts.last().alignment() = vit->alignment();
544 autoPtr<mapDistribute> mapDist =
551 mapDist().distribute(sizes);
552 mapDist().distribute(alignments);
559 DynamicList<Vb> verticesToInsert(
points.size());
564 verticesToInsert.append(farPts[ptI]);
570 verticesToInsert.append
581 verticesToInsert.last().targetCellSize() = sizes[pI];
582 verticesToInsert.last().alignment() = alignments[pI];
585 Info<<
nl <<
" Inserting distributed background tessellation..." <<
endl;
587 this->rangeInsertWithInfo
589 verticesToInsert.begin(),
590 verticesToInsert.end(),
594 sync(decomposition.procBounds());
596 Info<<
" Total number of vertices after redistribution " 608 Finite_vertices_iterator vit = finite_vertices_begin();
609 vit != finite_vertices_end();
613 alignmentsTmp[count++] = vit->alignment();
616 return alignmentsTmp;
622 Info<<
"Writing " << meshSubDir <<
endl;
628 Finite_cells_iterator cit = finite_cells_begin();
629 cit != finite_cells_end();
633 if (!cit->hasFarPoint() && !is_infinite(cit))
635 cit->cellIndex() = cellCount++;
648 const polyMesh& mesh =
meshPtr();
655 mesh.time().timeName(),
670 mesh.time().timeName(),
684 Finite_vertices_iterator vit = finite_vertices_begin();
685 vit != finite_vertices_end();
689 if (!vit->farPoint())
692 sizes[vertexMap[
labelPair(vit->index(), vit->procIndex())]] =
693 vit->targetCellSize();
695 alignments[vertexMap[
labelPair(vit->index(), vit->procIndex())]] =
732 const autoPtr<backgroundMeshDecomposition>& decomposition
737 scalar cellCount = 0;
741 Finite_cells_iterator cit = finite_cells_begin();
742 cit != finite_cells_end();
746 if (!cit->hasFarPoint() && !is_infinite(cit))
749 CGAL::Tetrahedron_3<baseK> tet
751 cit->vertex(0)->point(),
752 cit->vertex(1)->point(),
753 cit->vertex(2)->point(),
754 cit->vertex(3)->point()
762 && !decomposition().positionOnThisProcessor(centre)
768 scalar volume = CGAL::to_double(tet.volume());
770 scalar averagedPointCellSize = 0;
774 for (
label vI = 0; vI < 4; ++vI)
776 averagedPointCellSize += cit->vertex(vI)->targetCellSize();
780 averagedPointCellSize /= 4;
796 cellCount += volume/
pow(averagedPointCellSize, 3);
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
void writeTriangulation()
tetrahedron< point, const point & > tetPointRef
#define forAll(list, i)
Loop across all elements in list.
~cellShapeControlMesh()
Destructor.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
pointFromPoint topoint(const Point &P)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
tmp< pointField > cellCentres() const
Get the centres of all the tets.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
HashTable< label, labelPair, FixedList< label, 2 >::Hash<> > labelTolabelPairHashTable
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool distribute(const boundBox &bb)
void barycentricCoords(const Foam::point &pt, barycentric &bary, Cell_handle &ch) const
Calculate and return the barycentric coordinates for.
A bounding box defined in terms of the points at its extremities.
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
const dimensionedScalar & c
Speed of light in a vacuum.
tensorField dumpAlignments() const
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
static const pointMesh & New(const polyMesh &mesh)
vectorField pointField
pointField is a vectorField.
List< scalar > scalarList
A List of scalars.
void reset()
Clear the entire triangulation.
Pair< label > labelPair
Label pair.
PointFrompoint toPoint(const Foam::point &p)
List< label > labelList
A List of labels.
void distribute(const backgroundMeshDecomposition &decomposition)
errorManip< error > abort(error &err)
cellShapeControlMesh(const Time &runTime)
defineTypeNameAndDebug(combustionModel, 0)
word name(const complex &)
Return a string representation of a complex.
CGAL::Delaunay_triangulation_3< K, Tds, FastLocator > CellSizeDelaunay
static word meshSubDir
Return the mesh sub-directory name (usually "cellShapeControlMesh")
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
CellSizeDelaunay::Cell_handle Cell_handle
static bool & parRun()
Is this a parallel run?
dimensioned< scalar > mag(const dimensioned< Type > &)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const doubleScalar e
Elementary charge.
label estimateCellCount(const autoPtr< backgroundMeshDecomposition > &decomposition) const
A class for managing temporary objects.
IOField< triad > triadIOField
triadField with IO.
autoPtr< polyMesh > createMesh(const fileName &name, labelTolabelPairHashTable &vertexMap, labelList &cellMap, const bool writeDelaunayData=true) const
Create an fvMesh from the triangulation.