All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cellSizeAndAlignmentGrid.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-2021 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 Application
25  Test-distributedDelaunayMesh
26 
27 Description
28 
29 \*---------------------------------------------------------------------------*/
30 
32 
33 #include "indexedVertex.H"
34 #include "indexedCell.H"
35 
36 #include "argList.H"
37 #include "Time.H"
40 #include "searchableSurfaces.H"
41 #include "conformationSurfaces.H"
42 #include "PrintTable.H"
43 #include "Random.H"
44 #include "boundBox.H"
45 #include "point.H"
46 #include "cellShapeControlMesh.H"
47 #include "triadField.H"
48 #include "scalarIOField.H"
49 #include "pointIOField.H"
50 #include "triadIOField.H"
51 
52 using namespace Foam;
53 
54 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
55 // Main program:
56 
57 template<class Triangulation, class Type>
58 Foam::tmp<Foam::Field<Type>> filterFarPoints
59 (
60  const Triangulation& mesh,
61  const Field<Type>& field
62 )
63 {
64  tmp<Field<Type>> tNewField(new Field<Type>(field.size()));
65  Field<Type>& newField = tNewField.ref();
66 
67  label added = 0;
68  label count = 0;
69 
70  for
71  (
72  typename Triangulation::Finite_vertices_iterator vit =
73  mesh.finite_vertices_begin();
74  vit != mesh.finite_vertices_end();
75  ++vit
76  )
77  {
78  if (vit->real())
79  {
80  newField[added++] = field[count];
81  }
82 
83  count++;
84  }
85 
86  newField.resize(added);
87 
88  return tNewField;
89 }
90 
91 
92 template<class T>
94 (
95  const T& mesh,
96  labelListList& pointPoints
97 )
98 {
99  pointPoints.setSize(mesh.vertexCount());
100 
101  globalIndex globalIndexing(mesh.vertexCount());
102 
103  for
104  (
105  typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
106  vit != mesh.finite_vertices_end();
107  ++vit
108  )
109  {
110  if (!vit->real())
111  {
112  continue;
113  }
114 
115  std::list<typename T::Vertex_handle> adjVerts;
116  mesh.finite_adjacent_vertices(vit, std::back_inserter(adjVerts));
117 
118  DynamicList<label> indices(adjVerts.size());
119 
120  for
121  (
122  typename std::list<typename T::Vertex_handle>::const_iterator
123  adjVertI = adjVerts.begin();
124  adjVertI != adjVerts.end();
125  ++adjVertI
126  )
127  {
128  typename T::Vertex_handle vh = *adjVertI;
129 
130  if (!vh->farPoint())
131  {
132  indices.append
133  (
134  globalIndexing.toGlobal(vh->procIndex(), vh->index())
135  );
136  }
137  }
138 
139  pointPoints[vit->index()].transfer(indices);
140  }
141 
142  List<Map<label>> compactMap;
143 
145  (
146  new mapDistribute
147  (
148  globalIndexing,
149  pointPoints,
150  compactMap
151  )
152  );
153 }
154 
155 
156 template<class T>
157 Foam::tmp<Foam::triadField> buildAlignmentField(const T& mesh)
158 {
159  tmp<triadField> tAlignments
160  (
161  new triadField(mesh.vertexCount(), triad::unset)
162  );
163  triadField& alignments = tAlignments.ref();
164 
165  for
166  (
167  typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
168  vit != mesh.finite_vertices_end();
169  ++vit
170  )
171  {
172  if (!vit->real())
173  {
174  continue;
175  }
176 
177  alignments[vit->index()] = vit->alignment();
178  }
179 
180  return tAlignments;
181 }
182 
183 
184 template<class T>
185 Foam::tmp<Foam::pointField> buildPointField(const T& mesh)
186 {
187  tmp<pointField> tPoints
188  (
189  new pointField(mesh.vertexCount(), point(great, great, great))
190  );
191  pointField& points = tPoints.ref();
192 
193  for
194  (
195  typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
196  vit != mesh.finite_vertices_end();
197  ++vit
198  )
199  {
200  if (!vit->real())
201  {
202  continue;
203  }
204 
205  points[vit->index()] = topoint(vit->point());
206  }
207 
208  return tPoints;
209 }
210 
211 
212 void refine
213 (
214  cellShapeControlMesh& mesh,
215  const conformationSurfaces& geometryToConformTo,
216  const label maxRefinementIterations,
217  const scalar defaultCellSize
218 )
219 {
220  for (label iter = 0; iter < maxRefinementIterations; ++iter)
221  {
222  DynamicList<point> ptsToInsert;
223 
224  for
225  (
226  CellSizeDelaunay::Finite_cells_iterator cit =
227  mesh.finite_cells_begin();
228  cit != mesh.finite_cells_end();
229  ++cit
230  )
231  {
232  const point newPoint =
233  topoint
234  (
235  CGAL::centroid
236  (
237  cit->vertex(0)->point(),
238  cit->vertex(1)->point(),
239  cit->vertex(2)->point(),
240  cit->vertex(3)->point()
241  )
242  );
243 
244  if (geometryToConformTo.inside(newPoint))
245  {
246  ptsToInsert.append(newPoint);
247  }
248  }
249 
250  Info<< " Adding " << returnReduce(ptsToInsert.size(), sumOp<label>())
251  << endl;
252 
253  forAll(ptsToInsert, ptI)
254  {
255  mesh.insert
256  (
257  ptsToInsert[ptI],
258  defaultCellSize,
259  triad::unset,
261  );
262  }
263  }
264 }
265 
266 
267 int main(int argc, char *argv[])
268 {
269  #include "setRootCase.H"
270  #include "createTime.H"
271 
272  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
273 
274  label maxRefinementIterations = 2;
275  label maxSmoothingIterations = 200;
276  scalar minResidual = 0;
277  scalar defaultCellSize = 0.001;
278  scalar nearFeatDistSqrCoeff = 1e-8;
279 
280 
281  // Need to decouple vertex and cell type from this class?
282  // Vertex must have:
283  // + index
284  // + procIndex
285  // - type should be optional
287 
288  IOdictionary foamyHexMeshDict
289  (
290  IOobject
291  (
292  "foamyHexMeshDict",
293  runTime.system(),
294  runTime,
297  )
298  );
299 
301 
302  searchableSurfaces allGeometry
303  (
304  IOobject
305  (
306  "cvSearchableSurfaces",
307  runTime.constant(),
309  runTime,
312  ),
313  foamyHexMeshDict.subDict("geometry"),
314  foamyHexMeshDict.lookupOrDefault("singleRegionName", true)
315  );
316 
317  conformationSurfaces geometryToConformTo
318  (
319  runTime,
320  rndGen,
321  allGeometry,
322  foamyHexMeshDict.subDict("surfaceConformation")
323  );
324 
326  if (Pstream::parRun())
327  {
328  bMesh.set
329  (
331  (
332  runTime,
333  rndGen,
334  geometryToConformTo,
335  foamyHexMeshDict.subDict("backgroundMeshDecomposition")
336  )
337  );
338  }
339 
340  // Nice to have IO for the delaunay mesh
341  // IO depend on vertex type.
342  //
343  // Define a delaunay mesh as:
344  // + list of points of the triangulation
345  // + optionally a list of cells
346 
347  Info<< nl << "Loop over surfaces" << endl;
348 
349  forAll(geometryToConformTo.surfaces(), sI)
350  {
351  const label surfI = geometryToConformTo.surfaces()[sI];
352 
353  const searchableSurface& surface =
354  geometryToConformTo.geometry()[surfI];
355 
356  Info<< nl << "Inserting points from surface " << surface.name()
357  << " (" << surface.type() << ")" << endl;
358 
359  const tmp<pointField> tpoints(surface.points());
360  const pointField& points = tpoints();
361 
362  Info<< " Number of points = " << points.size() << endl;
363 
364  forAll(points, pI)
365  {
366  // Is the point in the extendedFeatureEdgeMesh? If so get the
367  // point normal, otherwise get the surface normal from
368  // searchableSurface
369 
370  pointIndexHit info;
371  label infoFeature;
372  geometryToConformTo.findFeaturePointNearest
373  (
374  points[pI],
375  nearFeatDistSqrCoeff,
376  info,
377  infoFeature
378  );
379 
380 
381  autoPtr<triad> pointAlignment;
382 
383  if (info.hit())
384  {
385  const extendedFeatureEdgeMesh& features =
386  geometryToConformTo.features()[infoFeature];
387 
388  vectorField norms = features.featurePointNormals(info.index());
389 
390  // Create a triad from these norms.
391  pointAlignment.set(new triad());
392  forAll(norms, nI)
393  {
394  pointAlignment() += norms[nI];
395  }
396 
397  pointAlignment().normalise();
398  pointAlignment().orthogonalise();
399  }
400  else
401  {
402  geometryToConformTo.findEdgeNearest
403  (
404  points[pI],
405  nearFeatDistSqrCoeff,
406  info,
407  infoFeature
408  );
409 
410  if (info.hit())
411  {
412  const extendedFeatureEdgeMesh& features =
413  geometryToConformTo.features()[infoFeature];
414 
415  vectorField norms = features.edgeNormals(info.index());
416 
417  // Create a triad from these norms.
418  pointAlignment.set(new triad());
419  forAll(norms, nI)
420  {
421  pointAlignment() += norms[nI];
422  }
423 
424  pointAlignment().normalise();
425  pointAlignment().orthogonalise();
426  }
427  else
428  {
429  pointField ptField(1, points[pI]);
430  scalarField distField(1, nearFeatDistSqrCoeff);
431  List<pointIndexHit> infoList(1, pointIndexHit());
432 
433  surface.findNearest(ptField, distField, infoList);
434 
435  vectorField normals(1);
436  surface.getNormal(infoList, normals);
437 
438  pointAlignment.set(new triad(normals[0]));
439  }
440  }
441 
442  if (Pstream::parRun())
443  {
444  if (bMesh().positionOnThisProcessor(points[pI]))
445  {
446  CellSizeDelaunay::Vertex_handle vh = mesh.insert
447  (
448  points[pI],
449  defaultCellSize,
450  pointAlignment(),
452  );
453  }
454  }
455  else
456  {
457  CellSizeDelaunay::Vertex_handle vh = mesh.insert
458  (
459  points[pI],
460  defaultCellSize,
461  pointAlignment(),
463  );
464  }
465  }
466  }
467 
468 
469  // Refine the mesh
470  refine
471  (
472  mesh,
473  geometryToConformTo,
474  maxRefinementIterations,
475  defaultCellSize
476  );
477 
478 
479  if (Pstream::parRun())
480  {
481  mesh.distribute(bMesh);
482  }
483 
484 
485  labelListList pointPoints;
486  autoPtr<mapDistribute> meshDistributor = buildMap(mesh, pointPoints);
487 
488 
489  triadField alignments(buildAlignmentField(mesh));
490  pointField points(buildPointField(mesh));
491 
492  mesh.printInfo(Info);
493 
494 
495  // Setup the sizes and alignments on each point
496  triadField fixedAlignments(mesh.vertexCount(), triad::unset);
497 
498  for
499  (
500  CellSizeDelaunay::Finite_vertices_iterator vit =
501  mesh.finite_vertices_begin();
502  vit != mesh.finite_vertices_end();
503  ++vit
504  )
505  {
506  if (vit->nearBoundary())
507  {
508  fixedAlignments[vit->index()] = vit->alignment();
509  }
510  }
511 
512  Info<< nl << "Smoothing alignments" << endl;
513 
514  for (label iter = 0; iter < maxSmoothingIterations; iter++)
515  {
516  Info<< "Iteration " << iter;
517 
518  meshDistributor().distribute(points);
519  meshDistributor().distribute(alignments);
520 
521  scalar residual = 0;
522 
523  triadField triadAv(alignments.size(), triad::unset);
524 
525  forAll(pointPoints, pI)
526  {
527  const labelList& pPoints = pointPoints[pI];
528 
529  if (pPoints.empty())
530  {
531  continue;
532  }
533 
534  const triad& oldTriad = alignments[pI];
535  triad& newTriad = triadAv[pI];
536 
537  // Enforce the boundary conditions
538  const triad& fixedAlignment = fixedAlignments[pI];
539 
540  forAll(pPoints, adjPointi)
541  {
542  const label adjPointIndex = pPoints[adjPointi];
543 
544  scalar dist = mag(points[pI] - points[adjPointIndex]);
545 
546 // dist = max(dist, small);
547 
548  triad tmpTriad = alignments[adjPointIndex];
549 
550  for (direction dir = 0; dir < 3; dir++)
551  {
552  if (tmpTriad.set(dir))
553  {
554  tmpTriad[dir] *= (1.0/dist);
555  }
556  }
557 
558  newTriad += tmpTriad;
559  }
560 
561  newTriad.normalise();
562  newTriad.orthogonalise();
563 // newTriad = newTriad.sortxyz();
564 
565  label nFixed = 0;
566 
567  forAll(fixedAlignment, dirI)
568  {
569  if (fixedAlignment.set(dirI))
570  {
571  nFixed++;
572  }
573  }
574 
575  if (nFixed == 1)
576  {
577  forAll(fixedAlignment, dirI)
578  {
579  if (fixedAlignment.set(dirI))
580  {
581  newTriad.align(fixedAlignment[dirI]);
582  }
583  }
584  }
585  else if (nFixed == 2)
586  {
587  forAll(fixedAlignment, dirI)
588  {
589  if (fixedAlignment.set(dirI))
590  {
591  newTriad[dirI] = fixedAlignment[dirI];
592  }
593  else
594  {
595  newTriad[dirI] = triad::unset[dirI];
596  }
597  }
598 
599  newTriad.orthogonalise();
600  }
601  else if (nFixed == 3)
602  {
603  forAll(fixedAlignment, dirI)
604  {
605  if (fixedAlignment.set(dirI))
606  {
607  newTriad[dirI] = fixedAlignment[dirI];
608  }
609  }
610  }
611 
612  for (direction dir = 0; dir < 3; ++dir)
613  {
614  if
615  (
616  newTriad.set(dir)
617  && oldTriad.set(dir)
618  //&& !fixedAlignment.set(dir)
619  )
620  {
621  scalar dotProd = (oldTriad[dir] & newTriad[dir]);
622  scalar diff = mag(dotProd) - 1.0;
623 
624  residual += mag(diff);
625  }
626  }
627  }
628 
629  forAll(alignments, pI)
630  {
631  alignments[pI] = triadAv[pI].sortxyz();
632  }
633 
634  reduce(residual, sumOp<scalar>());
635 
636  Info<< ", Residual = " << residual << endl;
637 
638  if (residual <= minResidual)
639  {
640  break;
641  }
642  }
643 
644 
645  // Write alignments to a .obj file
646  OFstream str(runTime.path()/"alignments.obj");
647 
648  forAll(alignments, pI)
649  {
650  const triad& tri = alignments[pI];
651 
652  if (tri.set())
653  {
654  forAll(tri, dirI)
655  {
656  meshTools::writeOBJ(str, points[pI], tri[dirI] + points[pI]);
657  }
658  }
659  }
660 
661 
662  // Remove the far points
663  pointIOField pointsIO
664  (
665  IOobject
666  (
667  "points",
668  runTime.constant(),
669  runTime,
672  ),
673  filterFarPoints(mesh, points)
674  );
675 
676  scalarField sizes(points.size(), defaultCellSize);
677  scalarIOField sizesIO
678  (
679  IOobject
680  (
681  "sizes",
682  runTime.constant(),
683  runTime,
686  ),
687  filterFarPoints(mesh, sizes)
688  );
689 
690  triadIOField alignmentsIO
691  (
692  IOobject
693  (
694  "alignments",
695  runTime.constant(),
696  runTime,
699  ),
700  filterFarPoints(mesh, alignments)
701  );
702 
703  pointsIO.write();
704  sizesIO.write();
705  alignmentsIO.write();
706 
707  Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
708  << " ClockTime = " << runTime.elapsedClockTime() << " s"
709  << nl << endl;
710 
711  Info<< "\nEnd\n" << endl;
712 
713  return 0;
714 }
715 
716 
717 // ************************************************************************* //
Field< bool > inside(const pointField &samplePts) const
Check if points are inside surfaces to conform to.
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.
Definition: triad.C:403
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:325
const searchableSurfaces & geometry() const
Return reference to the searchableSurfaces object containing all.
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
const word & name() const
Return name.
Definition: IOobject.H:303
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
pointFromPoint topoint(const Point &P)
const labelListList & edgeNormals() const
Return the indices of the normals that are adjacent to the.
void align(const vector &v)
Align this triad with the given vector v.
Definition: triad.C:255
const labelList & surfaces() const
Return the surface indices.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
Output to file stream.
Definition: OFstream.H:82
uint8_t direction
Definition: direction.H:45
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label count(const ListType &l, typename ListType::const_reference x)
Count the number of occurrences of a value in a list.
void findFeaturePointNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &fpHit, label &featureHit) const
Find the nearest point on any feature edge.
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...
Definition: PointIndexHit.H:53
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
void resize(const label)
Alias for setSize(const label)
Definition: ListI.H:138
bool set(const direction d) const
Is the vector in the direction d set.
Definition: triadI.H:61
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
void orthogonalise()
Orthogonalise this triad so that it is ortho-normal.
Definition: triad.C:108
Random rndGen(label(0))
void printInfo(Ostream &os) const
Write mesh statistics to stream.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
const PtrList< extendedFeatureEdgeMesh > & features() const
Return the object holding the feature points and edges.
dynamicFvMesh & mesh
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const =0
const pointField & points
static const triad unset
Definition: triad.H:99
Pre-declare SubField and related Field type.
Definition: Field.H:56
static const word & geometryDir()
Return the geometry directory name.
const word & constant() const
Return constant name.
Definition: TimePaths.H:123
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
bool hit() const
Is there a hit.
Container for searchableSurfaces.
void findEdgeNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &edgeHit, label &featureHit) const
Find the nearest point on any feature edge.
void distribute(const backgroundMeshDecomposition &decomposition)
Random number generator.
Definition: Random.H:57
const word & system() const
Return system name.
Definition: TimePaths.H:113
Field< triad > triadField
Specialisation of Field<T> for triad.
Definition: triadField.H:49
virtual tmp< pointField > points() const =0
Get the points that define the surface.
static const char nl
Definition: Ostream.H:260
void set(T *)
Set pointer to that given.
Definition: autoPtrI.H:99
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.
Definition: triad.H:65
void setSize(const label)
Reset size of List.
Definition: List.C:281
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
time_t elapsedClockTime() const
Returns wall-clock time from clock instantiation.
Definition: clock.C:122
Class containing processor-to-processor mapping information.
label vertexCount() const
Return the vertex count (the next unique vertex index)
double elapsedCpuTime() const
Return CPU time (in seconds) from the start.
Definition: cpuTime.C:67
vector point
Point is a vector.
Definition: point.H:41
void normalise()
Normalise each set axis vector to have a unit magnitude.
Definition: triadI.H:111
const polyBoundaryMesh & bMesh
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)
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...
Definition: PtrList.H:52
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
virtual bool write(const bool write=true) const
Write using setting from DB.
A class for managing temporary objects.
Definition: PtrList.H:53
fileName path() const
Return path.
Definition: TimePaths.H:139
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
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.
Definition: IOField.H:50
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
label index() const
Return index.
Namespace for OpenFOAM.