CV2D.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) 2013-2018 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 "CV2D.H"
27 #include "Random.H"
28 #include "transform.H"
29 #include "IFstream.H"
30 #include "uint.H"
31 
32 namespace Foam
33 {
34  defineTypeNameAndDebug(CV2D, 0);
35 }
36 
37 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
38 
39 void Foam::CV2D::insertBoundingBox()
40 {
41  Info<< "insertBoundingBox: creating bounding mesh" << endl;
42  scalar bigSpan = 10*meshControls().span();
43  insertPoint(point2D(-bigSpan, -bigSpan), Vb::FAR_POINT);
44  insertPoint(point2D(-bigSpan, bigSpan), Vb::FAR_POINT);
45  insertPoint(point2D(bigSpan, -bigSpan), Vb::FAR_POINT);
46  insertPoint(point2D(bigSpan, bigSpan), Vb::FAR_POINT);
47 }
48 
49 
50 void Foam::CV2D::fast_restore_Delaunay(Vertex_handle vh)
51 {
52  int i;
53  Face_handle f = vh->face(), next, start(f);
54 
55  do
56  {
57  i=f->index(vh);
58  if (!is_infinite(f))
59  {
60  if (!internal_flip(f, cw(i))) external_flip(f, i);
61  if (f->neighbor(i) == start) start = f;
62  }
63  f = f->neighbor(cw(i));
64  } while (f != start);
65 }
66 
67 
68 void Foam::CV2D::external_flip(Face_handle& f, int i)
69 {
70  Face_handle n = f->neighbor(i);
71 
72  if
73  (
74  CGAL::ON_POSITIVE_SIDE
75  != side_of_oriented_circle(n, f->vertex(i)->point())
76  ) return;
77 
78  flip(f, i);
79  i = n->index(f->vertex(i));
80  external_flip(n, i);
81 }
82 
83 
84 bool Foam::CV2D::internal_flip(Face_handle& f, int i)
85 {
86  Face_handle n = f->neighbor(i);
87 
88  if
89  (
90  CGAL::ON_POSITIVE_SIDE
91  != side_of_oriented_circle(n, f->vertex(i)->point())
92  )
93  {
94  return false;
95  }
96 
97  flip(f, i);
98 
99  return true;
100 }
101 
102 
103 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
104 
105 Foam::CV2D::CV2D
106 (
107  const Time& runTime,
108  const dictionary& cvMeshDict
109 )
110 :
111  Delaunay(),
112  runTime_(runTime),
113  rndGen_(64293*Pstream::myProcNo()),
114  allGeometry_
115  (
116  IOobject
117  (
118  "cvSearchableSurfaces",
119  runTime_.constant(),
120  "triSurface",
121  runTime_,
122  IOobject::MUST_READ,
123  IOobject::NO_WRITE
124  ),
125  cvMeshDict.subDict("geometry"),
126  cvMeshDict.lookupOrDefault("singleRegionName", true)
127  ),
128  qSurf_
129  (
130  runTime_,
131  rndGen_,
132  allGeometry_,
133  cvMeshDict.subDict("surfaceConformation")
134  ),
135  controls_(cvMeshDict, qSurf_.globalBounds()),
136  cellSizeControl_
137  (
138  runTime,
139  cvMeshDict.subDict("motionControl").subDict("shapeControlFunctions"),
140  qSurf_,
141  controls_.minCellSize()
142  ),
143  relaxationModel_
144  (
145  relaxationModel::New
146  (
147  cvMeshDict.subDict("motionControl"),
148  runTime
149  )
150  ),
151  z_
152  (
153  point
154  (
155  cvMeshDict.subDict("surfaceConformation").lookup("locationInMesh")
156  ).z()
157  ),
158  startOfInternalPoints_(0),
159  startOfSurfacePointPairs_(0),
160  startOfBoundaryConformPointPairs_(0),
161  featurePoints_()
162 {
163  Info<< meshControls() << endl;
164 
165  insertBoundingBox();
166  insertFeaturePoints();
167 }
168 
169 
170 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
171 
173 {}
174 
175 
176 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
177 
179 (
180  const point2DField& points,
181  const scalar nearness
182 )
183 {
184  Info<< "insertInitialPoints(const point2DField& points): ";
185 
186  startOfInternalPoints_ = number_of_vertices();
187  label nVert = startOfInternalPoints_;
188 
189  // Add the points and index them
190  forAll(points, i)
191  {
192  const point2D& p = points[i];
193 
194  if (qSurf_.wellInside(toPoint3D(p), nearness))
195  {
196  insert(toPoint(p))->index() = nVert++;
197  }
198  else
199  {
200  Warning
201  << "Rejecting point " << p << " outside surface" << endl;
202  }
203  }
204 
205  Info<< nVert << " vertices inserted" << endl;
206 
207  if (meshControls().objOutput())
208  {
209  // Checking validity of triangulation
210  assert(is_valid());
211 
212  writeTriangles("initial_triangles.obj", true);
213  writeFaces("initial_faces.obj", true);
214  }
215 }
216 
217 
218 void Foam::CV2D::insertPoints(const fileName& pointFileName)
219 {
220  IFstream pointsFile(pointFileName);
221 
222  if (pointsFile.good())
223  {
225  (
226  point2DField(pointsFile),
227  0.5*meshControls().minCellSize2()
228  );
229  }
230  else
231  {
233  << "Could not open pointsFile " << pointFileName
234  << exit(FatalError);
235  }
236 }
237 
238 
240 {
241  Info<< "insertInitialGrid: ";
242 
243  startOfInternalPoints_ = number_of_vertices();
244  label nVert = startOfInternalPoints_;
245 
246  scalar x0 = qSurf_.globalBounds().min().x();
247  scalar xR = qSurf_.globalBounds().max().x() - x0;
248  int ni = int(xR/meshControls().minCellSize()) + 1;
249  scalar deltax = xR/ni;
250 
251  scalar y0 = qSurf_.globalBounds().min().y();
252  scalar yR = qSurf_.globalBounds().max().y() - y0;
253  int nj = int(yR/meshControls().minCellSize()) + 1;
254  scalar deltay = yR/nj;
255 
256  Random rndGen(1321);
257  scalar pert = meshControls().randomPerturbation()*min(deltax, deltay);
258 
259  for (int i=0; i<ni; i++)
260  {
261  for (int j=0; j<nj; j++)
262  {
263  point p(x0 + i*deltax, y0 + j*deltay, 0);
264 
265  if (meshControls().randomiseInitialGrid())
266  {
267  p.x() += pert*(rndGen.scalar01() - 0.5);
268  p.y() += pert*(rndGen.scalar01() - 0.5);
269  }
270 
271  if (qSurf_.wellInside(p, 0.5*meshControls().minCellSize2()))
272  {
273  insert(Point(p.x(), p.y()))->index() = nVert++;
274  }
275  }
276  }
277 
278  Info<< nVert << " vertices inserted" << endl;
279 
280  if (meshControls().objOutput())
281  {
282  // Checking validity of triangulation
283  assert(is_valid());
284 
285  writeTriangles("initial_triangles.obj", true);
286  writeFaces("initial_faces.obj", true);
287  }
288 }
289 
290 
292 {
293  startOfSurfacePointPairs_ = number_of_vertices();
294 
295  if (meshControls().insertSurfaceNearestPointPairs())
296  {
297  insertSurfaceNearestPointPairs();
298  }
299 
300  write("nearest");
301 
302  // Insertion of point-pairs for near-points may cause protrusions
303  // so insertBoundaryConformPointPairs must be executed last
304  if (meshControls().insertSurfaceNearPointPairs())
305  {
306  insertSurfaceNearPointPairs();
307  }
308 
309  startOfBoundaryConformPointPairs_ = number_of_vertices();
310 }
311 
312 
314 {
315  if (!meshControls().insertSurfaceNearestPointPairs())
316  {
317  markNearBoundaryPoints();
318  }
319 
320  // Mark all the faces as SAVE_CHANGED
321  for
322  (
323  Triangulation::Finite_faces_iterator fit = finite_faces_begin();
324  fit != finite_faces_end();
325  fit++
326  )
327  {
328  fit->faceIndex() = Fb::SAVE_CHANGED;
329  }
330 
331  for (label iter=1; iter<=meshControls().maxBoundaryConformingIter(); iter++)
332  {
333  label nIntersections = insertBoundaryConformPointPairs
334  (
335  "surfaceIntersections_" + Foam::name(iter) + ".obj"
336  );
337 
338  if (nIntersections == 0)
339  {
340  break;
341  }
342  else
343  {
344  Info<< "BC iteration " << iter << ": "
345  << nIntersections << " point-pairs inserted" << endl;
346  }
347 
348  // Any faces changed by insertBoundaryConformPointPairs will now
349  // be marked CHANGED, mark those as SAVE_CHANGED and those that
350  // remained SAVE_CHANGED as UNCHANGED
351  for
352  (
353  Triangulation::Finite_faces_iterator fit = finite_faces_begin();
354  fit != finite_faces_end();
355  fit++
356  )
357  {
358  if (fit->faceIndex() == Fb::SAVE_CHANGED)
359  {
360  fit->faceIndex() = Fb::UNCHANGED;
361  }
362  else if (fit->faceIndex() == Fb::CHANGED)
363  {
364  fit->faceIndex() = Fb::SAVE_CHANGED;
365  }
366  }
367  }
368 
369  Info<< nl;
370 
371  write("boundary");
372 }
373 
374 
376 {
377  for
378  (
379  Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
380  vit != finite_vertices_end();
381  ++vit
382  )
383  {
384  if (vit->index() >= startOfSurfacePointPairs_)
385  {
386  remove(vit);
387  }
388  }
389 }
390 
391 
393 {
394  const scalar relaxation = relaxationModel_->relaxation();
395 
396  Info<< "Relaxation = " << relaxation << endl;
397 
398  Field<point2D> dualVertices(number_of_faces());
399 
400  label dualVerti = 0;
401 
402  // Find the dual point of each tetrahedron and assign it an index.
403  for
404  (
405  Triangulation::Finite_faces_iterator fit = finite_faces_begin();
406  fit != finite_faces_end();
407  ++fit
408  )
409  {
410  fit->faceIndex() = -1;
411 
412  if
413  (
414  fit->vertex(0)->internalOrBoundaryPoint()
415  || fit->vertex(1)->internalOrBoundaryPoint()
416  || fit->vertex(2)->internalOrBoundaryPoint()
417  )
418  {
419  fit->faceIndex() = dualVerti;
420 
421  dualVertices[dualVerti] = toPoint2D(circumcenter(fit));
422 
423  dualVerti++;
424  }
425  }
426 
427  dualVertices.setSize(dualVerti);
428 
429  Field<vector2D> displacementAccumulator
430  (
431  startOfSurfacePointPairs_,
433  );
434 
435  // Calculate target size and alignment for vertices
436  scalarField sizes
437  (
438  number_of_vertices(),
439  meshControls().minCellSize()
440  );
441 
442  Field<vector2D> alignments
443  (
444  number_of_vertices(),
445  vector2D(1, 0)
446  );
447 
448  for
449  (
450  Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
451  vit != finite_vertices_end();
452  ++vit
453  )
454  {
455  if (vit->internalOrBoundaryPoint())
456  {
457  point2D vert = toPoint2D(vit->point());
458 
459  // alignment and size determination
460  pointIndexHit pHit;
461  label hitSurface = -1;
462 
463  qSurf_.findSurfaceNearest
464  (
465  toPoint3D(vert),
466  meshControls().span2(),
467  pHit,
468  hitSurface
469  );
470 
471  if (pHit.hit())
472  {
473  vectorField norm(1);
474  allGeometry_[hitSurface].getNormal
475  (
476  List<pointIndexHit>(1, pHit),
477  norm
478  );
479 
480  alignments[vit->index()] = toPoint2D(norm[0]);
481 
482  sizes[vit->index()] =
483  cellSizeControl_.cellSize
484  (
485  toPoint3D(vit->point())
486  );
487  }
488  }
489  }
490 
491  // Info<< "Calculated alignments" << endl;
492 
493  scalar cosAlignmentAcceptanceAngle = 0.68;
494 
495  // Upper and lower edge length ratios for weight
496  scalar u = 1.0;
497  scalar l = 0.7;
498 
499  PackedBoolList pointToBeRetained(startOfSurfacePointPairs_, true);
500 
501  std::list<Point> pointsToInsert;
502 
503  for
504  (
505  Triangulation::Finite_edges_iterator eit = finite_edges_begin();
506  eit != finite_edges_end();
507  eit++
508  )
509  {
510  Vertex_handle vA = eit->first->vertex(cw(eit->second));
511  Vertex_handle vB = eit->first->vertex(ccw(eit->second));
512 
513  if (!vA->internalOrBoundaryPoint() || !vB->internalOrBoundaryPoint())
514  {
515  continue;
516  }
517 
518  const point2D& dualV1 = dualVertices[eit->first->faceIndex()];
519  const point2D& dualV2 =
520  dualVertices[eit->first->neighbor(eit->second)->faceIndex()];
521 
522  scalar dualEdgeLength = mag(dualV1 - dualV2);
523 
524  point2D dVA = toPoint2D(vA->point());
525  point2D dVB = toPoint2D(vB->point());
526 
527  Field<vector2D> alignmentDirsA(2);
528 
529  alignmentDirsA[0] = alignments[vA->index()];
530  alignmentDirsA[1] = vector2D
531  (
532  -alignmentDirsA[0].y(),
533  alignmentDirsA[0].x()
534  );
535 
536  Field<vector2D> alignmentDirsB(2);
537 
538  alignmentDirsB[0] = alignments[vB->index()];
539  alignmentDirsB[1] = vector2D
540  (
541  -alignmentDirsB[0].y(),
542  alignmentDirsB[0].x()
543  );
544 
545  Field<vector2D> alignmentDirs(alignmentDirsA);
546 
547  forAll(alignmentDirsA, aA)
548  {
549  const vector2D& a(alignmentDirsA[aA]);
550 
551  scalar maxDotProduct = 0.0;
552 
553  forAll(alignmentDirsB, aB)
554  {
555  const vector2D& b(alignmentDirsB[aB]);
556 
557  scalar dotProduct = a & b;
558 
559  if (mag(dotProduct) > maxDotProduct)
560  {
561  maxDotProduct = mag(dotProduct);
562 
563  alignmentDirs[aA] = a + sign(dotProduct)*b;
564 
565  alignmentDirs[aA] /= mag(alignmentDirs[aA]);
566  }
567  }
568  }
569 
570  vector2D rAB = dVA - dVB;
571 
572  scalar rABMag = mag(rAB);
573 
574  forAll(alignmentDirs, aD)
575  {
576  vector2D& alignmentDir = alignmentDirs[aD];
577 
578  if ((rAB & alignmentDir) < 0)
579  {
580  // swap the direction of the alignment so that has the
581  // same sense as rAB
582  alignmentDir *= -1;
583  }
584 
585  scalar alignmentDotProd = ((rAB/rABMag) & alignmentDir);
586 
587  if (alignmentDotProd > cosAlignmentAcceptanceAngle)
588  {
589  scalar targetFaceSize =
590  0.5*(sizes[vA->index()] + sizes[vB->index()]);
591 
592  // Test for changing aspect ratio on second alignment (first
593  // alignment is neartest surface normal)
594  // if (aD == 1)
595  // {
596  // targetFaceSize *= 2.0;
597  // }
598 
599  alignmentDir *= 0.5*targetFaceSize;
600 
601  vector2D delta = alignmentDir - 0.5*rAB;
602 
603  if (dualEdgeLength < 0.7*targetFaceSize)
604  {
605  delta *= 0;
606  }
607  else if (dualEdgeLength < targetFaceSize)
608  {
609  delta *=
610  (
611  dualEdgeLength
612  /(targetFaceSize*(u - l))
613  - 1/((u/l) - 1)
614  );
615  }
616 
617  if
618  (
619  vA->internalPoint()
620  && vB->internalPoint()
621  && rABMag > 1.75*targetFaceSize
622  && dualEdgeLength > 0.05*targetFaceSize
623  && alignmentDotProd > 0.93
624  )
625  {
626  // Point insertion
627  pointsToInsert.push_back(toPoint(0.5*(dVA + dVB)));
628  }
629  else if
630  (
631  (vA->internalPoint() || vB->internalPoint())
632  && rABMag < 0.65*targetFaceSize
633  )
634  {
635  // Point removal
636 
637  // Only insert a point at the midpoint of the short edge
638  // if neither attached point has already been identified
639  // to be removed.
640  if
641  (
642  pointToBeRetained[vA->index()] == true
643  && pointToBeRetained[vB->index()] == true
644  )
645  {
646  pointsToInsert.push_back(toPoint(0.5*(dVA + dVB)));
647  }
648 
649  if (vA->internalPoint())
650  {
651  pointToBeRetained[vA->index()] = false;
652  }
653 
654  if (vB->internalPoint())
655  {
656  pointToBeRetained[vB->index()] = false;
657  }
658  }
659  else
660  {
661  if (vA->internalPoint())
662  {
663  displacementAccumulator[vA->index()] += delta;
664  }
665 
666  if (vB->internalPoint())
667  {
668  displacementAccumulator[vB->index()] += -delta;
669  }
670  }
671  }
672  }
673  }
674 
675  vector2D totalDisp = sum(displacementAccumulator);
676  scalar totalDist = sum(mag(displacementAccumulator));
677 
678  // Relax the calculated displacement
679  displacementAccumulator *= relaxation;
680 
681  label numberOfNewPoints = pointsToInsert.size();
682 
683  for
684  (
685  Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
686  vit != finite_vertices_end();
687  ++vit
688  )
689  {
690  if (vit->internalPoint())
691  {
692  if (pointToBeRetained[vit->index()])
693  {
694  pointsToInsert.push_front
695  (
696  toPoint
697  (
698  toPoint2D(vit->point())
699  + displacementAccumulator[vit->index()]
700  )
701  );
702  }
703  }
704  }
705 
706  // Clear the triangulation and reinsert the bounding box and feature points.
707  // This is faster than removing and moving points.
708  this->clear();
709 
710  insertBoundingBox();
711 
712  reinsertFeaturePoints();
713 
714  startOfInternalPoints_ = number_of_vertices();
715 
716  label nVert = startOfInternalPoints_;
717 
718  Info<< "Inserting " << numberOfNewPoints << " new points" << endl;
719 
720  // Use the range insert as it is faster than individually inserting points.
721  insert(pointsToInsert.begin(), pointsToInsert.end());
722 
723  for
724  (
725  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
726  vit != finite_vertices_end();
727  ++vit
728  )
729  {
730  if
731  (
732  vit->type() == Vb::INTERNAL_POINT
733  && vit->index() == Vb::INTERNAL_POINT
734  )
735  {
736  vit->index() = nVert++;
737  }
738  }
739 
740  Info<< " Total displacement = " << totalDisp << nl
741  << " Total distance = " << totalDist << nl
742  << " Points added = " << pointsToInsert.size()
743  << endl;
744 
745  write("internal");
746 
748 
749  boundaryConform();
750 
751 
752  // Old Method
753  /*
754  for
755  (
756  Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
757  vit != finite_vertices_end();
758  ++vit
759  )
760  {
761  if (vit->internalPoint())
762  {
763  // Current dual-cell defining vertex ("centre")
764  point2DFromPoint defVert0 = toPoint2D(vit->point());
765 
766  Triangulation::Edge_circulator ec = incident_edges(vit);
767  Triangulation::Edge_circulator ecStart = ec;
768 
769  // Circulate around the edges to find the first which is not
770  // infinite
771  do
772  {
773  if (!is_infinite(ec)) break;
774  } while (++ec != ecStart);
775 
776  // Store the start-end of the first non-infinte edge
777  point2D de0 = toPoint2D(circumcenter(ec->first));
778 
779  // Keep track of the maximum edge length^2
780  scalar maxEdgeLen2 = 0.0;
781 
782  // Keep track of the index of the longest edge
783  label edgecd0i = -1;
784 
785  // Edge counter
786  label edgei = 0;
787 
788  do
789  {
790  if (!is_infinite(ec))
791  {
792  // Get the end of the current edge
793  point2D de1 = toPoint2D
794  (
795  circumcenter(ec->first->neighbor(ec->second))
796  );
797 
798  // Store the current edge vector
799  edges[edgei] = de1 - de0;
800 
801  // Store the edge mid-point in the vertices array
802  vertices[edgei] = 0.5*(de1 + de0);
803 
804  // Move the current edge end into the edge start for the
805  // next iteration
806  de0 = de1;
807 
808  // Keep track of the longest edge
809 
810  scalar edgeLen2 = magSqr(edges[edgei]);
811 
812  if (edgeLen2 > maxEdgeLen2)
813  {
814  maxEdgeLen2 = edgeLen2;
815  edgecd0i = edgei;
816  }
817 
818  edgei++;
819  }
820  } while (++ec != ecStart);
821 
822  // Initialise cd0 such that the mesh will align
823  // in in the x-y directions
824  vector2D cd0(1, 0);
825 
826  if (meshControls().relaxOrientation())
827  {
828  // Get the longest edge from the array and use as the primary
829  // direction of the coordinate system of the "square" cell
830  cd0 = edges[edgecd0i];
831  }
832 
833  if (meshControls().nearWallAlignedDist() > 0)
834  {
835  pointIndexHit pHit = qSurf_.tree().findNearest
836  (
837  toPoint3D(defVert0),
838  meshControls().nearWallAlignedDist2()
839  );
840 
841  if (pHit.hit())
842  {
843  cd0 = toPoint2D(faceNormals[pHit.index()]);
844  }
845  }
846 
847  // Rotate by 45deg needed to create an averaging procedure which
848  // encourages the cells to be square
849  cd0 = vector2D(cd0.x() + cd0.y(), cd0.y() - cd0.x());
850 
851  // Normalise the primary coordinate direction
852  cd0 /= mag(cd0);
853 
854  // Calculate the orthogonal coordinate direction
855  vector2D cd1(-cd0.y(), cd0.x());
856 
857 
858  // Restart the circulator
859  ec = ecStart;
860 
861  // ... and the counter
862  edgei = 0;
863 
864  // Initialise the displacement for the centre and sum-weights
865  vector2D disp = Zero;
866  scalar sumw = 0;
867 
868  do
869  {
870  if (!is_infinite(ec))
871  {
872  // Pick up the current edge
873  const vector2D& ei = edges[edgei];
874 
875  // Calculate the centre to edge-centre vector
876  vector2D deltai = vertices[edgei] - defVert0;
877 
878  // Set the weight for this edge contribution
879  scalar w = 1;
880 
881  if (meshControls().squares())
882  {
883  w = magSqr(deltai.x()*ei.y() - deltai.y()*ei.x());
884  // alternative weights
885  // w = mag(deltai.x()*ei.y() - deltai.y()*ei.x());
886  // w = magSqr(ei)*mag(deltai);
887 
888  // Use the following for an ~square mesh
889  // Find the coordinate contributions for this edge delta
890  scalar cd0deltai = cd0 & deltai;
891  scalar cd1deltai = cd1 & deltai;
892 
893  // Create a "square" displacement
894  if (mag(cd0deltai) > mag(cd1deltai))
895  {
896  disp += (w*cd0deltai)*cd0;
897  }
898  else
899  {
900  disp += (w*cd1deltai)*cd1;
901  }
902  }
903  else
904  {
905  // Use this for a hexagon/pentagon mesh
906  disp += w*deltai;
907  }
908 
909  // Sum the weights
910  sumw += w;
911  }
912  else
913  {
914  FatalErrorInFunction
915  << "Infinite triangle found in internal mesh"
916  << exit(FatalError);
917  }
918 
919  edgei++;
920 
921  } while (++ec != ecStart);
922 
923  // Calculate the average displacement
924  disp /= sumw;
925  totalDisp += disp;
926  totalDist += mag(disp);
927 
928  // Move the point by a fraction of the average displacement
929  movePoint(vit, defVert0 + relaxation*disp);
930  }
931  }
932 
933  Info << "\nTotal displacement = " << totalDisp
934  << " total distance = " << totalDist << endl;
935  */
936 }
937 
938 /*
939 void Foam::CV2D::moveInternalPoints(const point2DField& newPoints)
940 {
941  label pointi = 0;
942 
943  for
944  (
945  Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
946  vit != finite_vertices_end();
947  ++vit
948  )
949  {
950  if (vit->internalPoint())
951  {
952  movePoint(vit, newPoints[pointi++]);
953  }
954  }
955 }
956 */
957 
958 void Foam::CV2D::write() const
959 {
960  if (meshControls().objOutput())
961  {
962  writeFaces("allFaces.obj", false);
963  writeFaces("faces.obj", true);
964  writeTriangles("allTriangles.obj", false);
965  writeTriangles("triangles.obj", true);
966  writePatch("patch.pch");
967  }
968 }
969 
970 
971 void Foam::CV2D::write(const word& stage) const
972 {
973  if (meshControls().objOutput())
974  {
975  Foam::mkDir(stage + "Faces");
976  Foam::mkDir(stage + "Triangles");
977 
978  writeFaces
979  (
980  stage
981  + "Faces/allFaces_"
982  + runTime_.timeName()
983  + ".obj",
984  false
985  );
986 
988  (
989  stage
990  + "Triangles/allTriangles_"
991  + runTime_.timeName()
992  + ".obj",
993  false
994  );
995  }
996 }
997 
998 
999 // ************************************************************************* //
scalar delta
dimensionedScalar sign(const dimensionedScalar &ds)
CGAL::Delaunay_triangulation_3< K, Tds, CompactLocator > Delaunay
scalar cellSize(const point &pt) const
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
vector2D point2D
Point2D is a vector.
Definition: point2D.H:41
void writePatch(const fileName &fName) const
Write patch.
const treeBoundBox & globalBounds() const
Return the global bounds.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Field< bool > wellInside(const pointField &samplePts, const scalarField &testDistSqr) const
Check if point is inside surfaces to conform to by at least.
~CV2D()
Destructor.
tUEqn clear()
void insertSurfacePointPairs()
Insert all surface point-pairs from.
Vector2D< scalar > vector2D
vector2D obtained from generic Vector2D
Definition: vector2D.H:49
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
point toPoint3D(const point2D &) const
Definition: CV2DI.H:141
const cv2DControls & meshControls() const
Definition: CV2DI.H:118
void newPoints()
Move the internal points to the given new locations and update.
dimensionedScalar y0(const dimensionedScalar &ds)
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
const Cmpt & y() const
Definition: VectorI.H:81
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
Random rndGen(label(0))
const point2D & toPoint2D(const point &) const
Definition: CV2DI.H:124
scalar y
void insertGrid()
Create the initial mesh as a regular grid of points.
void writeFaces(const fileName &fName, bool internalOnly) const
Write dual faces as .obj file.
stressControl lookup("compactNormalStress") >> compactNormalStress
PointFromPoint2D toPoint(const point2D &) const
Definition: CV2DI.H:168
System uinteger.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
3D tensor transformation operations.
timeIndices insert(timeIndex, timeDirs[timeI].value())
void findSurfaceNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &surfHit, label &hitSurface) const
Find the nearest point to the sample and return it to the.
const Cmpt & x() const
Definition: VectorI.H:75
void insertPoints(const point2DField &points, const scalar nearness)
Create the initial mesh from the given internal points.
label maxBoundaryConformingIter() const
Return the maximum number of boundary conformation iterations.
Definition: cv2DControlsI.H:98
static const char nl
Definition: Ostream.H:265
void write() const
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
bool mkDir(const fileName &, mode_t=0777)
Make a directory and return an error if it could not be created.
Definition: POSIX.C:289
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
messageStream Warning
vector2DField point2DField
point2DField is a vector2DField.
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
void writeTriangles(const fileName &fName, bool internalOnly) const
Write triangles as .obj file.
vector point
Point is a vector.
Definition: point.H:41
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
volScalarField & p
scalar randomPerturbation() const
Return the random perturbation factor.
Definition: cv2DControlsI.H:92
void removeSurfacePointPairs()
Remove the point-pairs introduced by insertSurfacePointPairs.
void boundaryConform()
Insert point-pairs where there are protrusions into.
Namespace for OpenFOAM.