All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cellShapeControlMesh.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2012-2019 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "cellShapeControlMesh.H"
28 #include "pointIOField.H"
29 #include "scalarIOField.H"
30 #include "triadIOField.H"
31 #include "tetPointRef.H"
32 #include "plane.H"
33 #include "transform.H"
34 #include "meshTools.H"
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40 defineTypeNameAndDebug(cellShapeControlMesh, 0);
41 
42 word cellShapeControlMesh::meshSubDir = "cellShapeControlMesh";
43 }
44 
45 
46 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
47 
48 //Foam::tensor Foam::cellShapeControlMesh::requiredAlignment
49 //(
50 // const Foam::point& pt,
51 // const searchableSurfaces& allGeometry,
52 // const conformationSurfaces& geometryToConformTo
53 //) const
54 //{
55 // pointIndexHit surfHit;
56 // label hitSurface;
57 //
58 // geometryToConformTo.findSurfaceNearest
59 // (
60 // pt,
61 // sqr(great),
62 // surfHit,
63 // hitSurface
64 // );
65 //
66 // if (!surfHit.hit())
67 // {
68 // FatalErrorInFunction
69 // << "findSurfaceNearest did not find a hit across the surfaces."
70 // << exit(FatalError) << endl;
71 // }
72 //
73 // // Primary alignment
74 //
75 // vectorField norm(1);
76 //
77 // allGeometry[hitSurface].getNormal
78 // (
79 // List<pointIndexHit>(1, surfHit),
80 // norm
81 // );
82 //
83 // const vector np = norm[0];
84 //
85 // // Generate equally spaced 'spokes' in a circle normal to the
86 // // direction from the vertex to the closest point on the surface
87 // // and look for a secondary intersection.
88 //
89 // const vector d = surfHit.hitPoint() - pt;
90 //
91 // const tensor Rp = rotationTensor(vector(0,0,1), np);
92 //
93 // const label s = 36;//foamyHexMeshControls().alignmentSearchSpokes();
94 //
95 // scalar closestSpokeHitDistance = great;
96 //
97 // pointIndexHit closestSpokeHit;
98 //
99 // label closestSpokeSurface = -1;
100 //
101 // const scalar spanMag = geometryToConformTo.globalBounds().mag();
102 //
103 // for (label i = 0; i < s; i++)
104 // {
105 // vector spoke
106 // (
107 // Foam::cos(i*constant::mathematical::twoPi/s),
108 // Foam::sin(i*constant::mathematical::twoPi/s),
109 // 0
110 // );
111 //
112 // spoke *= spanMag;
113 //
114 // spoke = Rp & spoke;
115 //
116 // pointIndexHit spokeHit;
117 //
118 // label spokeSurface = -1;
119 //
120 // // internal spoke
121 //
122 // geometryToConformTo.findSurfaceNearestIntersection
123 // (
124 // pt,
125 // pt + spoke,
126 // spokeHit,
127 // spokeSurface
128 // );
129 //
130 // if (spokeHit.hit())
131 // {
132 // scalar spokeHitDistance = mag
133 // (
134 // spokeHit.hitPoint() - pt
135 // );
136 //
137 // if (spokeHitDistance < closestSpokeHitDistance)
138 // {
139 // closestSpokeHit = spokeHit;
140 // closestSpokeSurface = spokeSurface;
141 // closestSpokeHitDistance = spokeHitDistance;
142 // }
143 // }
144 //
145 // // external spoke
146 //
147 // Foam::point mirrorPt = pt + 2*d;
148 //
149 // geometryToConformTo.findSurfaceNearestIntersection
150 // (
151 // mirrorPt,
152 // mirrorPt + spoke,
153 // spokeHit,
154 // spokeSurface
155 // );
156 //
157 // if (spokeHit.hit())
158 // {
159 // scalar spokeHitDistance = mag
160 // (
161 // spokeHit.hitPoint() - mirrorPt
162 // );
163 //
164 // if (spokeHitDistance < closestSpokeHitDistance)
165 // {
166 // closestSpokeHit = spokeHit;
167 // closestSpokeSurface = spokeSurface;
168 // closestSpokeHitDistance = spokeHitDistance;
169 // }
170 // }
171 // }
172 //
173 // if (closestSpokeSurface == -1)
174 // {
180 //
181 // return I;
182 // }
183 //
184 // // Auxiliary alignment generated by spoke intersection normal.
185 //
186 // allGeometry[closestSpokeSurface].getNormal
187 // (
188 // List<pointIndexHit>(1, closestSpokeHit),
189 // norm
190 // );
191 //
192 // const vector& na = norm[0];
193 //
194 // // Secondary alignment
195 // vector ns = np ^ na;
196 //
197 // if (mag(ns) < small)
198 // {
199 // FatalErrorInFunction
200 // << "Parallel normals detected in spoke search." << nl
201 // << "point: " << pt << nl
202 // << "closest surface point: " << surfHit.hitPoint() << nl
203 // << "closest spoke hit: " << closestSpokeHit.hitPoint() << nl
204 // << "np: " << surfHit.hitPoint() + np << nl
205 // << "ns: " << closestSpokeHit.hitPoint() + na << nl
206 // << exit(FatalError);
207 // }
208 //
209 // ns /= mag(ns);
210 //
211 // tensor Rs = rotationTensor((Rp & vector(0,1,0)), ns);
212 //
213 // return (Rs & Rp);
214 //}
215 
216 
218 {
219  label nRemoved = 0;
220  for
221  (
222  CellSizeDelaunay::Finite_vertices_iterator vit =
223  finite_vertices_begin();
224  vit != finite_vertices_end();
225  ++vit
226  )
227  {
228  std::list<Vertex_handle> verts;
229  adjacent_vertices(vit, std::back_inserter(verts));
230 
231  bool removePt = true;
232  for
233  (
234  std::list<Vertex_handle>::iterator aVit = verts.begin();
235  aVit != verts.end();
236  ++aVit
237  )
238  {
239  Vertex_handle avh = *aVit;
240 
241  scalar diff =
242  mag(avh->targetCellSize() - vit->targetCellSize())
243  /max(vit->targetCellSize(), 1e-6);
244 
245  if (diff > 0.05)
246  {
247  removePt = false;
248  }
249  }
250 
251  if (removePt)
252  {
253  remove(vit);
254  nRemoved++;
255  }
256  }
257 
258  return nRemoved;
259 }
260 
261 
263 {
264  tmp<pointField> tcellCentres(new pointField(number_of_finite_cells()));
265  pointField& cellCentres = tcellCentres.ref();
266 
267  label count = 0;
268  for
269  (
270  CellSizeDelaunay::Finite_cells_iterator c = finite_cells_begin();
271  c != finite_cells_end();
272  ++c
273  )
274  {
275  if (c->hasFarPoint())
276  {
277  continue;
278  }
279 
280  scalarList bary;
282 
283  const Foam::point centre = topoint
284  (
285  CGAL::centroid<baseK>
286  (
287  c->vertex(0)->point(),
288  c->vertex(1)->point(),
289  c->vertex(2)->point(),
290  c->vertex(3)->point()
291  )
292  );
293 
294  cellCentres[count++] = centre;
295  }
296 
297  cellCentres.resize(count);
298 
299  return tcellCentres;
300 }
301 
302 
304 {
305  OFstream str
306  (
307  "refinementTriangulation_"
309  + ".obj"
310  );
311 
312  label count = 0;
313 
314  Info<< "Write refinementTriangulation" << endl;
315 
316  for
317  (
318  CellSizeDelaunay::Finite_edges_iterator e = finite_edges_begin();
319  e != finite_edges_end();
320  ++e
321  )
322  {
323  Cell_handle c = e->first;
324  Vertex_handle vA = c->vertex(e->second);
325  Vertex_handle vB = c->vertex(e->third);
326 
327  // Don't write far edges
328  if (vA->farPoint() || vB->farPoint())
329  {
330  continue;
331  }
332 
333  // Don't write unowned edges
334  if (vA->referred() && vB->referred())
335  {
336  continue;
337  }
338 
339  pointFromPoint p1 = topoint(vA->point());
340  pointFromPoint p2 = topoint(vB->point());
341 
342  meshTools::writeOBJ(str, p1, p2, count);
343  }
344 
345  if (is_valid())
346  {
347  Info<< " Triangulation is valid" << endl;
348  }
349  else
350  {
352  << "Triangulation is not valid"
353  << abort(FatalError);
354  }
355 }
356 
357 
358 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
359 
361 :
362  DistributedDelaunayMesh<CellSizeDelaunay>
363  (
364  runTime,
365  meshSubDir
366  ),
367  runTime_(runTime)
368 {
369  if (this->vertexCount())
370  {
371  fvMesh mesh
372  (
373  IOobject
374  (
375  meshSubDir,
376  runTime.timeName(),
377  runTime,
380  )
381  );
382 
383  if (mesh.nPoints() == this->vertexCount())
384  {
385  pointScalarField sizes
386  (
387  IOobject
388  (
389  "sizes",
390  runTime.timeName(),
391  meshSubDir,
392  runTime,
395  false
396  ),
397  pointMesh::New(mesh)
398  );
399 
400  triadIOField alignments
401  (
402  IOobject
403  (
404  "alignments",
405  mesh.time().timeName(),
406  meshSubDir,
407  mesh.time(),
410  false
411  )
412  );
413 
414  if
415  (
416  sizes.size() == this->vertexCount()
417  && alignments.size() == this->vertexCount()
418  )
419  {
420  for
421  (
422  Finite_vertices_iterator vit = finite_vertices_begin();
423  vit != finite_vertices_end();
424  ++vit
425  )
426  {
427  vit->targetCellSize() = sizes[vit->index()];
428  vit->alignment() = alignments[vit->index()];
429  }
430  }
431  else
432  {
434  << "Cell size point field is not the same size as the "
435  << "mesh."
436  << abort(FatalError);
437  }
438  }
439  }
440 }
441 
442 
443 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
444 
446 {}
447 
448 
449 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
450 
452 (
453  const Foam::point& pt,
454  barycentric& bary,
455  Cell_handle& ch
456 ) const
457 {
458  // Use the previous cell handle as a hint on where to start searching
459  // Giving a hint causes strange errors...
460  ch = locate(toPoint(pt));
461 
462  if (dimension() > 2 && !is_infinite(ch))
463  {
464  oldCellHandle_ = ch;
465 
466  tetPointRef tet
467  (
468  topoint(ch->vertex(0)->point()),
469  topoint(ch->vertex(1)->point()),
470  topoint(ch->vertex(2)->point()),
471  topoint(ch->vertex(3)->point())
472  );
473 
474  bary = tet.pointToBarycentric(pt);
475  }
476 }
477 
478 
480 {
481  DynamicList<Foam::point> pts(number_of_vertices());
482 
483  for
484  (
485  Finite_vertices_iterator vit = finite_vertices_begin();
486  vit != finite_vertices_end();
487  ++vit
488  )
489  {
490  if (vit->real())
491  {
492  pts.append(topoint(vit->point()));
493  }
494  }
495 
496  boundBox bb(pts);
497 
498  return bb;
499 }
500 
501 
503 (
504  const backgroundMeshDecomposition& decomposition
505 )
506 {
507  DynamicList<Foam::point> points(number_of_vertices());
508  DynamicList<scalar> sizes(number_of_vertices());
509  DynamicList<tensor> alignments(number_of_vertices());
510 
511  DynamicList<Vb> farPts(8);
512 
513  for
514  (
515  Finite_vertices_iterator vit = finite_vertices_begin();
516  vit != finite_vertices_end();
517  ++vit
518  )
519  {
520  if (vit->real())
521  {
522  points.append(topoint(vit->point()));
523  sizes.append(vit->targetCellSize());
524  alignments.append(vit->alignment());
525  }
526  else if (vit->farPoint())
527  {
528  farPts.append
529  (
530  Vb
531  (
532  vit->point(),
533  -1,
534  Vb::vtFar,
536  )
537  );
538 
539  farPts.last().targetCellSize() = vit->targetCellSize();
540  farPts.last().alignment() = vit->alignment();
541  }
542  }
543 
544  autoPtr<mapDistribute> mapDist =
546  (
547  decomposition,
548  points
549  );
550 
551  mapDist().distribute(sizes);
552  mapDist().distribute(alignments);
553 
554  // Reset the entire tessellation
556 
557 
558  // Internal points have to be inserted first
559  DynamicList<Vb> verticesToInsert(points.size());
560 
561 
562  forAll(farPts, ptI)
563  {
564  verticesToInsert.append(farPts[ptI]);
565  }
566 
567 
568  forAll(points, pI)
569  {
570  verticesToInsert.append
571  (
572  Vb
573  (
574  toPoint(points[pI]),
575  -1,
578  )
579  );
580 
581  verticesToInsert.last().targetCellSize() = sizes[pI];
582  verticesToInsert.last().alignment() = alignments[pI];
583  }
584 
585  Info<< nl << " Inserting distributed background tessellation..." << endl;
586 
587  this->rangeInsertWithInfo
588  (
589  verticesToInsert.begin(),
590  verticesToInsert.end(),
591  true
592  );
593 
594  sync(decomposition.procBounds());
595 
596  Info<< " Total number of vertices after redistribution "
597  << returnReduce(label(number_of_vertices()), sumOp<label>()) << endl;
598 }
599 
600 
602 {
603  tensorField alignmentsTmp(number_of_vertices(), Zero);
604 
605  label count = 0;
606  for
607  (
608  Finite_vertices_iterator vit = finite_vertices_begin();
609  vit != finite_vertices_end();
610  ++vit
611  )
612  {
613  alignmentsTmp[count++] = vit->alignment();
614  }
615 
616  return alignmentsTmp;
617 }
618 
619 
621 {
622  Info<< "Writing " << meshSubDir << endl;
623 
624  // Reindex the cells
625  label cellCount = 0;
626  for
627  (
628  Finite_cells_iterator cit = finite_cells_begin();
629  cit != finite_cells_end();
630  ++cit
631  )
632  {
633  if (!cit->hasFarPoint() && !is_infinite(cit))
634  {
635  cit->cellIndex() = cellCount++;
636  }
637  }
638 
640  labelList cellMap;
641 
642  autoPtr<polyMesh> meshPtr = DelaunayMesh<CellSizeDelaunay>::createMesh
643  (
644  meshSubDir,
645  vertexMap,
646  cellMap
647  );
648  const polyMesh& mesh = meshPtr();
649 
650  pointScalarField sizes
651  (
652  IOobject
653  (
654  "sizes",
655  mesh.time().timeName(),
656  meshSubDir,
657  mesh.time(),
660  ),
661  pointMesh::New(mesh),
662  scalar(0)
663  );
664 
665  triadIOField alignments
666  (
667  IOobject
668  (
669  "alignments",
670  mesh.time().timeName(),
671  meshSubDir,
672  mesh.time(),
675  ),
676  sizes.size()
677  );
678 
679  // Write alignments
680 // OFstream str(runTime_.path()/"alignments.obj");
681 
682  for
683  (
684  Finite_vertices_iterator vit = finite_vertices_begin();
685  vit != finite_vertices_end();
686  ++vit
687  )
688  {
689  if (!vit->farPoint())
690  {
691  // Populate sizes
692  sizes[vertexMap[labelPair(vit->index(), vit->procIndex())]] =
693  vit->targetCellSize();
694 
695  alignments[vertexMap[labelPair(vit->index(), vit->procIndex())]] =
696  vit->alignment();
697 
698 // // Write alignments
699 // const tensor& alignment = vit->alignment();
700 // pointFromPoint pt = topoint(vit->point());
701 //
702 // if
703 // (
704 // alignment.x() == triad::unset[0]
705 // || alignment.y() == triad::unset[0]
706 // || alignment.z() == triad::unset[0]
707 // )
708 // {
709 // Info<< "Bad alignment = " << vit->info();
710 //
711 // vit->alignment() = tensor::I;
712 //
713 // Info<< "New alignment = " << vit->info();
714 //
715 // continue;
716 // }
717 //
718 // meshTools::writeOBJ(str, pt, alignment.x() + pt);
719 // meshTools::writeOBJ(str, pt, alignment.y() + pt);
720 // meshTools::writeOBJ(str, pt, alignment.z() + pt);
721  }
722  }
723 
724  mesh.write();
725  sizes.write();
726  alignments.write();
727 }
728 
729 
731 (
732  const autoPtr<backgroundMeshDecomposition>& decomposition
733 ) const
734 {
735  // Loop over all the tets and estimate the cell count in each one
736 
737  scalar cellCount = 0;
738 
739  for
740  (
741  Finite_cells_iterator cit = finite_cells_begin();
742  cit != finite_cells_end();
743  ++cit
744  )
745  {
746  if (!cit->hasFarPoint() && !is_infinite(cit))
747  {
748  // TODO: Check if tet centre is on the processor..
749  CGAL::Tetrahedron_3<baseK> tet
750  (
751  cit->vertex(0)->point(),
752  cit->vertex(1)->point(),
753  cit->vertex(2)->point(),
754  cit->vertex(3)->point()
755  );
756 
757  pointFromPoint centre = topoint(CGAL::centroid(tet));
758 
759  if
760  (
762  && !decomposition().positionOnThisProcessor(centre)
763  )
764  {
765  continue;
766  }
767 
768  scalar volume = CGAL::to_double(tet.volume());
769 
770  scalar averagedPointCellSize = 0;
771  // scalar averagedPointCellSize = 1;
772 
773  // Get an average volume by averaging the cell size of the vertices
774  for (label vI = 0; vI < 4; ++vI)
775  {
776  averagedPointCellSize += cit->vertex(vI)->targetCellSize();
777  // averagedPointCellSize *= cit->vertex(vI)->targetCellSize();
778  }
779 
780  averagedPointCellSize /= 4;
781  // averagedPointCellSize = ::sqrt(averagedPointCellSize);
782 
783 // if (averagedPointCellSize < small)
784 // {
785 // Pout<< "Volume = " << volume << endl;
786 //
787 // for (label vI = 0; vI < 4; ++vI)
788 // {
789 // Pout<< "Point " << vI
790 // << ", point = " << topoint(cit->vertex(vI)->point())
791 // << ", size = " << cit->vertex(vI)->targetCellSize()
792 // << endl;
793 // }
794 // }
795 
796  cellCount += volume/pow(averagedPointCellSize, 3);
797  }
798  }
799 
800  return cellCount;
801 }
802 
803 
804 // ************************************************************************* //
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
Definition: pointFields.H:49
tetrahedron< point, const point & > tetPointRef
Definition: tetPointRef.H:44
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
~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.
Definition: label.H:59
pointFromPoint topoint(const Point &P)
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
tmp< pointField > cellCentres() const
Get the centres of all the tets.
static fvMesh * meshPtr
Definition: globalFoam.H:52
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
HashTable< label, labelPair, FixedList< label, 2 >::Hash<> > labelTolabelPairHashTable
Definition: DelaunayMesh.H:90
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
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.
Definition: boundBox.H:58
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
Definition: barycentric.H:45
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...
Definition: indexedVertex.H:51
static const pointMesh & New(const polyMesh &mesh)
Definition: MeshObject.C:44
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
dynamicFvMesh & mesh
const pointField & points
3D tensor transformation operations.
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
void reset()
Clear the entire triangulation.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
PointFrompoint toPoint(const Foam::point &p)
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
void distribute(const backgroundMeshDecomposition &decomposition)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
cellShapeControlMesh(const Time &runTime)
static const char nl
Definition: Ostream.H:260
defineTypeNameAndDebug(combustionModel, 0)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
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?
Definition: UPstream.H:399
boundBox bounds() const
messageStream Info
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.
Definition: doubleScalar.H:105
label estimateCellCount(const autoPtr< backgroundMeshDecomposition > &decomposition) const
A class for managing temporary objects.
Definition: PtrList.H:53
IOField< triad > triadIOField
triadField with IO.
Definition: triadIOField.H:42
autoPtr< polyMesh > createMesh(const fileName &name, labelTolabelPairHashTable &vertexMap, labelList &cellMap, const bool writeDelaunayData=true) const
Create an fvMesh from the triangulation.
Namespace for OpenFOAM.