cellShapeControlMesh.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2012-2016 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 "tetrahedron.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 
360 Foam::cellShapeControlMesh::cellShapeControlMesh(const Time& runTime)
361 :
362  DistributedDelaunayMesh<CellSizeDelaunay>
363  (
364  runTime,
365  meshSubDir
366  ),
367  runTime_(runTime),
368  defaultCellSize_(0.0)
369 {
370  if (this->vertexCount())
371  {
372  fvMesh mesh
373  (
374  IOobject
375  (
376  meshSubDir,
377  runTime.timeName(),
378  runTime,
381  )
382  );
383 
384  if (mesh.nPoints() == this->vertexCount())
385  {
386  pointScalarField sizes
387  (
388  IOobject
389  (
390  "sizes",
391  runTime.timeName(),
392  meshSubDir,
393  runTime,
396  false
397  ),
398  pointMesh::New(mesh)
399  );
400 
401  triadIOField alignments
402  (
403  IOobject
404  (
405  "alignments",
406  mesh.time().timeName(),
407  meshSubDir,
408  mesh.time(),
411  false
412  )
413  );
414 
415  if
416  (
417  sizes.size() == this->vertexCount()
418  && alignments.size() == this->vertexCount()
419  )
420  {
421  for
422  (
423  Finite_vertices_iterator vit = finite_vertices_begin();
424  vit != finite_vertices_end();
425  ++vit
426  )
427  {
428  vit->targetCellSize() = sizes[vit->index()];
429  vit->alignment() = alignments[vit->index()];
430  }
431  }
432  else
433  {
435  << "Cell size point field is not the same size as the "
436  << "mesh."
437  << abort(FatalError);
438  }
439  }
440  }
441 }
442 
443 
444 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
445 
447 {}
448 
449 
450 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
451 
453 (
454  const Foam::point& pt,
455  scalarList& bary,
456  Cell_handle& ch
457 ) const
458 {
459  // Use the previous cell handle as a hint on where to start searching
460  // Giving a hint causes strange errors...
461  ch = locate(toPoint(pt));
462 
463  if (dimension() > 2 && !is_infinite(ch))
464  {
465  oldCellHandle_ = ch;
466 
467  tetPointRef tet
468  (
469  topoint(ch->vertex(0)->point()),
470  topoint(ch->vertex(1)->point()),
471  topoint(ch->vertex(2)->point()),
472  topoint(ch->vertex(3)->point())
473  );
474 
475  tet.barycentric(pt, bary);
476  }
477 }
478 
479 
481 {
482  DynamicList<Foam::point> pts(number_of_vertices());
483 
484  for
485  (
486  Finite_vertices_iterator vit = finite_vertices_begin();
487  vit != finite_vertices_end();
488  ++vit
489  )
490  {
491  if (vit->real())
492  {
493  pts.append(topoint(vit->point()));
494  }
495  }
496 
497  boundBox bb(pts);
498 
499  return bb;
500 }
501 
502 
504 (
505  const backgroundMeshDecomposition& decomposition
506 )
507 {
508  DynamicList<Foam::point> points(number_of_vertices());
509  DynamicList<scalar> sizes(number_of_vertices());
510  DynamicList<tensor> alignments(number_of_vertices());
511 
512  DynamicList<Vb> farPts(8);
513 
514  for
515  (
516  Finite_vertices_iterator vit = finite_vertices_begin();
517  vit != finite_vertices_end();
518  ++vit
519  )
520  {
521  if (vit->real())
522  {
523  points.append(topoint(vit->point()));
524  sizes.append(vit->targetCellSize());
525  alignments.append(vit->alignment());
526  }
527  else if (vit->farPoint())
528  {
529  farPts.append
530  (
531  Vb
532  (
533  vit->point(),
534  -1,
535  Vb::vtFar,
537  )
538  );
539 
540  farPts.last().targetCellSize() = vit->targetCellSize();
541  farPts.last().alignment() = vit->alignment();
542  }
543  }
544 
545  autoPtr<mapDistribute> mapDist =
547  (
548  decomposition,
549  points
550  );
551 
552  mapDist().distribute(sizes);
553  mapDist().distribute(alignments);
554 
555  // Reset the entire tessellation
557 
558 
559  // Internal points have to be inserted first
560  DynamicList<Vb> verticesToInsert(points.size());
561 
562 
563  forAll(farPts, ptI)
564  {
565  verticesToInsert.append(farPts[ptI]);
566  }
567 
568 
569  forAll(points, pI)
570  {
571  verticesToInsert.append
572  (
573  Vb
574  (
575  toPoint(points[pI]),
576  -1,
579  )
580  );
581 
582  verticesToInsert.last().targetCellSize() = sizes[pI];
583  verticesToInsert.last().alignment() = alignments[pI];
584  }
585 
586  Info<< nl << " Inserting distributed background tessellation..." << endl;
587 
588  this->rangeInsertWithInfo
589  (
590  verticesToInsert.begin(),
591  verticesToInsert.end(),
592  true
593  );
594 
595  sync(decomposition.procBounds());
596 
597  Info<< " Total number of vertices after redistribution "
598  << returnReduce(label(number_of_vertices()), sumOp<label>()) << endl;
599 }
600 
601 
603 {
604  tensorField alignmentsTmp(number_of_vertices(), Zero);
605 
606  label count = 0;
607  for
608  (
609  Finite_vertices_iterator vit = finite_vertices_begin();
610  vit != finite_vertices_end();
611  ++vit
612  )
613  {
614  alignmentsTmp[count++] = vit->alignment();
615  }
616 
617  return alignmentsTmp;
618 }
619 
620 
622 {
623  Info<< "Writing " << meshSubDir << endl;
624 
625  // Reindex the cells
626  label cellCount = 0;
627  for
628  (
629  Finite_cells_iterator cit = finite_cells_begin();
630  cit != finite_cells_end();
631  ++cit
632  )
633  {
634  if (!cit->hasFarPoint() && !is_infinite(cit))
635  {
636  cit->cellIndex() = cellCount++;
637  }
638  }
639 
641  labelList cellMap;
642 
643  autoPtr<polyMesh> meshPtr = DelaunayMesh<CellSizeDelaunay>::createMesh
644  (
645  meshSubDir,
646  vertexMap,
647  cellMap
648  );
649  const polyMesh& mesh = meshPtr();
650 
651  pointScalarField sizes
652  (
653  IOobject
654  (
655  "sizes",
656  mesh.time().timeName(),
657  meshSubDir,
658  mesh.time(),
661  ),
662  pointMesh::New(mesh),
663  scalar(0)
664  );
665 
666  triadIOField alignments
667  (
668  IOobject
669  (
670  "alignments",
671  mesh.time().timeName(),
672  meshSubDir,
673  mesh.time(),
676  ),
677  sizes.size()
678  );
679 
680  // Write alignments
681 // OFstream str(runTime_.path()/"alignments.obj");
682 
683  for
684  (
685  Finite_vertices_iterator vit = finite_vertices_begin();
686  vit != finite_vertices_end();
687  ++vit
688  )
689  {
690  if (!vit->farPoint())
691  {
692  // Populate sizes
693  sizes[vertexMap[labelPair(vit->index(), vit->procIndex())]] =
694  vit->targetCellSize();
695 
696  alignments[vertexMap[labelPair(vit->index(), vit->procIndex())]] =
697  vit->alignment();
698 
699 // // Write alignments
700 // const tensor& alignment = vit->alignment();
701 // pointFromPoint pt = topoint(vit->point());
702 //
703 // if
704 // (
705 // alignment.x() == triad::unset[0]
706 // || alignment.y() == triad::unset[0]
707 // || alignment.z() == triad::unset[0]
708 // )
709 // {
710 // Info<< "Bad alignment = " << vit->info();
711 //
712 // vit->alignment() = tensor::I;
713 //
714 // Info<< "New alignment = " << vit->info();
715 //
716 // continue;
717 // }
718 //
719 // meshTools::writeOBJ(str, pt, alignment.x() + pt);
720 // meshTools::writeOBJ(str, pt, alignment.y() + pt);
721 // meshTools::writeOBJ(str, pt, alignment.z() + pt);
722  }
723  }
724 
725  mesh.write();
726  sizes.write();
727  alignments.write();
728 }
729 
730 
732 (
733  const autoPtr<backgroundMeshDecomposition>& decomposition
734 ) const
735 {
736  // Loop over all the tets and estimate the cell count in each one
737 
738  scalar cellCount = 0;
739 
740  for
741  (
742  Finite_cells_iterator cit = finite_cells_begin();
743  cit != finite_cells_end();
744  ++cit
745  )
746  {
747  if (!cit->hasFarPoint() && !is_infinite(cit))
748  {
749  // TODO: Check if tet centre is on the processor..
750  CGAL::Tetrahedron_3<baseK> tet
751  (
752  cit->vertex(0)->point(),
753  cit->vertex(1)->point(),
754  cit->vertex(2)->point(),
755  cit->vertex(3)->point()
756  );
757 
758  pointFromPoint centre = topoint(CGAL::centroid(tet));
759 
760  if
761  (
763  && !decomposition().positionOnThisProcessor(centre)
764  )
765  {
766  continue;
767  }
768 
769  scalar volume = CGAL::to_double(tet.volume());
770 
771  scalar averagedPointCellSize = 0;
772  //scalar averagedPointCellSize = 1;
773 
774  // Get an average volume by averaging the cell size of the vertices
775  for (label vI = 0; vI < 4; ++vI)
776  {
777  averagedPointCellSize += cit->vertex(vI)->targetCellSize();
778  //averagedPointCellSize *= cit->vertex(vI)->targetCellSize();
779  }
780 
781  averagedPointCellSize /= 4;
782  //averagedPointCellSize = ::sqrt(averagedPointCellSize);
783 
784 // if (averagedPointCellSize < SMALL)
785 // {
786 // Pout<< "Volume = " << volume << endl;
787 //
788 // for (label vI = 0; vI < 4; ++vI)
789 // {
790 // Pout<< "Point " << vI
791 // << ", point = " << topoint(cit->vertex(vI)->point())
792 // << ", size = " << cit->vertex(vI)->targetCellSize()
793 // << endl;
794 // }
795 // }
796 
797  cellCount += volume/pow(averagedPointCellSize, 3);
798  }
799  }
800 
801  return cellCount;
802 }
803 
804 
805 // ************************************************************************* //
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
Definition: pointFields.H:49
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
~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)
const double e
Elementary charge.
Definition: doubleFloat.H:78
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
autoPtr< polyMesh > createMesh(const fileName &name, labelTolabelPairHashTable &vertexMap, labelList &cellMap, const bool writeDelaunayData=true) const
Create an fvMesh from the triangulation.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
tetrahedron< point, const point & > tetPointRef
Definition: tetrahedron.H:78
static fvMesh * meshPtr
Definition: globalFoam.H:52
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)
Definition: UPstream.H:417
HashTable< label, labelPair, FixedList< label, 2 >::Hash<> > labelTolabelPairHashTable
Definition: DelaunayMesh.H:90
boundBox bounds() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
bool distribute(const boundBox &bb)
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
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)
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 barycentricCoords(const Foam::point &pt, scalarList &bary, Cell_handle &ch) const
Calculate and return the barycentric coordinates for.
void reset()
Clear the entire triangulation.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
label estimateCellCount(const autoPtr< backgroundMeshDecomposition > &decomposition) const
PointFrompoint toPoint(const Foam::point &p)
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:91
void distribute(const backgroundMeshDecomposition &decomposition)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static const char nl
Definition: Ostream.H:262
tensorField dumpAlignments() const
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:393
const dimensionedScalar c
Speed of light in a vacuum.
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)
A class for managing temporary objects.
Definition: PtrList.H:54
IOField< triad > triadIOField
triadField with IO.
Definition: triadIOField.H:42
Namespace for OpenFOAM.