snappySnapDriverFeature.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2015 OpenFOAM Foundation
9  Copyright (C) 2015-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "snappySnapDriver.H"
30 #include "polyTopoChange.H"
31 #include "syncTools.H"
32 #include "fvMesh.H"
33 #include "OBJstream.H"
34 #include "motionSmoother.H"
35 #include "refinementSurfaces.H"
36 #include "refinementFeatures.H"
37 #include "unitConversion.H"
38 #include "plane.H"
39 #include "featureEdgeMesh.H"
40 #include "treeDataPoint.H"
41 #include "indexedOctree.H"
42 #include "snapParameters.H"
43 #include "PatchTools.H"
44 #include "pyramidPointFaceRef.H"
45 #include "localPointRegion.H"
46 
47 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
48 
49 namespace Foam
50 {
51  template<class T>
53  {
54  public:
55 
56  void operator()(List<T>& x, const List<T>& y) const
57  {
58  label sz = x.size();
59  x.setSize(sz+y.size());
60  forAll(y, i)
61  {
62  x[sz++] = y[i];
63  }
64  }
65  };
66 }
67 
68 
69 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
70 
71 bool Foam::snappySnapDriver::isFeaturePoint
72 (
73  const scalar featureCos,
74  const indirectPrimitivePatch& pp,
75  const bitSet& isFeatureEdge,
76  const label pointi
77 ) const
78 {
79  const pointField& points = pp.localPoints();
80  const edgeList& edges = pp.edges();
81  const labelList& pEdges = pp.pointEdges()[pointi];
82 
83  label nFeatEdges = 0;
84 
85  forAll(pEdges, i)
86  {
87  if (isFeatureEdge[pEdges[i]])
88  {
89  nFeatEdges++;
90 
91  for (label j = i+1; j < pEdges.size(); j++)
92  {
93  if (isFeatureEdge[pEdges[j]])
94  {
95  const edge& ei = edges[pEdges[i]];
96  const edge& ej = edges[pEdges[j]];
97 
98  const point& p = points[pointi];
99  const point& pi = points[ei.otherVertex(pointi)];
100  const point& pj = points[ej.otherVertex(pointi)];
101 
102  vector vi = p-pi;
103  scalar viMag = mag(vi);
104 
105  vector vj = pj-p;
106  scalar vjMag = mag(vj);
107 
108  if
109  (
110  viMag > SMALL
111  && vjMag > SMALL
112  && ((vi/viMag & vj/vjMag) < featureCos)
113  )
114  {
115  return true;
116  }
117  }
118  }
119  }
120  }
121 
122  if (nFeatEdges == 1)
123  {
124  // End of feature-edge string
125  return true;
126  }
127 
128  return false;
129 }
130 
131 
132 void Foam::snappySnapDriver::smoothAndConstrain
133 (
134  const bitSet& isPatchMasterEdge,
135  const indirectPrimitivePatch& pp,
136  const labelList& meshEdges,
137  const List<pointConstraint>& constraints,
138  vectorField& disp
139 ) const
140 {
141  const fvMesh& mesh = meshRefiner_.mesh();
142 
143  for (label avgIter = 0; avgIter < 20; avgIter++)
144  {
145  // Calculate average displacement of neighbours
146  // - unconstrained (i.e. surface) points use average of all
147  // neighbouring points
148  // - from testing it has been observed that it is not beneficial
149  // to have edge constrained points use average of all edge or point
150  // constrained neighbours since they're already attracted to
151  // the nearest point on the feature.
152  // Having them attract to point-constrained neighbours does not
153  // make sense either since there is usually just one of them so
154  // it severely distorts it.
155  // - same for feature points. They are already attracted to the
156  // nearest feature point.
157 
158  vectorField dispSum(pp.nPoints(), Zero);
159  labelList dispCount(pp.nPoints(), Zero);
160 
161  const labelListList& pointEdges = pp.pointEdges();
162  const edgeList& edges = pp.edges();
163 
164  forAll(pointEdges, pointi)
165  {
166  const labelList& pEdges = pointEdges[pointi];
167 
168  label nConstraints = constraints[pointi].first();
169 
170  if (nConstraints <= 1)
171  {
172  forAll(pEdges, i)
173  {
174  label edgei = pEdges[i];
175 
176  if (isPatchMasterEdge[edgei])
177  {
178  label nbrPointi = edges[edgei].otherVertex(pointi);
179  if (constraints[nbrPointi].first() >= nConstraints)
180  {
181  dispSum[pointi] += disp[nbrPointi];
182  dispCount[pointi]++;
183  }
184  }
185  }
186  }
187  }
188 
190  (
191  mesh,
192  pp.meshPoints(),
193  dispSum,
194  plusEqOp<point>(),
195  vector::zero,
197  );
199  (
200  mesh,
201  pp.meshPoints(),
202  dispCount,
203  plusEqOp<label>(),
204  label(0),
206  );
207 
208  // Constraints
209  forAll(constraints, pointi)
210  {
211  if (dispCount[pointi] > 0)
212  {
213  // Mix my displacement with neighbours' displacement
214  disp[pointi] =
215  0.5
216  *(disp[pointi] + dispSum[pointi]/dispCount[pointi]);
217  }
218  }
219  }
220 }
221 
222 
223 void Foam::snappySnapDriver::calcNearestFace
224 (
225  const label iter,
226  const indirectPrimitivePatch& pp,
227  const scalarField& faceSnapDist,
228  vectorField& faceDisp,
229  vectorField& faceSurfaceNormal,
230  labelList& faceSurfaceGlobalRegion
231  //vectorField& faceRotation
232 ) const
233 {
234  const fvMesh& mesh = meshRefiner_.mesh();
235  const refinementSurfaces& surfaces = meshRefiner_.surfaces();
236 
237  // Displacement and orientation per pp face.
238  faceDisp.setSize(pp.size());
239  faceDisp = Zero;
240  faceSurfaceNormal.setSize(pp.size());
241  faceSurfaceNormal = Zero;
242  faceSurfaceGlobalRegion.setSize(pp.size());
243  faceSurfaceGlobalRegion = -1;
244 
245  // Divide surfaces into zoned and unzoned
246  const labelList zonedSurfaces =
247  surfaceZonesInfo::getNamedSurfaces(surfaces.surfZones());
248  const labelList unzonedSurfaces =
249  surfaceZonesInfo::getUnnamedSurfaces(surfaces.surfZones());
250 
251  // Per pp face the current surface snapped to
252  labelList snapSurf(pp.size(), -1);
253 
254 
255  // Do zoned surfaces
256  // ~~~~~~~~~~~~~~~~~
257  // Zoned faces only attract to corresponding surface
258 
259  // Extract faces per zone
260  const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
261 
262  forAll(zonedSurfaces, i)
263  {
264  label zoneSurfi = zonedSurfaces[i];
265 
266  const wordList& faceZoneNames = surfZones[zoneSurfi].faceZoneNames();
267 
268  // Get indices of faces on pp that are also in zone
269  DynamicList<label> ppFaces;
270  DynamicList<label> meshFaces;
271  forAll(faceZoneNames, fzi)
272  {
273  const word& faceZoneName = faceZoneNames[fzi];
274  label zonei = mesh.faceZones().findZoneID(faceZoneName);
275  if (zonei == -1)
276  {
278  << "Problem. Cannot find zone " << faceZoneName
279  << exit(FatalError);
280  }
281  const faceZone& fZone = mesh.faceZones()[zonei];
282  const bitSet isZonedFace(mesh.nFaces(), fZone);
283 
284  ppFaces.reserve(ppFaces.capacity()+fZone.size());
285  meshFaces.reserve(meshFaces.capacity()+fZone.size());
286 
287  forAll(pp.addressing(), i)
288  {
289  if (isZonedFace[pp.addressing()[i]])
290  {
291  snapSurf[i] = zoneSurfi;
292  ppFaces.append(i);
293  meshFaces.append(pp.addressing()[i]);
294  }
295  }
296 
297  //Pout<< "For faceZone " << fZone.name()
298  // << " found " << ppFaces.size() << " out of " << pp.size()
299  // << endl;
300  }
301 
302  pointField fc
303  (
305  (
306  IndirectList<face>(mesh.faces(), meshFaces),
307  mesh.points()
308  ).faceCentres()
309  );
310 
311  List<pointIndexHit> hitInfo;
312  labelList hitSurface;
313  labelList hitRegion;
314  vectorField hitNormal;
315  surfaces.findNearestRegion
316  (
317  labelList(1, zoneSurfi),
318  fc,
319  sqr(scalarField(faceSnapDist, ppFaces)),// sqr of attract dist
320  hitSurface,
321  hitInfo,
322  hitRegion,
323  hitNormal
324  );
325 
326  forAll(hitInfo, hiti)
327  {
328  if (hitInfo[hiti].hit())
329  {
330  label facei = ppFaces[hiti];
331  faceDisp[facei] = hitInfo[hiti].hitPoint() - fc[hiti];
332  faceSurfaceNormal[facei] = hitNormal[hiti];
333  faceSurfaceGlobalRegion[facei] = surfaces.globalRegion
334  (
335  hitSurface[hiti],
336  hitRegion[hiti]
337  );
338  }
339  else
340  {
341  static label nWarn = 0;
342 
343  if (nWarn < 100)
344  {
346  << "Did not find surface near face centre " << fc[hiti]
347  << endl;
348  nWarn++;
349  if (nWarn == 100)
350  {
352  << "Reached warning limit " << nWarn
353  << ". Suppressing further warnings." << endl;
354  }
355  }
356  }
357  }
358  }
359 
360 
361  // Do unzoned surfaces
362  // ~~~~~~~~~~~~~~~~~~~
363  // Unzoned faces attract to any unzoned surface
364 
365  DynamicList<label> ppFaces(pp.size());
366  DynamicList<label> meshFaces(pp.size());
367  forAll(pp.addressing(), i)
368  {
369  if (snapSurf[i] == -1)
370  {
371  ppFaces.append(i);
372  meshFaces.append(pp.addressing()[i]);
373  }
374  }
375  //Pout<< "Found " << ppFaces.size() << " unzoned faces out of "
376  // << pp.size() << endl;
377 
378  pointField fc
379  (
381  (
382  IndirectList<face>(mesh.faces(), meshFaces),
383  mesh.points()
384  ).faceCentres()
385  );
386 
387  List<pointIndexHit> hitInfo;
388  labelList hitSurface;
389  labelList hitRegion;
390  vectorField hitNormal;
391  surfaces.findNearestRegion
392  (
393  unzonedSurfaces,
394  fc,
395  sqr(scalarField(faceSnapDist, ppFaces)),// sqr of attract dist
396  hitSurface,
397  hitInfo,
398  hitRegion,
399  hitNormal
400  );
401 
402  forAll(hitInfo, hiti)
403  {
404  if (hitInfo[hiti].hit())
405  {
406  label facei = ppFaces[hiti];
407  faceDisp[facei] = hitInfo[hiti].hitPoint() - fc[hiti];
408  faceSurfaceNormal[facei] = hitNormal[hiti];
409  faceSurfaceGlobalRegion[facei] = surfaces.globalRegion
410  (
411  hitSurface[hiti],
412  hitRegion[hiti]
413  );
414  }
415  else
416  {
417  static label nWarn = 0;
418 
419  if (nWarn < 100)
420  {
422  << "Did not find surface near face centre " << fc[hiti]
423  << endl;
424 
425  nWarn++;
426  if (nWarn == 100)
427  {
429  << "Reached warning limit " << nWarn
430  << ". Suppressing further warnings." << endl;
431  }
432  }
433  }
434  }
435 
436 
439  //
441  //faceRotation.setSize(pp.size());
442  //faceRotation = Zero;
443  //
444  //forAll(faceRotation, facei)
445  //{
446  // // Note: extend to >180 degrees checking
447  // faceRotation[facei] =
448  // pp.faceNormals()[facei]
449  // ^ faceSurfaceNormal[facei];
450  //}
451  //
452  //if (debug&meshRefinement::ATTRACTION)
453  //{
454  // dumpMove
455  // (
456  // mesh.time().path()
457  // / "faceDisp_" + name(iter) + ".obj",
458  // pp.faceCentres(),
459  // pp.faceCentres() + faceDisp
460  // );
461  // dumpMove
462  // (
463  // mesh.time().path()
464  // / "faceRotation_" + name(iter) + ".obj",
465  // pp.faceCentres(),
466  // pp.faceCentres() + faceRotation
467  // );
468  //}
469 }
470 
471 
472 // Collect (possibly remote) per point data of all surrounding faces
473 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
474 // - faceSurfaceNormal
475 // - faceDisp
476 // - faceCentres&faceNormal
477 void Foam::snappySnapDriver::calcNearestFacePointProperties
478 (
479  const label iter,
480  const indirectPrimitivePatch& pp,
481 
482  const vectorField& faceDisp,
483  const vectorField& faceSurfaceNormal,
484  const labelList& faceSurfaceGlobalRegion,
485 
486  List<List<point>>& pointFaceSurfNormals,
487  List<List<point>>& pointFaceDisp,
488  List<List<point>>& pointFaceCentres,
489  List<labelList>& pointFacePatchID
490 ) const
491 {
492  const fvMesh& mesh = meshRefiner_.mesh();
493 
494  const bitSet isMasterFace(syncTools::getMasterFaces(mesh));
495 
496 
497  // For now just get all surrounding face data. Expensive - should just
498  // store and sync data on coupled points only
499  // (see e.g PatchToolsNormals.C)
500 
501  pointFaceSurfNormals.setSize(pp.nPoints());
502  pointFaceDisp.setSize(pp.nPoints());
503  pointFaceCentres.setSize(pp.nPoints());
504  pointFacePatchID.setSize(pp.nPoints());
505 
506  // Fill local data
507  forAll(pp.pointFaces(), pointi)
508  {
509  const labelList& pFaces = pp.pointFaces()[pointi];
510 
511  // Count valid face normals
512  label nFaces = 0;
513  forAll(pFaces, i)
514  {
515  label facei = pFaces[i];
516  if (isMasterFace[facei] && faceSurfaceGlobalRegion[facei] != -1)
517  {
518  nFaces++;
519  }
520  }
521 
522 
523  List<point>& pNormals = pointFaceSurfNormals[pointi];
524  pNormals.setSize(nFaces);
525  List<point>& pDisp = pointFaceDisp[pointi];
526  pDisp.setSize(nFaces);
527  List<point>& pFc = pointFaceCentres[pointi];
528  pFc.setSize(nFaces);
529  labelList& pFid = pointFacePatchID[pointi];
530  pFid.setSize(nFaces);
531 
532  nFaces = 0;
533  forAll(pFaces, i)
534  {
535  label facei = pFaces[i];
536  label globalRegioni = faceSurfaceGlobalRegion[facei];
537 
538  if (isMasterFace[facei] && globalRegioni != -1)
539  {
540  pNormals[nFaces] = faceSurfaceNormal[facei];
541  pDisp[nFaces] = faceDisp[facei];
542  pFc[nFaces] = pp.faceCentres()[facei];
543  pFid[nFaces] = globalToMasterPatch_[globalRegioni];
544  nFaces++;
545  }
546  }
547  }
548 
549 
550  // Collect additionally 'normal' boundary faces for boundaryPoints of pp
551  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
552  // points on the boundary of pp should pick up non-pp normals
553  // as well for the feature-reconstruction to behave correctly.
554  // (the movement is already constrained outside correctly so it
555  // is only that the unconstrained attraction vector is calculated
556  // correctly)
557  {
558  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
559  labelList patchID(pbm.patchID());
560 
561  // Unmark all non-coupled boundary faces
562  forAll(pbm, patchi)
563  {
564  const polyPatch& pp = pbm[patchi];
565 
566  if (pp.coupled() || isA<emptyPolyPatch>(pp))
567  {
568  forAll(pp, i)
569  {
570  label meshFacei = pp.start()+i;
571  patchID[meshFacei-mesh.nInternalFaces()] = -1;
572  }
573  }
574  }
575 
576  // Remove any meshed faces
577  forAll(pp.addressing(), i)
578  {
579  label meshFacei = pp.addressing()[i];
580  patchID[meshFacei-mesh.nInternalFaces()] = -1;
581  }
582 
583 
584 
585  // See if edge of pp uses any non-meshed boundary faces. If so add the
586  // boundary face as additional constraint. Note that we account for
587  // both 'real' boundary edges and boundary edge of baffles
588 
589  const labelList bafflePair
590  (
592  );
593 
594 
595  // Mark all points on 'boundary' edges
596  bitSet isBoundaryPoint(pp.nPoints());
597 
598  const labelListList& edgeFaces = pp.edgeFaces();
599  const edgeList& edges = pp.edges();
600 
601  forAll(edgeFaces, edgei)
602  {
603  const edge& e = edges[edgei];
604  const labelList& eFaces = edgeFaces[edgei];
605 
606  if (eFaces.size() == 1)
607  {
608  // 'real' boundary edge
609  isBoundaryPoint.set(e[0]);
610  isBoundaryPoint.set(e[1]);
611  }
612  else if (eFaces.size() == 2 && bafflePair[eFaces[0]] == eFaces[1])
613  {
614  // 'baffle' boundary edge
615  isBoundaryPoint.set(e[0]);
616  isBoundaryPoint.set(e[1]);
617  }
618  }
619 
620 
621  // Construct labelList equivalent of meshPointMap
622  labelList meshToPatchPoint(mesh.nPoints(), -1);
623  forAll(pp.meshPoints(), pointi)
624  {
625  meshToPatchPoint[pp.meshPoints()[pointi]] = pointi;
626  }
627 
628  forAll(patchID, bFacei)
629  {
630  label patchi = patchID[bFacei];
631 
632  if (patchi != -1)
633  {
634  label facei = mesh.nInternalFaces()+bFacei;
635  const face& f = mesh.faces()[facei];
636 
637  forAll(f, fp)
638  {
639  label pointi = meshToPatchPoint[f[fp]];
640 
641  if (pointi != -1 && isBoundaryPoint.test(pointi))
642  {
643  List<point>& pNormals = pointFaceSurfNormals[pointi];
644  List<point>& pDisp = pointFaceDisp[pointi];
645  List<point>& pFc = pointFaceCentres[pointi];
646  labelList& pFid = pointFacePatchID[pointi];
647 
648  const point& pt = mesh.points()[f[fp]];
649  vector fn = mesh.faceAreas()[facei];
650 
651  pNormals.append(fn/mag(fn));
652  pDisp.append(mesh.faceCentres()[facei]-pt);
653  pFc.append(mesh.faceCentres()[facei]);
654  pFid.append(patchi);
655  }
656  }
657  }
658  }
659  }
660 
662  (
663  mesh,
664  pp.meshPoints(),
665  pointFaceSurfNormals,
666  listPlusEqOp<point>(),
667  List<point>(),
669  );
671  (
672  mesh,
673  pp.meshPoints(),
674  pointFaceDisp,
675  listPlusEqOp<point>(),
676  List<point>(),
678  );
679 
680  {
681  // Make into displacement before synchronising to avoid any problems
682  // with parallel cyclics
683  pointField localPoints(pp.points(), pp.meshPoints());
684  forAll(pointFaceCentres, pointi)
685  {
686  const point& pt = pp.points()[pp.meshPoints()[pointi]];
687 
688  List<point>& pFc = pointFaceCentres[pointi];
689  for (point& p : pFc)
690  {
691  p -= pt;
692  }
693  }
695  (
696  mesh,
697  pp.meshPoints(),
698  pointFaceCentres,
699  listPlusEqOp<point>(),
700  List<point>(),
702  );
703  forAll(pointFaceCentres, pointi)
704  {
705  const point& pt = pp.points()[pp.meshPoints()[pointi]];
706 
707  List<point>& pFc = pointFaceCentres[pointi];
708  for (point& p : pFc)
709  {
710  p += pt;
711  }
712  }
713  }
714 
716  (
717  mesh,
718  pp.meshPoints(),
719  pointFacePatchID,
720  listPlusEqOp<label>(),
721  List<label>()
722  );
723 
724 
725  // Sort the data according to the face centres. This is only so we get
726  // consistent behaviour serial and parallel.
727  labelList visitOrder;
728  forAll(pointFaceDisp, pointi)
729  {
730  List<point>& pNormals = pointFaceSurfNormals[pointi];
731  List<point>& pDisp = pointFaceDisp[pointi];
732  List<point>& pFc = pointFaceCentres[pointi];
733  labelList& pFid = pointFacePatchID[pointi];
734 
735  sortedOrder(mag(pFc)(), visitOrder);
736 
737  pNormals = List<point>(pNormals, visitOrder);
738  pDisp = List<point>(pDisp, visitOrder);
739  pFc = List<point>(pFc, visitOrder);
740  pFid = labelUIndList(pFid, visitOrder)();
741  }
742 }
743 
744 
745 // Gets passed in offset to nearest point on feature edge. Calculates
746 // if the point has a different number of faces on either side of the feature
747 // and if so attracts the point to that non-dominant plane.
748 void Foam::snappySnapDriver::correctAttraction
749 (
750  const DynamicList<point>& surfacePoints,
751  const DynamicList<label>& surfaceCounts,
752  const point& edgePt,
753  const vector& edgeNormal, // normalised normal
754  const point& pt,
755 
756  vector& edgeOffset // offset from pt to point on edge
757 ) const
758 {
759  // Tangential component along edge
760  scalar tang = ((pt-edgePt)&edgeNormal);
761 
762  labelList order(sortedOrder(surfaceCounts));
763 
764  if (order[0] < order[1])
765  {
766  // There is a non-dominant plane. Use the point on the plane to
767  // attract to.
768  vector attractD = surfacePoints[order[0]]-edgePt;
769  // Tangential component along edge
770  scalar tang2 = (attractD&edgeNormal);
771  // Normal component
772  attractD -= tang2*edgeNormal;
773  // Calculate fraction of normal distances
774  scalar magAttractD = mag(attractD);
775  scalar fraction = magAttractD/(magAttractD+mag(edgeOffset));
776 
777  point linePt =
778  edgePt
779  + ((1.0-fraction)*tang2 + fraction*tang)*edgeNormal;
780  edgeOffset = linePt-pt;
781  }
782 }
783 
784 
785 Foam::pointIndexHit Foam::snappySnapDriver::findMultiPatchPoint
786 (
787  const point& pt,
788  const labelList& patchIDs,
789  const List<point>& faceCentres
790 ) const
791 {
792  // Determine if multiple patchIDs
793  if (patchIDs.size())
794  {
795  label patch0 = patchIDs[0];
796 
797  for (label i = 1; i < patchIDs.size(); i++)
798  {
799  if (patchIDs[i] != patch0)
800  {
801  return pointIndexHit(true, pt, labelMax);
802  }
803  }
804  }
805  return pointIndexHit(false, Zero, labelMax);
806 }
807 
808 
809 Foam::label Foam::snappySnapDriver::findNormal
810 (
811  const scalar featureCos,
812  const vector& n,
813  const DynamicList<vector>& surfaceNormals
814 ) const
815 {
816  label index = -1;
817 
818  forAll(surfaceNormals, j)
819  {
820  scalar cosAngle = (n&surfaceNormals[j]);
821 
822  if
823  (
824  (cosAngle >= featureCos)
825  || (cosAngle < (-1+0.001)) // triangle baffles
826  )
827  {
828  index = j;
829  break;
830  }
831  }
832  return index;
833 }
834 
835 
836 // Detect multiple patches. Returns pointIndexHit:
837 // - false, index=-1 : single patch
838 // - true , index=0 : multiple patches but on different normals planes
839 // (so geometric feature edge is also a region edge)
840 // - true , index=1 : multiple patches on same normals plane i.e. flat region
841 // edge
842 Foam::pointIndexHit Foam::snappySnapDriver::findMultiPatchPoint
843 (
844  const point& pt,
845  const labelList& patchIDs,
846  const DynamicList<vector>& surfaceNormals,
847  const labelList& faceToNormalBin
848 ) const
849 {
850  if (patchIDs.empty())
851  {
852  return pointIndexHit(false, pt, -1);
853  }
854 
855  // Detect single patch situation (to avoid allocation)
856  label patch0 = patchIDs[0];
857 
858  for (label i = 1; i < patchIDs.size(); i++)
859  {
860  if (patchIDs[i] != patch0)
861  {
862  patch0 = -1;
863  break;
864  }
865  }
866 
867  if (patch0 >= 0)
868  {
869  // Single patch
870  return pointIndexHit(false, pt, -1);
871  }
872  else
873  {
874  if (surfaceNormals.size() == 1)
875  {
876  // Same normals plane, flat region edge.
877  return pointIndexHit(true, pt, 1);
878  }
879  else
880  {
881  // Detect per normals bin
882  labelList normalToPatch(surfaceNormals.size(), -1);
883  forAll(faceToNormalBin, i)
884  {
885  if (faceToNormalBin[i] != -1)
886  {
887  label& patch = normalToPatch[faceToNormalBin[i]];
888  if (patch == -1)
889  {
890  // First occurrence
891  patch = patchIDs[i];
892  }
893  else if (patch == -2)
894  {
895  // Already marked as being on multiple patches
896  }
897  else if (patch != patchIDs[i])
898  {
899  // Mark as being on multiple patches
900  patch = -2;
901  }
902  }
903  }
904 
905  forAll(normalToPatch, normali)
906  {
907  if (normalToPatch[normali] == -2)
908  {
909  // Multiple patches on same normals plane, flat region
910  // edge
911  return pointIndexHit(true, pt, 1);
912  }
913  }
914 
915  // All patches on either side of geometric feature anyway
916  return pointIndexHit(true, pt, 0);
917  }
918  }
919 }
920 
921 
922 void Foam::snappySnapDriver::writeStats
923 (
924  const indirectPrimitivePatch& pp,
925  const bitSet& isPatchMasterPoint,
926  const List<pointConstraint>& patchConstraints
927 ) const
928 {
929  label nMasterPoints = 0;
930  label nPlanar = 0;
931  label nEdge = 0;
932  label nPoint = 0;
933 
934  forAll(patchConstraints, pointi)
935  {
936  if (isPatchMasterPoint[pointi])
937  {
938  nMasterPoints++;
939 
940  if (patchConstraints[pointi].first() == 1)
941  {
942  nPlanar++;
943  }
944  else if (patchConstraints[pointi].first() == 2)
945  {
946  nEdge++;
947  }
948  else if (patchConstraints[pointi].first() == 3)
949  {
950  nPoint++;
951  }
952  }
953  }
954 
955  reduce(nMasterPoints, sumOp<label>());
956  reduce(nPlanar, sumOp<label>());
957  reduce(nEdge, sumOp<label>());
958  reduce(nPoint, sumOp<label>());
959  Info<< "total master points :" << nMasterPoints
960  << " of which attracted to :" << nl
961  << " feature point : " << nPoint << nl
962  << " feature edge : " << nEdge << nl
963  << " nearest surface : " << nPlanar << nl
964  << " rest : " << nMasterPoints-nPoint-nEdge-nPlanar
965  << nl
966  << endl;
967 }
968 
969 
970 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
971 (
972  const label iter,
973  const scalar featureCos,
974 
975  const indirectPrimitivePatch& pp,
976  const scalarField& snapDist,
977  const vectorField& nearestDisp,
978  const label pointi,
979 
980  const List<List<point>>& pointFaceSurfNormals,
981  const List<List<point>>& pointFaceDisp,
982  const List<List<point>>& pointFaceCentres,
983  const labelListList& pointFacePatchID,
984 
985  DynamicList<point>& surfacePoints,
986  DynamicList<vector>& surfaceNormals,
987  labelList& faceToNormalBin,
988 
989  vector& patchAttraction,
990  pointConstraint& patchConstraint
991 ) const
992 {
993  patchAttraction = Zero;
994  patchConstraint = pointConstraint();
995 
996  const List<point>& pfSurfNormals = pointFaceSurfNormals[pointi];
997  const List<point>& pfDisp = pointFaceDisp[pointi];
998  const List<point>& pfCentres = pointFaceCentres[pointi];
999 
1000  // Bin according to surface normal
1001  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1002 
1003  //- Bins of differing normals:
1004  // - one normal : flat(tish) surface
1005  // - two normals : geometric feature edge
1006  // - three normals: geometric feature point
1007  // - four normals : too complex a feature
1008  surfacePoints.clear();
1009  surfaceNormals.clear();
1010 
1011  //- From face to above normals bin
1012  faceToNormalBin.setSize(pfDisp.size());
1013  faceToNormalBin = -1;
1014 
1015  forAll(pfSurfNormals, i)
1016  {
1017  const point& fc = pfCentres[i];
1018  const vector& fSNormal = pfSurfNormals[i];
1019  const vector& fDisp = pfDisp[i];
1020 
1021  // What to do with very far attraction? For now just ignore the face
1022  if (magSqr(fDisp) < sqr(snapDist[pointi]) && mag(fSNormal) > VSMALL)
1023  {
1024  const point pt = fc + fDisp;
1025 
1026  // Do we already have surface normal?
1027  faceToNormalBin[i] = findNormal
1028  (
1029  featureCos,
1030  fSNormal,
1031  surfaceNormals
1032  );
1033 
1034  if (faceToNormalBin[i] != -1)
1035  {
1036  // Same normal
1037  }
1038  else
1039  {
1040  // Now check if the planes go through the same edge or point
1041 
1042  if (surfacePoints.size() <= 1)
1043  {
1044  surfacePoints.append(pt);
1045  faceToNormalBin[i] = surfaceNormals.size();
1046  surfaceNormals.append(fSNormal);
1047  }
1048  else if (surfacePoints.size() == 2)
1049  {
1050  plane pl0(surfacePoints[0], surfaceNormals[0]);
1051  plane pl1(surfacePoints[1], surfaceNormals[1]);
1052  plane::ray r(pl0.planeIntersect(pl1));
1053  vector featureNormal = r.dir() / mag(r.dir());
1054 
1055  if (mag(fSNormal&featureNormal) >= 0.001)
1056  {
1057  // Definitely makes a feature point
1058  surfacePoints.append(pt);
1059  faceToNormalBin[i] = surfaceNormals.size();
1060  surfaceNormals.append(fSNormal);
1061  }
1062  }
1063  else if (surfacePoints.size() == 3)
1064  {
1065  // Have already feature point. See if this new plane is
1066  // the same point or not.
1067  plane pl0(surfacePoints[0], surfaceNormals[0]);
1068  plane pl1(surfacePoints[1], surfaceNormals[1]);
1069  plane pl2(surfacePoints[2], surfaceNormals[2]);
1070  point p012(pl0.planePlaneIntersect(pl1, pl2));
1071 
1072  plane::ray r(pl0.planeIntersect(pl1));
1073  vector featureNormal = r.dir() / mag(r.dir());
1074  if (mag(fSNormal&featureNormal) >= 0.001)
1075  {
1076  plane pl3(pt, fSNormal);
1077  point p013(pl0.planePlaneIntersect(pl1, pl3));
1078 
1079  if (mag(p012-p013) > snapDist[pointi])
1080  {
1081  // Different feature point
1082  surfacePoints.append(pt);
1083  faceToNormalBin[i] = surfaceNormals.size();
1084  surfaceNormals.append(fSNormal);
1085  }
1086  }
1087  }
1088  }
1089  }
1090  }
1091 
1092 
1093  const point& pt = pp.localPoints()[pointi];
1094 
1095  // Check the number of directions
1096  if (surfaceNormals.size() == 1)
1097  {
1098  // Normal distance to plane
1099  vector d =
1100  ((surfacePoints[0]-pt) & surfaceNormals[0])
1101  *surfaceNormals[0];
1102 
1103  // Trim to snap distance
1104  if (magSqr(d) > sqr(snapDist[pointi]))
1105  {
1106  d *= Foam::sqrt(sqr(snapDist[pointi])/magSqr(d));
1107  }
1108 
1109  patchAttraction = d;
1110 
1111  // Store constraints
1112  patchConstraint.applyConstraint(surfaceNormals[0]);
1113  }
1114  else if (surfaceNormals.size() == 2)
1115  {
1116  plane pl0(surfacePoints[0], surfaceNormals[0]);
1117  plane pl1(surfacePoints[1], surfaceNormals[1]);
1118  plane::ray r(pl0.planeIntersect(pl1));
1119  vector n = r.dir() / mag(r.dir());
1120 
1121  // Get nearest point on infinite ray
1122  vector d = r.refPoint()-pt;
1123  d -= (d&n)*n;
1124 
1125  // Trim to snap distance
1126  if (magSqr(d) > sqr(snapDist[pointi]))
1127  {
1128  d *= Foam::sqrt(sqr(snapDist[pointi])/magSqr(d));
1129  }
1130 
1131  patchAttraction = d;
1132 
1133  // Store constraints
1134  patchConstraint.applyConstraint(surfaceNormals[0]);
1135  patchConstraint.applyConstraint(surfaceNormals[1]);
1136  }
1137  else if (surfaceNormals.size() == 3)
1138  {
1139  // Calculate point from the faces.
1140  plane pl0(surfacePoints[0], surfaceNormals[0]);
1141  plane pl1(surfacePoints[1], surfaceNormals[1]);
1142  plane pl2(surfacePoints[2], surfaceNormals[2]);
1143  point cornerPt(pl0.planePlaneIntersect(pl1, pl2));
1144  vector d = cornerPt - pt;
1145 
1146  // Trim to snap distance
1147  if (magSqr(d) > sqr(snapDist[pointi]))
1148  {
1149  d *= Foam::sqrt(sqr(snapDist[pointi])/magSqr(d));
1150  }
1151 
1152  patchAttraction = d;
1153 
1154  // Store constraints
1155  patchConstraint.applyConstraint(surfaceNormals[0]);
1156  patchConstraint.applyConstraint(surfaceNormals[1]);
1157  patchConstraint.applyConstraint(surfaceNormals[2]);
1158  }
1159 }
1160 
1161 
1162 // Special version that calculates attraction in one go
1163 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
1164 (
1165  const label iter,
1166  const scalar featureCos,
1167 
1168  const indirectPrimitivePatch& pp,
1169  const scalarField& snapDist,
1170  const vectorField& nearestDisp,
1171 
1172  const List<List<point>>& pointFaceSurfNormals,
1173  const List<List<point>>& pointFaceDisp,
1174  const List<List<point>>& pointFaceCentres,
1175  const labelListList& pointFacePatchID,
1176 
1177  vectorField& patchAttraction,
1178  List<pointConstraint>& patchConstraints
1179 ) const
1180 {
1181  autoPtr<OBJstream> feStr;
1182  autoPtr<OBJstream> fpStr;
1184  {
1185  feStr.reset
1186  (
1187  new OBJstream
1188  (
1189  meshRefiner_.mesh().time().path()
1190  / "implicitFeatureEdge_" + name(iter) + ".obj"
1191  )
1192  );
1193  Info<< "Dumping implicit feature-edge direction to "
1194  << feStr().name() << endl;
1195 
1196  fpStr.reset
1197  (
1198  new OBJstream
1199  (
1200  meshRefiner_.mesh().time().path()
1201  / "implicitFeaturePoint_" + name(iter) + ".obj"
1202  )
1203  );
1204  Info<< "Dumping implicit feature-point direction to "
1205  << fpStr().name() << endl;
1206  }
1207 
1208 
1209  DynamicList<point> surfacePoints(4);
1210  DynamicList<vector> surfaceNormals(4);
1211  labelList faceToNormalBin;
1212 
1213  forAll(pp.localPoints(), pointi)
1214  {
1215  vector attraction = Zero;
1216  pointConstraint constraint;
1217 
1218  featureAttractionUsingReconstruction
1219  (
1220  iter,
1221  featureCos,
1222 
1223  pp,
1224  snapDist,
1225  nearestDisp,
1226 
1227  pointi,
1228 
1229  pointFaceSurfNormals,
1230  pointFaceDisp,
1231  pointFaceCentres,
1232  pointFacePatchID,
1233 
1234  surfacePoints,
1235  surfaceNormals,
1236  faceToNormalBin,
1237 
1238  attraction,
1239  constraint
1240  );
1241 
1242  if
1243  (
1244  (constraint.first() > patchConstraints[pointi].first())
1245  || (
1246  (constraint.first() == patchConstraints[pointi].first())
1247  && (magSqr(attraction) < magSqr(patchAttraction[pointi]))
1248  )
1249  )
1250  {
1251  patchAttraction[pointi] = attraction;
1252  patchConstraints[pointi] = constraint;
1253 
1254  const point& pt = pp.localPoints()[pointi];
1255 
1256  if (feStr && patchConstraints[pointi].first() == 2)
1257  {
1258  feStr().write(linePointRef(pt, pt+patchAttraction[pointi]));
1259  }
1260  else if (fpStr && patchConstraints[pointi].first() == 3)
1261  {
1262  fpStr().write(linePointRef(pt, pt+patchAttraction[pointi]));
1263  }
1264  }
1265  }
1266 }
1267 
1268 
1269 void Foam::snappySnapDriver::stringFeatureEdges
1270 (
1271  const label iter,
1272  const scalar featureCos,
1273 
1274  const indirectPrimitivePatch& pp,
1275  const scalarField& snapDist,
1276 
1277  const vectorField& rawPatchAttraction,
1278  const List<pointConstraint>& rawPatchConstraints,
1279 
1280  vectorField& patchAttraction,
1281  List<pointConstraint>& patchConstraints
1282 ) const
1283 {
1284  // Snap edges to feature edges
1285  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1286  // Walk existing edges and snap remaining ones (that are marked as
1287  // feature edges in rawPatchConstraints)
1288 
1289  // What this does is fill in any faces where not all points
1290  // on the face are being attracted:
1291  /*
1292  +
1293  / \
1294  / \
1295  ---+ +---
1296  \ /
1297  \ /
1298  +
1299  */
1300  // so the top and bottom will never get attracted since the nearest
1301  // back from the feature edge will always be one of the left or right
1302  // points since the face is diamond like. So here we walk the feature edges
1303  // and add any non-attracted points.
1304 
1305 
1306  while (true)
1307  {
1308  label nChanged = 0;
1309 
1310  const labelListList& pointEdges = pp.pointEdges();
1311  forAll(pointEdges, pointi)
1312  {
1313  if (patchConstraints[pointi].first() == 2)
1314  {
1315  const point& pt = pp.localPoints()[pointi];
1316  const labelList& pEdges = pointEdges[pointi];
1317  const vector& featVec = patchConstraints[pointi].second();
1318 
1319  // Detect whether there are edges in both directions.
1320  // (direction along the feature edge that is)
1321  bool hasPos = false;
1322  bool hasNeg = false;
1323 
1324  forAll(pEdges, pEdgei)
1325  {
1326  const edge& e = pp.edges()[pEdges[pEdgei]];
1327  label nbrPointi = e.otherVertex(pointi);
1328 
1329  if (patchConstraints[nbrPointi].first() > 1)
1330  {
1331  const point& nbrPt = pp.localPoints()[nbrPointi];
1332  const point featPt =
1333  nbrPt + patchAttraction[nbrPointi];
1334  const scalar cosAngle = (featVec & (featPt-pt));
1335 
1336  if (cosAngle > 0)
1337  {
1338  hasPos = true;
1339  }
1340  else
1341  {
1342  hasNeg = true;
1343  }
1344  }
1345  }
1346 
1347  if (!hasPos || !hasNeg)
1348  {
1349  //Pout<< "**Detected feature string end at "
1350  // << pp.localPoints()[pointi] << endl;
1351 
1352  // No string. Assign best choice on either side
1353  label bestPosPointi = -1;
1354  scalar minPosDistSqr = GREAT;
1355  label bestNegPointi = -1;
1356  scalar minNegDistSqr = GREAT;
1357 
1358  forAll(pEdges, pEdgei)
1359  {
1360  const edge& e = pp.edges()[pEdges[pEdgei]];
1361  label nbrPointi = e.otherVertex(pointi);
1362 
1363  if
1364  (
1365  patchConstraints[nbrPointi].first() <= 1
1366  && rawPatchConstraints[nbrPointi].first() > 1
1367  )
1368  {
1369  const vector& nbrFeatVec =
1370  rawPatchConstraints[pointi].second();
1371 
1372  if (mag(featVec&nbrFeatVec) > featureCos)
1373  {
1374  // nbrPointi attracted to sameish feature
1375  // Note: also check on position.
1376 
1377  scalar d2 = magSqr
1378  (
1379  rawPatchAttraction[nbrPointi]
1380  );
1381 
1382  const point featPt =
1383  pp.localPoints()[nbrPointi]
1384  + rawPatchAttraction[nbrPointi];
1385  const scalar cosAngle =
1386  (featVec & (featPt-pt));
1387 
1388  if (cosAngle > 0)
1389  {
1390  if (!hasPos && d2 < minPosDistSqr)
1391  {
1392  minPosDistSqr = d2;
1393  bestPosPointi = nbrPointi;
1394  }
1395  }
1396  else
1397  {
1398  if (!hasNeg && d2 < minNegDistSqr)
1399  {
1400  minNegDistSqr = d2;
1401  bestNegPointi = nbrPointi;
1402  }
1403  }
1404  }
1405  }
1406  }
1407 
1408  if (bestPosPointi != -1)
1409  {
1410  // Use reconstructed-feature attraction. Use only
1411  // part of it since not sure...
1412  //const point& bestPt =
1413  // pp.localPoints()[bestPosPointi];
1414  //Pout<< "**Overriding point " << bestPt
1415  // << " on reconstructed feature edge at "
1416  // << rawPatchAttraction[bestPosPointi]+bestPt
1417  // << " to attracted-to-feature-edge." << endl;
1418  patchAttraction[bestPosPointi] =
1419  0.5*rawPatchAttraction[bestPosPointi];
1420  patchConstraints[bestPosPointi] =
1421  rawPatchConstraints[bestPosPointi];
1422 
1423  nChanged++;
1424  }
1425  if (bestNegPointi != -1)
1426  {
1427  // Use reconstructed-feature attraction. Use only
1428  // part of it since not sure...
1429  //const point& bestPt =
1430  // pp.localPoints()[bestNegPointi];
1431  //Pout<< "**Overriding point " << bestPt
1432  // << " on reconstructed feature edge at "
1433  // << rawPatchAttraction[bestNegPointi]+bestPt
1434  // << " to attracted-to-feature-edge." << endl;
1435  patchAttraction[bestNegPointi] =
1436  0.5*rawPatchAttraction[bestNegPointi];
1437  patchConstraints[bestNegPointi] =
1438  rawPatchConstraints[bestNegPointi];
1439 
1440  nChanged++;
1441  }
1442  }
1443  }
1444  }
1445 
1446 
1447  reduce(nChanged, sumOp<label>());
1448  Info<< "Stringing feature edges : changed " << nChanged << " points"
1449  << endl;
1450  if (nChanged == 0)
1451  {
1452  break;
1453  }
1454  }
1455 }
1456 
1457 
1458 void Foam::snappySnapDriver::releasePointsNextToMultiPatch
1459 (
1460  const label iter,
1461  const scalar featureCos,
1462 
1463  const indirectPrimitivePatch& pp,
1464  const scalarField& snapDist,
1465 
1466  const List<List<point>>& pointFaceCentres,
1467  const labelListList& pointFacePatchID,
1468 
1469  const vectorField& rawPatchAttraction,
1470  const List<pointConstraint>& rawPatchConstraints,
1471 
1472  vectorField& patchAttraction,
1473  List<pointConstraint>& patchConstraints
1474 ) const
1475 {
1476  autoPtr<OBJstream> multiPatchStr;
1478  {
1479  multiPatchStr.reset
1480  (
1481  new OBJstream
1482  (
1483  meshRefiner_.mesh().time().path()
1484  / "multiPatch_" + name(iter) + ".obj"
1485  )
1486  );
1487  Info<< "Dumping removed constraints due to same-face"
1488  << " multi-patch points to "
1489  << multiPatchStr().name() << endl;
1490  }
1491 
1492 
1493  // 1. Mark points on multiple patches
1494  bitSet isMultiPatchPoint(pp.size());
1495 
1496  forAll(pointFacePatchID, pointi)
1497  {
1498  pointIndexHit multiPatchPt = findMultiPatchPoint
1499  (
1500  pp.localPoints()[pointi],
1501  pointFacePatchID[pointi],
1502  pointFaceCentres[pointi]
1503  );
1504  isMultiPatchPoint.set(pointi, multiPatchPt.hit());
1505  }
1506 
1507  // 2. Make sure multi-patch points are also attracted
1508  forAll(isMultiPatchPoint, pointi)
1509  {
1510  if (isMultiPatchPoint.test(pointi))
1511  {
1512  if
1513  (
1514  patchConstraints[pointi].first() <= 1
1515  && rawPatchConstraints[pointi].first() > 1
1516  )
1517  {
1518  patchAttraction[pointi] = rawPatchAttraction[pointi];
1519  patchConstraints[pointi] = rawPatchConstraints[pointi];
1520 
1521  //if (multiPatchStr)
1522  //{
1523  // Pout<< "Adding constraint on multiPatchPoint:"
1524  // << pp.localPoints()[pointi]
1525  // << " constraint:" << patchConstraints[pointi]
1526  // << " attraction:" << patchAttraction[pointi]
1527  // << endl;
1528  //}
1529  }
1530  }
1531  }
1532 
1533  // Up to here it is all parallel ok.
1534 
1535 
1536  // 3. Knock out any attraction on faces with multi-patch points
1537  label nChanged = 0;
1538  forAll(pp.localFaces(), facei)
1539  {
1540  const face& f = pp.localFaces()[facei];
1541 
1542  label nMultiPatchPoints = 0;
1543  forAll(f, fp)
1544  {
1545  label pointi = f[fp];
1546  if
1547  (
1548  isMultiPatchPoint.test(pointi)
1549  && patchConstraints[pointi].first() > 1
1550  )
1551  {
1552  ++nMultiPatchPoints;
1553  }
1554  }
1555 
1556  if (nMultiPatchPoints > 0)
1557  {
1558  forAll(f, fp)
1559  {
1560  label pointi = f[fp];
1561  if
1562  (
1563  !isMultiPatchPoint.test(pointi)
1564  && patchConstraints[pointi].first() > 1
1565  )
1566  {
1567  //Pout<< "Knocking out constraint"
1568  // << " on non-multiPatchPoint:"
1569  // << pp.localPoints()[pointi] << endl;
1570  patchAttraction[pointi] = Zero;
1571  patchConstraints[pointi] = pointConstraint();
1572  nChanged++;
1573 
1574  if (multiPatchStr)
1575  {
1576  multiPatchStr().write(pp.localPoints()[pointi]);
1577  }
1578  }
1579  }
1580  }
1581  }
1582 
1583  reduce(nChanged, sumOp<label>());
1584  Info<< "Removing constraints near multi-patch points : changed "
1585  << nChanged << " points" << endl;
1586 }
1587 
1588 
1589 Foam::labelPair Foam::snappySnapDriver::findDiagonalAttraction
1590 (
1591  const indirectPrimitivePatch& pp,
1592  const vectorField& patchAttraction,
1593  const List<pointConstraint>& patchConstraints,
1594  const label facei
1595 ) const
1596 {
1597  const face& f = pp.localFaces()[facei];
1598  // For now just detect any attraction. Improve this to look at
1599  // actual attraction position and orientation
1600 
1601  labelPair attractIndices(-1, -1);
1602 
1603  if (f.size() >= 4)
1604  {
1605  for (label startFp = 0; startFp < f.size()-2; startFp++)
1606  {
1607  label minFp = f.rcIndex(startFp);
1608 
1609  for
1610  (
1611  label endFp = f.fcIndex(f.fcIndex(startFp));
1612  endFp < f.size() && endFp != minFp;
1613  endFp++
1614  )
1615  {
1616  if
1617  (
1618  patchConstraints[f[startFp]].first() >= 2
1619  && patchConstraints[f[endFp]].first() >= 2
1620  )
1621  {
1622  attractIndices = labelPair(startFp, endFp);
1623  break;
1624  }
1625  }
1626  }
1627  }
1628  return attractIndices;
1629 }
1630 
1631 
1632 bool Foam::snappySnapDriver::isSplitAlignedWithFeature
1633 (
1634  const scalar featureCos,
1635  const point& p0,
1636  const pointConstraint& pc0,
1637  const point& p1,
1638  const pointConstraint& pc1
1639 ) const
1640 {
1641  vector d(p1-p0);
1642  scalar magD = mag(d);
1643  if (magD < VSMALL)
1644  {
1645  // Two diagonal points already colocated?
1646  return false;
1647  }
1648  else
1649  {
1650  d /= magD;
1651 
1652  // Is diagonal d aligned with at least one of the feature
1653  // edges?
1654 
1655  if (pc0.first() == 2 && mag(d & pc0.second()) > featureCos)
1656  {
1657  return true;
1658  }
1659  else if (pc1.first() == 2 && mag(d & pc1.second()) > featureCos)
1660  {
1661  return true;
1662  }
1663  else
1664  {
1665  return false;
1666  }
1667  }
1668 }
1669 
1670 
1671 // Is situation very concave
1672 bool Foam::snappySnapDriver::isConcave
1673 (
1674  const point& c0,
1675  const vector& area0,
1676  const point& c1,
1677  const vector& area1,
1678  const scalar concaveCos
1679 ) const
1680 {
1681  vector n0 = area0;
1682  scalar magN0 = mag(n0);
1683  if (magN0 < VSMALL)
1684  {
1685  // Zero area face. What to return? For now disable splitting.
1686  return true;
1687  }
1688  n0 /= magN0;
1689 
1690  // Distance from c1 to plane of face0
1691  scalar d = (c1-c0)&n0;
1692 
1693  if (d <= 0)
1694  {
1695  // Convex (face1 centre on 'inside' of face0)
1696  return false;
1697  }
1698  else
1699  {
1700  // Is a bit or very concave?
1701  vector n1 = area1;
1702  scalar magN1 = mag(n1);
1703  if (magN1 < VSMALL)
1704  {
1705  // Zero area face. See above
1706  return true;
1707  }
1708  n1 /= magN1;
1709 
1710  if ((n0&n1) < concaveCos)
1711  {
1712  return true;
1713  }
1714  else
1715  {
1716  return false;
1717  }
1718  }
1719 }
1720 
1721 
1722 Foam::labelPair Foam::snappySnapDriver::findDiagonalAttraction
1723 (
1724  const scalar featureCos,
1725  const scalar concaveCos,
1726  const scalar minAreaRatio,
1727  const indirectPrimitivePatch& pp,
1728  const vectorField& patchAttr,
1729  const List<pointConstraint>& patchConstraints,
1730  const vectorField& nearestAttr,
1731  const vectorField& nearestNormal,
1732  const label facei,
1733 
1734  DynamicField<point>& points0,
1735  DynamicField<point>& points1
1736 ) const
1737 {
1738  const face& localF = pp.localFaces()[facei];
1739 
1740  labelPair attractIndices(-1, -1);
1741 
1742  if (localF.size() >= 4)
1743  {
1744  const pointField& localPts = pp.localPoints();
1745 
1749  //const polyMesh& mesh = meshRefiner_.mesh();
1750  //label meshFacei = pp.addressing()[facei];
1751  //const face& meshF = mesh.faces()[meshFacei];
1752  //label celli = mesh.faceOwner()[meshFacei];
1753  //const labelList& cPoints = mesh.cellPoints(celli);
1754  //
1755  //point cc(mesh.points()[meshF[0]]);
1756  //for (label i = 1; i < meshF.size(); i++)
1757  //{
1758  // cc += mesh.points()[meshF[i]]+patchAttr[localF[i]];
1759  //}
1760  //forAll(cPoints, i)
1761  //{
1762  // label pointi = cPoints[i];
1763  // if (!meshF.found(pointi))
1764  // {
1765  // cc += mesh.points()[pointi];
1766  // }
1767  //}
1768  //cc /= cPoints.size();
1770  //
1771  //const scalar vol = pyrVol(pp, patchAttr, localF, cc);
1772  //const scalar area = localF.mag(localPts);
1773 
1774 
1775 
1776  // Try all diagonal cuts
1777  // ~~~~~~~~~~~~~~~~~~~~~
1778 
1779  face f0(3);
1780  face f1(3);
1781 
1782  for (label startFp = 0; startFp < localF.size()-2; startFp++)
1783  {
1784  label minFp = localF.rcIndex(startFp);
1785 
1786  for
1787  (
1788  label endFp = localF.fcIndex(localF.fcIndex(startFp));
1789  endFp < localF.size() && endFp != minFp;
1790  endFp++
1791  )
1792  {
1793  label startPti = localF[startFp];
1794  label endPti = localF[endFp];
1795 
1796  const pointConstraint& startPc = patchConstraints[startPti];
1797  const pointConstraint& endPc = patchConstraints[endPti];
1798 
1799  if (startPc.first() >= 2 && endPc.first() >= 2)
1800  {
1801  if (startPc.first() == 2 || endPc.first() == 2)
1802  {
1803  // Check if
1804  // - sameish feature edge normal
1805  // - diagonal aligned with feature edge normal
1806  point start = localPts[startPti]+patchAttr[startPti];
1807  point end = localPts[endPti]+patchAttr[endPti];
1808 
1809  if
1810  (
1811  !isSplitAlignedWithFeature
1812  (
1813  featureCos,
1814  start,
1815  startPc,
1816  end,
1817  endPc
1818  )
1819  )
1820  {
1821  // Attract to different features. No need to
1822  // introduce split
1823  continue;
1824  }
1825  }
1826 
1827 
1828 
1829  // Form two faces
1830  // ~~~~~~~~~~~~~~
1831  // Predict position of faces. End points of the faces
1832  // attract to the feature
1833  // and all the other points just attract to the nearest
1834 
1835  // face0
1836 
1837  f0.setSize(endFp-startFp+1);
1838  label i0 = 0;
1839  for (label fp = startFp; fp <= endFp; fp++)
1840  {
1841  f0[i0++] = localF[fp];
1842  }
1843 
1844  // Get compact face and points
1845  const face compact0(identity(f0.size()));
1846  points0.clear();
1847  points0.append(localPts[f0[0]] + patchAttr[f0[0]]);
1848  for (label fp=1; fp < f0.size()-1; fp++)
1849  {
1850  label pi = f0[fp];
1851  points0.append(localPts[pi] + nearestAttr[pi]);
1852  }
1853  points0.append
1854  (
1855  localPts[f0.last()] + patchAttr[f0.last()]
1856  );
1857 
1858 
1859  // face1
1860 
1861  f1.setSize(localF.size()+2-f0.size());
1862  label i1 = 0;
1863  for
1864  (
1865  label fp = endFp;
1866  fp != startFp;
1867  fp = localF.fcIndex(fp)
1868  )
1869  {
1870  f1[i1++] = localF[fp];
1871  }
1872  f1[i1++] = localF[startFp];
1873 
1874 
1875  // Get compact face and points
1876  const face compact1(identity(f1.size()));
1877  points1.clear();
1878  points1.append(localPts[f1[0]] + patchAttr[f1[0]]);
1879  for (label fp=1; fp < f1.size()-1; fp++)
1880  {
1881  label pi = f1[fp];
1882  points1.append(localPts[pi] + nearestAttr[pi]);
1883  }
1884  points1.append
1885  (
1886  localPts[f1.last()] + patchAttr[f1.last()]
1887  );
1888 
1889 
1890 
1891  // Avoid splitting concave faces
1892  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1893 
1894  if
1895  (
1896  isConcave
1897  (
1898  compact0.centre(points0),
1899  compact0.areaNormal(points0),
1900  compact1.centre(points1),
1901  compact1.areaNormal(points1),
1902  concaveCos
1903  )
1904  )
1905  {
1907  //Pout<< "# f0:" << f0 << endl;
1908  //forAll(p, i)
1909  //{
1910  // meshTools::writeOBJ(Pout, points0[i]);
1911  //}
1912  //Pout<< "# f1:" << f1 << endl;
1913  //forAll(p, i)
1914  //{
1915  // meshTools::writeOBJ(Pout, points1[i]);
1916  //}
1917  }
1918  else
1919  {
1920  // Existing areas
1921  const scalar area0 = f0.mag(localPts);
1922  const scalar area1 = f1.mag(localPts);
1923 
1924  if
1925  (
1926  area0/area1 >= minAreaRatio
1927  && area1/area0 >= minAreaRatio
1928  )
1929  {
1930  attractIndices = labelPair(startFp, endFp);
1931  }
1932  }
1933  }
1934  }
1935  }
1936  }
1937  return attractIndices;
1938 }
1939 
1940 
1941 void Foam::snappySnapDriver::splitDiagonals
1942 (
1943  const scalar featureCos,
1944  const scalar concaveCos,
1945  const scalar minAreaRatio,
1946 
1947  const indirectPrimitivePatch& pp,
1948  const vectorField& nearestAttraction,
1949  const vectorField& nearestNormal,
1950 
1951  vectorField& patchAttraction,
1952  List<pointConstraint>& patchConstraints,
1953  DynamicList<label>& splitFaces,
1954  DynamicList<labelPair>& splits
1955 ) const
1956 {
1957  const labelList& bFaces = pp.addressing();
1958 
1959  splitFaces.clear();
1960  splitFaces.setCapacity(bFaces.size());
1961  splits.clear();
1962  splits.setCapacity(bFaces.size());
1963 
1964 
1965  // Work arrays for storing points of face
1966  DynamicField<point> facePoints0;
1967  DynamicField<point> facePoints1;
1968 
1969  forAll(bFaces, facei)
1970  {
1971  const labelPair split
1972  (
1973  findDiagonalAttraction
1974  (
1975  featureCos,
1976  concaveCos,
1977  minAreaRatio,
1978 
1979  pp,
1980  patchAttraction,
1981  patchConstraints,
1982 
1983  nearestAttraction,
1984  nearestNormal,
1985  facei,
1986 
1987  facePoints0,
1988  facePoints1
1989  )
1990  );
1991 
1992  if (split != labelPair(-1, -1))
1993  {
1994  splitFaces.append(bFaces[facei]);
1995  splits.append(split);
1996 
1997  const face& f = pp.localFaces()[facei];
1998 
1999  // Knock out other attractions on face
2000  forAll(f, fp)
2001  {
2002  // Knock out any other constraints
2003  if
2004  (
2005  fp != split[0]
2006  && fp != split[1]
2007  && patchConstraints[f[fp]].first() >= 2
2008  )
2009  {
2010  patchConstraints[f[fp]] = pointConstraint
2011  (
2012  Tuple2<label, vector>
2013  (
2014  1,
2015  nearestNormal[f[fp]]
2016  )
2017  );
2018  patchAttraction[f[fp]] = nearestAttraction[f[fp]];
2019  }
2020  }
2021  }
2022  }
2023 }
2024 
2025 
2026 void Foam::snappySnapDriver::avoidDiagonalAttraction
2027 (
2028  const label iter,
2029  const scalar featureCos,
2030 
2031  const indirectPrimitivePatch& pp,
2032 
2033  vectorField& patchAttraction,
2034  List<pointConstraint>& patchConstraints
2035 ) const
2036 {
2037  forAll(pp.localFaces(), facei)
2038  {
2039  const face& f = pp.localFaces()[facei];
2040 
2041  labelPair diag = findDiagonalAttraction
2042  (
2043  pp,
2044  patchAttraction,
2045  patchConstraints,
2046  facei
2047  );
2048 
2049  if (diag[0] != -1 && diag[1] != -1)
2050  {
2051  // Found two diagonal points that being attracted.
2052  // For now just attract my one to the average of those.
2053  const label i0 = f[diag[0]];
2054  const point pt0 =
2055  pp.localPoints()[i0]+patchAttraction[i0];
2056  const label i1 = f[diag[1]];
2057  const point pt1 =
2058  pp.localPoints()[i1]+patchAttraction[i1];
2059  const point mid = 0.5*(pt0+pt1);
2060 
2061  const scalar cosAngle = mag
2062  (
2063  patchConstraints[i0].second()
2064  & patchConstraints[i1].second()
2065  );
2066 
2067  //Pout<< "Found diagonal attraction at indices:"
2068  // << diag[0]
2069  // << " and " << diag[1]
2070  // << " with cosAngle:" << cosAngle
2071  // << " mid:" << mid << endl;
2072 
2073  if (cosAngle > featureCos)
2074  {
2075  // The two feature edges are roughly in the same direction.
2076  // Add the nearest of the other points in the face as
2077  // attractor
2078  label minFp = -1;
2079  scalar minDistSqr = GREAT;
2080  forAll(f, fp)
2081  {
2082  label pointi = f[fp];
2083  if (patchConstraints[pointi].first() <= 1)
2084  {
2085  const point& pt = pp.localPoints()[pointi];
2086  scalar distSqr = magSqr(mid-pt);
2087  if (distSqr < minDistSqr)
2088  {
2089  minFp = fp;
2090  }
2091  }
2092  }
2093  if (minFp != -1)
2094  {
2095  label minPointi = f[minFp];
2096  patchAttraction[minPointi] =
2097  mid-pp.localPoints()[minPointi];
2098  patchConstraints[minPointi] = patchConstraints[f[diag[0]]];
2099  }
2100  }
2101  else
2102  {
2103  //Pout<< "Diagonal attractors at" << nl
2104  // << " pt0:" << pt0
2105  // << " constraint:"
2106  // << patchConstraints[i0].second() << nl
2107  // << " pt1:" << pt1
2108  // << " constraint:"
2109  // << patchConstraints[i1].second() << nl
2110  // << " make too large an angle:"
2111  // << mag
2112  // (
2113  // patchConstraints[i0].second()
2114  // & patchConstraints[i1].second()
2115  // )
2116  // << endl;
2117  }
2118  }
2119  }
2120 }
2121 
2122 
2124 Foam::snappySnapDriver::findNearFeatureEdge
2125 (
2126  const bool isRegionEdge,
2127 
2128  const indirectPrimitivePatch& pp,
2129  const scalarField& snapDist,
2130  const label pointi,
2131  const point& estimatedPt,
2132 
2133  List<List<DynamicList<point>>>& edgeAttractors,
2134  List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2135  vectorField& patchAttraction,
2136  List<pointConstraint>& patchConstraints
2137 ) const
2138 {
2139  const refinementFeatures& features = meshRefiner_.features();
2140 
2141  labelList nearEdgeFeat;
2142  List<pointIndexHit> nearEdgeInfo;
2143  vectorField nearNormal;
2144 
2145  if (isRegionEdge)
2146  {
2147  features.findNearestRegionEdge
2148  (
2149  pointField(1, estimatedPt),
2150  scalarField(1, sqr(snapDist[pointi])),
2151  nearEdgeFeat,
2152  nearEdgeInfo,
2153  nearNormal
2154  );
2155  }
2156  else
2157  {
2158  features.findNearestEdge
2159  (
2160  pointField(1, estimatedPt),
2161  scalarField(1, sqr(snapDist[pointi])),
2162  nearEdgeFeat,
2163  nearEdgeInfo,
2164  nearNormal
2165  );
2166  }
2167 
2168  const pointIndexHit& nearInfo = nearEdgeInfo[0];
2169  label feati = nearEdgeFeat[0];
2170 
2171  if (nearInfo.hit())
2172  {
2173  // So we have a point on the feature edge. Use this
2174  // instead of our estimate from planes.
2175  edgeAttractors[feati][nearInfo.index()].append
2176  (
2177  nearInfo.hitPoint()
2178  );
2179  pointConstraint c(Tuple2<label, vector>(2, nearNormal[0]));
2180  edgeConstraints[feati][nearInfo.index()].append(c);
2181 
2182  // Store for later use
2183  patchAttraction[pointi] = nearInfo.hitPoint()-pp.localPoints()[pointi];
2184  patchConstraints[pointi] = c;
2185  }
2186  return Tuple2<label, pointIndexHit>(feati, nearInfo);
2187 }
2188 
2189 
2191 Foam::snappySnapDriver::findNearFeaturePoint
2192 (
2193  const bool isRegionPoint,
2194 
2195  const indirectPrimitivePatch& pp,
2196  const scalarField& snapDist,
2197  const label pointi,
2198  const point& estimatedPt,
2199 
2200  // Feature-point to pp point
2201  List<labelList>& pointAttractor,
2202  List<List<pointConstraint>>& pointConstraints,
2203  // Feature-edge to pp point
2204  List<List<DynamicList<point>>>& edgeAttractors,
2205  List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2206  // pp point to nearest feature
2207  vectorField& patchAttraction,
2208  List<pointConstraint>& patchConstraints
2209 ) const
2210 {
2211  const refinementFeatures& features = meshRefiner_.features();
2212 
2213  // Search for for featurePoints only! This ignores any non-feature points.
2214 
2215  labelList nearFeat;
2216  List<pointIndexHit> nearInfo;
2217  features.findNearestPoint
2218  (
2219  pointField(1, estimatedPt),
2220  scalarField(1, sqr(snapDist[pointi])),
2221  nearFeat,
2222  nearInfo
2223  );
2224 
2225  label feati = nearFeat[0];
2226 
2227  if (feati != -1)
2228  {
2229  const point& pt = pp.localPoints()[pointi];
2230 
2231  label featPointi = nearInfo[0].index();
2232  const point& featPt = nearInfo[0].hitPoint();
2233  scalar distSqr = magSqr(featPt-pt);
2234 
2235  // Check if already attracted
2236  label oldPointi = pointAttractor[feati][featPointi];
2237 
2238  if (oldPointi != -1)
2239  {
2240  // Check distance
2241  if (distSqr >= magSqr(featPt-pp.localPoints()[oldPointi]))
2242  {
2243  // oldPointi nearest. Keep.
2244  feati = -1;
2245  featPointi = -1;
2246  }
2247  else
2248  {
2249  // Current pointi nearer.
2250  pointAttractor[feati][featPointi] = pointi;
2251  pointConstraints[feati][featPointi].first() = 3;
2252  pointConstraints[feati][featPointi].second() = Zero;
2253 
2254  // Store for later use
2255  patchAttraction[pointi] = featPt-pt;
2256  patchConstraints[pointi] = pointConstraints[feati][featPointi];
2257 
2258  // Reset oldPointi to nearest on feature edge
2259  patchAttraction[oldPointi] = Zero;
2260  patchConstraints[oldPointi] = pointConstraint();
2261 
2262  findNearFeatureEdge
2263  (
2264  isRegionPoint, // search region edges only
2265 
2266  pp,
2267  snapDist,
2268  oldPointi,
2269  pp.localPoints()[oldPointi],
2270 
2271  edgeAttractors,
2272  edgeConstraints,
2273  patchAttraction,
2274  patchConstraints
2275  );
2276  }
2277  }
2278  else
2279  {
2280  // Current pointi nearer.
2281  pointAttractor[feati][featPointi] = pointi;
2282  pointConstraints[feati][featPointi].first() = 3;
2283  pointConstraints[feati][featPointi].second() = Zero;
2284 
2285  // Store for later use
2286  patchAttraction[pointi] = featPt-pt;
2287  patchConstraints[pointi] = pointConstraints[feati][featPointi];
2288  }
2289  }
2290 
2291  return Tuple2<label, pointIndexHit>(feati, nearInfo[0]);
2292 }
2293 
2294 
2295 // Determines for every pp point - that is on multiple faces that form
2296 // a feature - the nearest feature edge/point.
2297 void Foam::snappySnapDriver::determineFeatures
2298 (
2299  const label iter,
2300  const scalar featureCos,
2301  const bool multiRegionFeatureSnap,
2302 
2303  const indirectPrimitivePatch& pp,
2304  const scalarField& snapDist,
2305  const vectorField& nearestDisp,
2306 
2307  const List<List<point>>& pointFaceSurfNormals,
2308  const List<List<point>>& pointFaceDisp,
2309  const List<List<point>>& pointFaceCentres,
2310  const labelListList& pointFacePatchID,
2311 
2312  // Feature-point to pp point
2313  List<labelList>& pointAttractor,
2314  List<List<pointConstraint>>& pointConstraints,
2315  // Feature-edge to pp point
2316  List<List<DynamicList<point>>>& edgeAttractors,
2317  List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2318  // pp point to nearest feature
2319  vectorField& patchAttraction,
2320  List<pointConstraint>& patchConstraints
2321 ) const
2322 {
2323  autoPtr<OBJstream> featureEdgeStr;
2324  autoPtr<OBJstream> missedEdgeStr;
2325  autoPtr<OBJstream> featurePointStr;
2326  autoPtr<OBJstream> missedMP0Str;
2327  autoPtr<OBJstream> missedMP1Str;
2328 
2330  {
2331  featureEdgeStr.reset
2332  (
2333  new OBJstream
2334  (
2335  meshRefiner_.mesh().time().path()
2336  / "featureEdge_" + name(iter) + ".obj"
2337  )
2338  );
2339  Info<< "Dumping feature-edge sampling to "
2340  << featureEdgeStr().name() << endl;
2341 
2342  missedEdgeStr.reset
2343  (
2344  new OBJstream
2345  (
2346  meshRefiner_.mesh().time().path()
2347  / "missedFeatureEdge_" + name(iter) + ".obj"
2348  )
2349  );
2350  Info<< "Dumping feature-edges that are too far away to "
2351  << missedEdgeStr().name() << endl;
2352 
2353  featurePointStr.reset
2354  (
2355  new OBJstream
2356  (
2357  meshRefiner_.mesh().time().path()
2358  / "featurePoint_" + name(iter) + ".obj"
2359  )
2360  );
2361  Info<< "Dumping feature-point sampling to "
2362  << featurePointStr().name() << endl;
2363 
2364  missedMP0Str.reset
2365  (
2366  new OBJstream
2367  (
2368  meshRefiner_.mesh().time().path()
2369  / "missedFeatureEdgeFromMPEdge_" + name(iter) + ".obj"
2370  )
2371  );
2372  Info<< "Dumping region-edges that are too far away to "
2373  << missedMP0Str().name() << endl;
2374 
2375  missedMP1Str.reset
2376  (
2377  new OBJstream
2378  (
2379  meshRefiner_.mesh().time().path()
2380  / "missedFeatureEdgeFromMPPoint_" + name(iter) + ".obj"
2381  )
2382  );
2383  Info<< "Dumping region-points that are too far away to "
2384  << missedMP1Str().name() << endl;
2385  }
2386 
2387 
2388  DynamicList<point> surfacePoints(4);
2389  DynamicList<vector> surfaceNormals(4);
2390  labelList faceToNormalBin;
2391 
2392  forAll(pp.localPoints(), pointi)
2393  {
2394  const point& pt = pp.localPoints()[pointi];
2395 
2396 
2397  // Determine the geometric planes the point is (approximately) on.
2398  // This is returned as a
2399  // - attraction vector
2400  // - and a constraint
2401  // (1: attract to surface, constraint is normal of plane
2402  // 2: attract to feature line, constraint is feature line direction
2403  // 3: attract to feature point, constraint is zero)
2404 
2405  vector attraction = Zero;
2406  pointConstraint constraint;
2407 
2408  featureAttractionUsingReconstruction
2409  (
2410  iter,
2411  featureCos,
2412 
2413  pp,
2414  snapDist,
2415  nearestDisp,
2416 
2417  pointi,
2418 
2419  pointFaceSurfNormals,
2420  pointFaceDisp,
2421  pointFaceCentres,
2422  pointFacePatchID,
2423 
2424  surfacePoints,
2425  surfaceNormals,
2426  faceToNormalBin,
2427 
2428  attraction,
2429  constraint
2430  );
2431 
2432  // Now combine the reconstruction with the current state of the
2433  // point. The logic is quite complicated:
2434  // - the new constraint (from reconstruction) will only win if
2435  // - the constraint is higher (feature-point wins from feature-edge
2436  // etc.)
2437  // - or the constraint is the same but the attraction distance is less
2438  //
2439  // - then this will be combined with explicit searching on the
2440  // features and optionally the analysis of the patches using the
2441  // point. This analysis can do three thing:
2442  // - the point is not on multiple patches
2443  // - the point is on multiple patches but these are also
2444  // different planes (so the region feature is also a geometric
2445  // feature)
2446  // - the point is on multiple patches some of which are on
2447  // the same plane. This is the problem one - do we assume it is
2448  // an additional constraint (feat edge upgraded to region point,
2449  // see below)?
2450  //
2451  // Reconstruction MultiRegionFeatureSnap Attraction
2452  // ------- ---------------------- -----------
2453  // surface false surface
2454  // surface true region edge
2455  // feat edge false feat edge
2456  // feat edge true and no planar regions feat edge
2457  // feat edge true and yes planar regions region point
2458  // feat point false feat point
2459  // feat point true region point
2460 
2461 
2462  if
2463  (
2464  (constraint.first() > patchConstraints[pointi].first())
2465  || (
2466  (constraint.first() == patchConstraints[pointi].first())
2467  && (magSqr(attraction) < magSqr(patchAttraction[pointi]))
2468  )
2469  )
2470  {
2471  patchAttraction[pointi] = attraction;
2472  patchConstraints[pointi] = constraint;
2473 
2474  // Check the number of directions
2475  if (patchConstraints[pointi].first() == 1)
2476  {
2477  // Flat surface. Check for different patchIDs
2478  if (multiRegionFeatureSnap)
2479  {
2480  const point estimatedPt(pt + nearestDisp[pointi]);
2481  pointIndexHit multiPatchPt
2482  (
2483  findMultiPatchPoint
2484  (
2485  estimatedPt,
2486  pointFacePatchID[pointi],
2487  surfaceNormals,
2488  faceToNormalBin
2489  )
2490  );
2491 
2492  if (multiPatchPt.hit())
2493  {
2494  // Behave like when having two surface normals so
2495  // attract to nearest feature edge (with a guess for
2496  // the multipatch point as starting point)
2497  Tuple2<label, pointIndexHit> nearInfo =
2498  findNearFeatureEdge
2499  (
2500  true, // isRegionEdge
2501  pp,
2502  snapDist,
2503  pointi,
2504  multiPatchPt.hitPoint(), // estimatedPt
2505 
2506  edgeAttractors,
2507  edgeConstraints,
2508 
2509  patchAttraction,
2510  patchConstraints
2511  );
2512 
2513  const pointIndexHit& info = nearInfo.second();
2514  if (info.hit())
2515  {
2516  // Dump
2517  if (featureEdgeStr)
2518  {
2519  featureEdgeStr().write
2520  (
2521  linePointRef(pt, info.hitPoint())
2522  );
2523  }
2524  }
2525  else
2526  {
2527  if (missedEdgeStr)
2528  {
2529  missedEdgeStr().write
2530  (
2531  linePointRef(pt, multiPatchPt.hitPoint())
2532  );
2533  }
2534  }
2535  }
2536  }
2537  }
2538  else if (patchConstraints[pointi].first() == 2)
2539  {
2540  // Mark point on the nearest feature edge. Note that we
2541  // only search within the surrounding since the plane
2542  // reconstruction might find a feature where there isn't one.
2543  const point estimatedPt(pt + patchAttraction[pointi]);
2544 
2545  Tuple2<label, pointIndexHit> nearInfo(-1, pointIndexHit());
2546 
2547  // Geometric feature edge. Check for different patchIDs
2548  bool hasSnapped = false;
2549  if (multiRegionFeatureSnap)
2550  {
2551  pointIndexHit multiPatchPt
2552  (
2553  findMultiPatchPoint
2554  (
2555  estimatedPt,
2556  pointFacePatchID[pointi],
2557  surfaceNormals,
2558  faceToNormalBin
2559  )
2560  );
2561  if (multiPatchPt.hit())
2562  {
2563  if (multiPatchPt.index() == 0)
2564  {
2565  // Region edge is also a geometric feature edge
2566  nearInfo = findNearFeatureEdge
2567  (
2568  true, // isRegionEdge
2569  pp,
2570  snapDist,
2571  pointi,
2572  estimatedPt,
2573 
2574  edgeAttractors,
2575  edgeConstraints,
2576 
2577  patchAttraction,
2578  patchConstraints
2579  );
2580  hasSnapped = true;
2581 
2582  // Debug: dump missed feature point
2583  if (missedMP0Str && !nearInfo.second().hit())
2584  {
2585  missedMP0Str().write
2586  (
2587  linePointRef(pt, estimatedPt)
2588  );
2589  }
2590  }
2591  else
2592  {
2593  // One of planes of feature contains multiple
2594  // regions. We assume (contentious!) that the
2595  // separation between
2596  // the regions is not aligned with the geometric
2597  // feature so is an additional constraint on the
2598  // point -> is region-feature-point.
2599  nearInfo = findNearFeaturePoint
2600  (
2601  true, // isRegionPoint
2602  pp,
2603  snapDist,
2604  pointi,
2605  estimatedPt,
2606 
2607  // Feature-point to pp point
2608  pointAttractor,
2609  pointConstraints,
2610  // Feature-edge to pp point
2611  edgeAttractors,
2612  edgeConstraints,
2613  // pp point to nearest feature
2614  patchAttraction,
2615  patchConstraints
2616  );
2617  hasSnapped = true;
2618 
2619  // More contentious: if we don't find
2620  // a near feature point we will never find the
2621  // attraction to a feature edge either since
2622  // the edgeAttractors/edgeConstraints do not get
2623  // filled and we're using reverse attraction
2624  // Note that we're in multiRegionFeatureSnap which
2625  // where findMultiPatchPoint can decide the
2626  // wrong thing. So: if failed finding a near
2627  // feature point try for a feature edge
2628  if (!nearInfo.second().hit())
2629  {
2630  nearInfo = findNearFeatureEdge
2631  (
2632  true, // isRegionEdge
2633  pp,
2634  snapDist,
2635  pointi,
2636  estimatedPt,
2637 
2638  // Feature-edge to pp point
2639  edgeAttractors,
2640  edgeConstraints,
2641  // pp point to nearest feature
2642  patchAttraction,
2643  patchConstraints
2644  );
2645  }
2646 
2647  // Debug: dump missed feature point
2648  if (missedMP1Str && !nearInfo.second().hit())
2649  {
2650  missedMP1Str().write
2651  (
2652  linePointRef(pt, estimatedPt)
2653  );
2654  }
2655  }
2656  }
2657  }
2658 
2659  if (!hasSnapped)
2660  {
2661  // Determine nearest point on feature edge. Store
2662  // constraint
2663  // (calculated from feature edge, alternative would be to
2664  // use constraint calculated from both surfaceNormals)
2665  nearInfo = findNearFeatureEdge
2666  (
2667  false, // isRegionPoint
2668  pp,
2669  snapDist,
2670  pointi,
2671  estimatedPt,
2672 
2673  edgeAttractors,
2674  edgeConstraints,
2675 
2676  patchAttraction,
2677  patchConstraints
2678  );
2679  hasSnapped = true;
2680  }
2681 
2682  // Dump to obj
2683  const pointIndexHit& info = nearInfo.second();
2684  if (info.hit())
2685  {
2686  if
2687  (
2688  featurePointStr
2689  && patchConstraints[pointi].first() == 3
2690  )
2691  {
2692  featurePointStr().write
2693  (
2694  linePointRef(pt, info.hitPoint())
2695  );
2696  }
2697  else if
2698  (
2699  featureEdgeStr
2700  && patchConstraints[pointi].first() == 2
2701  )
2702  {
2703  featureEdgeStr().write
2704  (
2705  linePointRef(pt, info.hitPoint())
2706  );
2707  }
2708  }
2709  else
2710  {
2711  if (missedEdgeStr)
2712  {
2713  missedEdgeStr().write
2714  (
2715  linePointRef(pt, estimatedPt)
2716  );
2717  }
2718  }
2719  }
2720  else if (patchConstraints[pointi].first() == 3)
2721  {
2722  // Mark point on the nearest feature point.
2723  const point estimatedPt(pt + patchAttraction[pointi]);
2724 
2725  Tuple2<label, pointIndexHit> nearInfo(-1, pointIndexHit());
2726 
2727  if (multiRegionFeatureSnap)
2728  {
2729  pointIndexHit multiPatchPt
2730  (
2731  findMultiPatchPoint
2732  (
2733  estimatedPt,
2734  pointFacePatchID[pointi],
2735  surfaceNormals,
2736  faceToNormalBin
2737  )
2738  );
2739  if (multiPatchPt.hit())
2740  {
2741  // Multiple regions
2742  nearInfo = findNearFeaturePoint
2743  (
2744  true, // isRegionPoint
2745  pp,
2746  snapDist,
2747  pointi,
2748  estimatedPt,
2749 
2750  // Feature-point to pp point
2751  pointAttractor,
2752  pointConstraints,
2753  // Feature-edge to pp point
2754  edgeAttractors,
2755  edgeConstraints,
2756  // pp point to nearest feature
2757  patchAttraction,
2758  patchConstraints
2759  );
2760  }
2761  else
2762  {
2763  nearInfo = findNearFeaturePoint
2764  (
2765  false, // isRegionPoint
2766  pp,
2767  snapDist,
2768  pointi,
2769  estimatedPt,
2770 
2771  // Feature-point to pp point
2772  pointAttractor,
2773  pointConstraints,
2774  // Feature-edge to pp point
2775  edgeAttractors,
2776  edgeConstraints,
2777  // pp point to nearest feature
2778  patchAttraction,
2779  patchConstraints
2780  );
2781  }
2782  }
2783  else
2784  {
2785  // No multi-patch snapping
2786  nearInfo = findNearFeaturePoint
2787  (
2788  false, // isRegionPoint
2789  pp,
2790  snapDist,
2791  pointi,
2792  estimatedPt,
2793 
2794  // Feature-point to pp point
2795  pointAttractor,
2796  pointConstraints,
2797  // Feature-edge to pp point
2798  edgeAttractors,
2799  edgeConstraints,
2800  // pp point to nearest feature
2801  patchAttraction,
2802  patchConstraints
2803  );
2804  }
2805 
2806  const pointIndexHit& info = nearInfo.second();
2807  if (featurePointStr && info.hit())
2808  {
2809  featurePointStr().write
2810  (
2811  linePointRef(pt, info.hitPoint())
2812  );
2813  }
2814  }
2815  }
2816  }
2817 }
2818 
2819 
2820 // Baffle handling
2821 // ~~~~~~~~~~~~~~~
2822 // Override pointAttractor, edgeAttractor, patchAttration etc. to
2823 // implement 'baffle' handling.
2824 // Baffle: the mesh pp point originates from a loose standing
2825 // baffle.
2826 // Sampling the surface with the surrounding face-centres only picks up
2827 // a single triangle normal so above determineFeatures will not have
2828 // detected anything. So explicitly pick up feature edges on the pp
2829 // (after duplicating points & smoothing so will already have been
2830 // expanded) and match these to the features.
2831 void Foam::snappySnapDriver::determineBaffleFeatures
2832 (
2833  const label iter,
2834  const bool baffleFeaturePoints,
2835  const scalar featureCos,
2836 
2837  const indirectPrimitivePatch& pp,
2838  const scalarField& snapDist,
2839 
2840  // Feature-point to pp point
2841  List<labelList>& pointAttractor,
2842  List<List<pointConstraint>>& pointConstraints,
2843  // Feature-edge to pp point
2844  List<List<DynamicList<point>>>& edgeAttractors,
2845  List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2846  // pp point to nearest feature
2847  vectorField& patchAttraction,
2848  List<pointConstraint>& patchConstraints
2849 ) const
2850 {
2851  const fvMesh& mesh = meshRefiner_.mesh();
2852  const refinementFeatures& features = meshRefiner_.features();
2853 
2854  // Calculate edge-faces
2855  List<List<point>> edgeFaceNormals(pp.nEdges());
2856 
2857  // Fill local data
2858  forAll(pp.edgeFaces(), edgei)
2859  {
2860  const labelList& eFaces = pp.edgeFaces()[edgei];
2861  List<point>& eFc = edgeFaceNormals[edgei];
2862  eFc.setSize(eFaces.size());
2863  forAll(eFaces, i)
2864  {
2865  label facei = eFaces[i];
2866  eFc[i] = pp.faceNormals()[facei];
2867  }
2868  }
2869 
2870  {
2871  // Precalculate mesh edges for pp.edges.
2872  const labelList meshEdges
2873  (
2874  pp.meshEdges(mesh.edges(), mesh.pointEdges())
2875  );
2876  // Collect all coupled edges. Does not filter duplicates/order
2878  (
2879  mesh,
2880  meshEdges,
2881  edgeFaceNormals,
2882  listPlusEqOp<point>(),
2883  List<point>(),
2885  noOp()
2886  );
2887  }
2888 
2889  // Detect baffle edges. Assume initial mesh will have 0,90 or 180
2890  // (baffle) degree angles so smoothing should make 0,90
2891  // to be less than 90. Choose reasonable value
2892  const scalar baffleFeatureCos = Foam::cos(degToRad(110.0));
2893 
2894 
2895  autoPtr<OBJstream> baffleEdgeStr;
2897  {
2898  baffleEdgeStr.reset
2899  (
2900  new OBJstream
2901  (
2902  meshRefiner_.mesh().time().path()
2903  / "baffleEdge_" + name(iter) + ".obj"
2904  )
2905  );
2906  Info<< nl << "Dumping baffle-edges to "
2907  << baffleEdgeStr().name() << endl;
2908  }
2909 
2910 
2911  // Is edge on baffle
2912  bitSet isBaffleEdge(pp.nEdges());
2913  label nBaffleEdges = 0;
2914  // Is point on
2915  // 0 : baffle-edge (0)
2916  // 1 : baffle-feature-point (1)
2917  // -1 : rest
2918  labelList pointStatus(pp.nPoints(), -1);
2919 
2920  forAll(edgeFaceNormals, edgei)
2921  {
2922  const List<point>& efn = edgeFaceNormals[edgei];
2923 
2924  if (efn.size() == 2 && (efn[0]&efn[1]) < baffleFeatureCos)
2925  {
2926  isBaffleEdge.set(edgei);
2927  ++nBaffleEdges;
2928  const edge& e = pp.edges()[edgei];
2929  pointStatus[e[0]] = 0;
2930  pointStatus[e[1]] = 0;
2931 
2932  if (baffleEdgeStr)
2933  {
2934  const point& p0 = pp.localPoints()[e[0]];
2935  const point& p1 = pp.localPoints()[e[1]];
2936  baffleEdgeStr().write(linePointRef(p0, p1));
2937  }
2938  }
2939  }
2940 
2941  reduce(nBaffleEdges, sumOp<label>());
2942 
2943  Info<< "Detected " << nBaffleEdges
2944  << " baffle edges out of "
2945  << returnReduce(pp.nEdges(), sumOp<label>())
2946  << " edges." << endl;
2947 
2948 
2949  //- Baffle edges will be too ragged to sensibly determine feature points
2950  //forAll(pp.pointEdges(), pointi)
2951  //{
2952  // if
2953  // (
2954  // isFeaturePoint
2955  // (
2956  // featureCos,
2957  // pp,
2958  // isBaffleEdge,
2959  // pointi
2960  // )
2961  // )
2962  // {
2963  // //Pout<< "Detected feature point:" << pp.localPoints()[pointi]
2964  // // << endl;
2965  // //-TEMPORARILY DISABLED:
2966  // //pointStatus[pointi] = 1;
2967  // }
2968  //}
2969 
2970 
2971  label nBafflePoints = 0;
2972  forAll(pointStatus, pointi)
2973  {
2974  if (pointStatus[pointi] != -1)
2975  {
2976  nBafflePoints++;
2977  }
2978  }
2979  reduce(nBafflePoints, sumOp<label>());
2980 
2981 
2982  label nPointAttract = 0;
2983  label nEdgeAttract = 0;
2984 
2985  forAll(pointStatus, pointi)
2986  {
2987  const point& pt = pp.localPoints()[pointi];
2988 
2989  if (pointStatus[pointi] == 0) // baffle edge
2990  {
2991  // 1: attract to near feature edge first
2992 
2993  Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
2994  (
2995  false, // isRegionPoint?
2996  pp,
2997  snapDist,
2998  pointi,
2999  pt,
3000 
3001  edgeAttractors,
3002  edgeConstraints,
3003  patchAttraction,
3004  patchConstraints
3005  );
3006 
3007 
3008  //- MEJ:
3009  // 2: optionally override with nearest feature point.
3010  // On baffles we don't have enough normals to construct a feature
3011  // point so assume all feature edges are close to feature points
3012  if (nearInfo.second().hit())
3013  {
3014  nEdgeAttract++;
3015 
3016  if (baffleFeaturePoints)
3017  {
3018  nearInfo = findNearFeaturePoint
3019  (
3020  false, // isRegionPoint,
3021 
3022  pp,
3023  snapDist,
3024  pointi,
3025  pt, // estimatedPt,
3026 
3027  // Feature-point to pp point
3028  pointAttractor,
3029  pointConstraints,
3030  // Feature-edge to pp point
3031  edgeAttractors,
3032  edgeConstraints,
3033  // pp point to nearest feature
3034  patchAttraction,
3035  patchConstraints
3036  );
3037 
3038  if (nearInfo.first() != -1)
3039  {
3040  nEdgeAttract--;
3041  nPointAttract++;
3042  }
3043  }
3044  }
3045  }
3046  else if (pointStatus[pointi] == 1) // baffle point
3047  {
3048  labelList nearFeat;
3049  List<pointIndexHit> nearInfo;
3050  features.findNearestPoint
3051  (
3052  pointField(1, pt),
3053  scalarField(1, sqr(snapDist[pointi])),
3054  nearFeat,
3055  nearInfo
3056  );
3057 
3058  label feati = nearFeat[0];
3059 
3060  if (feati != -1)
3061  {
3062  nPointAttract++;
3063 
3064  label featPointi = nearInfo[0].index();
3065  const point& featPt = nearInfo[0].hitPoint();
3066  scalar distSqr = magSqr(featPt-pt);
3067 
3068  // Check if already attracted
3069  label oldPointi = pointAttractor[feati][featPointi];
3070 
3071  if
3072  (
3073  oldPointi == -1
3074  || (
3075  distSqr
3076  < magSqr(featPt-pp.localPoints()[oldPointi])
3077  )
3078  )
3079  {
3080  pointAttractor[feati][featPointi] = pointi;
3081  pointConstraints[feati][featPointi].first() = 3;
3082  pointConstraints[feati][featPointi].second() = Zero;
3083 
3084  // Store for later use
3085  patchAttraction[pointi] = featPt-pt;
3086  patchConstraints[pointi] =
3087  pointConstraints[feati][featPointi];
3088 
3089  if (oldPointi != -1)
3090  {
3091  // The current point is closer so wins. Reset
3092  // the old point to attract to nearest edge
3093  // instead.
3094  findNearFeatureEdge
3095  (
3096  false, // isRegionPoint
3097  pp,
3098  snapDist,
3099  oldPointi,
3100  pp.localPoints()[oldPointi],
3101 
3102  edgeAttractors,
3103  edgeConstraints,
3104  patchAttraction,
3105  patchConstraints
3106  );
3107  }
3108  }
3109  else
3110  {
3111  // Make it fall through to check below
3112  feati = -1;
3113  }
3114  }
3115 
3116  // Not found a feature point or another point is already
3117  // closer to that feature
3118  if (feati == -1)
3119  {
3120  //Pout<< "*** Falling back to finding nearest feature"
3121  // << " edge"
3122  // << " for baffle-feature-point " << pt
3123  // << endl;
3124 
3125  Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
3126  (
3127  false, // isRegionPoint
3128  pp,
3129  snapDist,
3130  pointi,
3131  pt, // starting point
3132 
3133  edgeAttractors,
3134  edgeConstraints,
3135  patchAttraction,
3136  patchConstraints
3137  );
3138 
3139  if (nearInfo.first() != -1)
3140  {
3141  nEdgeAttract++;
3142  }
3143  }
3144  }
3145  }
3146 
3147  reduce(nPointAttract, sumOp<label>());
3148  reduce(nEdgeAttract, sumOp<label>());
3149 
3150  Info<< "Baffle points : " << nBafflePoints
3151  << " of which attracted to :" << nl
3152  << " feature point : " << nPointAttract << nl
3153  << " feature edge : " << nEdgeAttract << nl
3154  << " rest : " << nBafflePoints-nPointAttract-nEdgeAttract
3155  << nl
3156  << endl;
3157 }
3158 
3159 
3160 void Foam::snappySnapDriver::reverseAttractMeshPoints
3161 (
3162  const label iter,
3163 
3164  const indirectPrimitivePatch& pp,
3165  const scalarField& snapDist,
3166 
3167  // Feature-point to pp point
3168  const List<labelList>& pointAttractor,
3169  const List<List<pointConstraint>>& pointConstraints,
3170  // Feature-edge to pp point
3171  const List<List<DynamicList<point>>>& edgeAttractors,
3172  const List<List<DynamicList<pointConstraint>>>& edgeConstraints,
3173 
3174  const vectorField& rawPatchAttraction,
3175  const List<pointConstraint>& rawPatchConstraints,
3176 
3177  // pp point to nearest feature
3178  vectorField& patchAttraction,
3179  List<pointConstraint>& patchConstraints
3180 ) const
3181 {
3182  const refinementFeatures& features = meshRefiner_.features();
3183 
3184  // Find nearest mesh point to feature edge
3185  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3186  // Reverse lookup : go through all edgeAttractors and find the
3187  // nearest point on pp
3188 
3189  // Get search domain and extend it a bit
3190  treeBoundBox bb(pp.localPoints());
3191  {
3192  // Random number generator. Bit dodgy since not exactly random ;-)
3193  Random rndGen(65431);
3194 
3195  // Slightly extended bb. Slightly off-centred just so on symmetric
3196  // geometry there are less face/edge aligned items.
3197  bb = bb.extend(rndGen, 1e-4);
3198  bb.min() -= point::uniform(ROOTVSMALL);
3199  bb.max() += point::uniform(ROOTVSMALL);
3200  }
3201 
3202  // Collect candidate points for attraction
3203  DynamicList<label> attractPoints(pp.nPoints());
3204  {
3205  const fvMesh& mesh = meshRefiner_.mesh();
3206 
3207  boolList isFeatureEdgeOrPoint(pp.nPoints(), false);
3208  label nFeats = 0;
3209  forAll(rawPatchConstraints, pointi)
3210  {
3211  if (rawPatchConstraints[pointi].first() >= 2)
3212  {
3213  isFeatureEdgeOrPoint[pointi] = true;
3214  nFeats++;
3215  }
3216  }
3217 
3218  Info<< "Initially selected " << returnReduce(nFeats, sumOp<label>())
3219  << " mesh points out of "
3220  << returnReduce(pp.nPoints(), sumOp<label>())
3221  << " for reverse attraction." << endl;
3222 
3223  // Make sure is synchronised (note: check if constraint is already
3224  // synced in which case this is not needed here)
3226  (
3227  mesh,
3228  pp.meshPoints(),
3229  isFeatureEdgeOrPoint,
3230  orEqOp<bool>(), // combine op
3231  false
3232  );
3233 
3234  for (label nGrow = 0; nGrow < 1; nGrow++)
3235  {
3236  boolList newIsFeatureEdgeOrPoint(isFeatureEdgeOrPoint);
3237 
3238  forAll(pp.localFaces(), facei)
3239  {
3240  const face& f = pp.localFaces()[facei];
3241 
3242  forAll(f, fp)
3243  {
3244  if (isFeatureEdgeOrPoint[f[fp]])
3245  {
3246  // Mark all points on face
3247  forAll(f, fp)
3248  {
3249  newIsFeatureEdgeOrPoint[f[fp]] = true;
3250  }
3251  break;
3252  }
3253  }
3254  }
3255 
3256  isFeatureEdgeOrPoint = newIsFeatureEdgeOrPoint;
3257 
3259  (
3260  mesh,
3261  pp.meshPoints(),
3262  isFeatureEdgeOrPoint,
3263  orEqOp<bool>(), // combine op
3264  false
3265  );
3266  }
3267 
3268 
3269  // Collect attractPoints
3270  forAll(isFeatureEdgeOrPoint, pointi)
3271  {
3272  if (isFeatureEdgeOrPoint[pointi])
3273  {
3274  attractPoints.append(pointi);
3275  }
3276  }
3277 
3278  Info<< "Selected "
3279  << returnReduce(attractPoints.size(), sumOp<label>())
3280  << " mesh points out of "
3281  << returnReduce(pp.nPoints(), sumOp<label>())
3282  << " for reverse attraction." << endl;
3283  }
3284 
3285 
3286  indexedOctree<treeDataPoint> ppTree
3287  (
3288  treeDataPoint(pp.localPoints(), attractPoints),
3289  bb, // overall search domain
3290  8, // maxLevel
3291  10, // leafsize
3292  3.0 // duplicity
3293  );
3294 
3295  // Per mesh point the point on nearest feature edge.
3296  patchAttraction.setSize(pp.nPoints());
3297  patchAttraction = Zero;
3298  patchConstraints.setSize(pp.nPoints());
3299  patchConstraints = pointConstraint();
3300 
3301  forAll(edgeAttractors, feati)
3302  {
3303  const List<DynamicList<point>>& edgeAttr = edgeAttractors[feati];
3304  const List<DynamicList<pointConstraint>>& edgeConstr =
3305  edgeConstraints[feati];
3306 
3307  forAll(edgeAttr, featEdgei)
3308  {
3309  const DynamicList<point>& attr = edgeAttr[featEdgei];
3310  forAll(attr, i)
3311  {
3312  // Find nearest pp point
3313  const point& featPt = attr[i];
3314  pointIndexHit nearInfo = ppTree.findNearest
3315  (
3316  featPt,
3317  sqr(GREAT)
3318  );
3319 
3320  if (nearInfo.hit())
3321  {
3322  label pointi =
3323  ppTree.shapes().pointLabels()[nearInfo.index()];
3324  const point attraction = featPt-pp.localPoints()[pointi];
3325 
3326  // Check if this point is already being attracted. If so
3327  // override it only if nearer.
3328  if
3329  (
3330  patchConstraints[pointi].first() <= 1
3331  || magSqr(attraction) < magSqr(patchAttraction[pointi])
3332  )
3333  {
3334  patchAttraction[pointi] = attraction;
3335  patchConstraints[pointi] = edgeConstr[featEdgei][i];
3336  }
3337  }
3338  else
3339  {
3340  static label nWarn = 0;
3341 
3342  if (nWarn < 100)
3343  {
3345  << "Did not find pp point near " << featPt
3346  << endl;
3347  nWarn++;
3348  if (nWarn == 100)
3349  {
3351  << "Reached warning limit " << nWarn
3352  << ". Suppressing further warnings." << endl;
3353  }
3354  }
3355 
3356  }
3357  }
3358  }
3359  }
3360 
3361 
3362  // Different procs might have different patchAttraction,patchConstraints
3363  // however these only contain geometric information, no topology
3364  // so as long as we synchronise after overriding with feature points
3365  // there is no problem, just possibly a small error.
3366 
3367 
3368  // Find nearest mesh point to feature point
3369  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3370  // (overrides attraction to feature edge)
3371  forAll(pointAttractor, feati)
3372  {
3373  const labelList& pointAttr = pointAttractor[feati];
3374  const List<pointConstraint>& pointConstr = pointConstraints[feati];
3375 
3376  forAll(pointAttr, featPointi)
3377  {
3378  if (pointAttr[featPointi] != -1)
3379  {
3380  const point& featPt = features[feati].points()
3381  [
3382  featPointi
3383  ];
3384 
3385  // Find nearest pp point
3386  pointIndexHit nearInfo = ppTree.findNearest
3387  (
3388  featPt,
3389  sqr(GREAT)
3390  );
3391 
3392  if (nearInfo.hit())
3393  {
3394  label pointi =
3395  ppTree.shapes().pointLabels()[nearInfo.index()];
3396 
3397  const point& pt = pp.localPoints()[pointi];
3398  const point attraction = featPt-pt;
3399 
3400  // - already attracted to feature edge : point always wins
3401  // - already attracted to feature point: nearest wins
3402 
3403  if (patchConstraints[pointi].first() <= 1)
3404  {
3405  patchAttraction[pointi] = attraction;
3406  patchConstraints[pointi] = pointConstr[featPointi];
3407  }
3408  else if (patchConstraints[pointi].first() == 2)
3409  {
3410  patchAttraction[pointi] = attraction;
3411  patchConstraints[pointi] = pointConstr[featPointi];
3412  }
3413  else if (patchConstraints[pointi].first() == 3)
3414  {
3415  // Only if nearer
3416  if
3417  (
3418  magSqr(attraction)
3419  < magSqr(patchAttraction[pointi])
3420  )
3421  {
3422  patchAttraction[pointi] = attraction;
3423  patchConstraints[pointi] =
3424  pointConstr[featPointi];
3425  }
3426  }
3427  }
3428  }
3429  }
3430  }
3431 }
3432 
3433 
3434 void Foam::snappySnapDriver::featureAttractionUsingFeatureEdges
3435 (
3436  const label iter,
3437  const bool multiRegionFeatureSnap,
3438 
3439  const bool detectBaffles,
3440  const bool baffleFeaturePoints,
3441 
3442  const bool releasePoints,
3443  const bool stringFeatures,
3444  const bool avoidDiagonal,
3445 
3446  const scalar featureCos,
3447 
3448  const indirectPrimitivePatch& pp,
3449  const scalarField& snapDist,
3450  const vectorField& nearestDisp,
3451  const vectorField& nearestNormal,
3452 
3453  const List<List<point>>& pointFaceSurfNormals,
3454  const List<List<point>>& pointFaceDisp,
3455  const List<List<point>>& pointFaceCentres,
3456  const labelListList& pointFacePatchID,
3457 
3458  vectorField& patchAttraction,
3459  List<pointConstraint>& patchConstraints
3460 ) const
3461 {
3462  const refinementFeatures& features = meshRefiner_.features();
3463  const fvMesh& mesh = meshRefiner_.mesh();
3464 
3465  const bitSet isPatchMasterPoint
3466  (
3468  (
3469  mesh,
3470  pp.meshPoints()
3471  )
3472  );
3473 
3474 
3475  // Collect ordered attractions on feature edges
3476  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3477 
3478  // Per feature, per feature-edge a list of attraction points and their
3479  // originating vertex.
3480  List<List<DynamicList<point>>> edgeAttractors(features.size());
3481  List<List<DynamicList<pointConstraint>>> edgeConstraints
3482  (
3483  features.size()
3484  );
3485  forAll(features, feati)
3486  {
3487  label nFeatEdges = features[feati].edges().size();
3488  edgeAttractors[feati].setSize(nFeatEdges);
3489  edgeConstraints[feati].setSize(nFeatEdges);
3490  }
3491 
3492  // Per feature, per feature-point the pp point that is attracted to it.
3493  // This list is only used to subset the feature-points that are actually
3494  // used.
3495  List<labelList> pointAttractor(features.size());
3496  List<List<pointConstraint>> pointConstraints(features.size());
3497  forAll(features, feati)
3498  {
3499  label nFeatPoints = features[feati].points().size();
3500  pointAttractor[feati].setSize(nFeatPoints, -1);
3501  pointConstraints[feati].setSize(nFeatPoints);
3502  }
3503 
3504  // Reverse: from pp point to nearest feature
3505  vectorField rawPatchAttraction(pp.nPoints(), Zero);
3506  List<pointConstraint> rawPatchConstraints(pp.nPoints());
3507 
3508  determineFeatures
3509  (
3510  iter,
3511  featureCos,
3512  multiRegionFeatureSnap,
3513 
3514  pp,
3515  snapDist, // per point max distance and nearest surface
3516  nearestDisp,
3517 
3518  pointFaceSurfNormals, // per face nearest surface
3519  pointFaceDisp,
3520  pointFaceCentres,
3521  pointFacePatchID,
3522 
3523  // Feature-point to pp point
3524  pointAttractor,
3525  pointConstraints,
3526  // Feature-edge to pp point
3527  edgeAttractors,
3528  edgeConstraints,
3529  // pp point to nearest feature
3530  rawPatchAttraction,
3531  rawPatchConstraints
3532  );
3533 
3534  // Print a bit about the attraction from patch point to feature
3535  if (debug)
3536  {
3537  Info<< "Raw geometric feature analysis : ";
3538  writeStats(pp, isPatchMasterPoint, rawPatchConstraints);
3539  }
3540 
3541  // Baffle handling
3542  // ~~~~~~~~~~~~~~~
3543  // Override pointAttractor, edgeAttractor, rawPatchAttration etc. to
3544  // implement 'baffle' handling.
3545  // Baffle: the mesh pp point originates from a loose standing
3546  // baffle.
3547  // Sampling the surface with the surrounding face-centres only picks up
3548  // a single triangle normal so above determineFeatures will not have
3549  // detected anything. So explicitly pick up feature edges on the pp
3550  // (after duplicating points & smoothing so will already have been
3551  // expanded) and match these to the features.
3552  if (detectBaffles)
3553  {
3554  determineBaffleFeatures
3555  (
3556  iter,
3557  baffleFeaturePoints,
3558  featureCos,
3559 
3560  pp,
3561  snapDist,
3562 
3563  // Feature-point to pp point
3564  pointAttractor,
3565  pointConstraints,
3566  // Feature-edge to pp point
3567  edgeAttractors,
3568  edgeConstraints,
3569  // pp point to nearest feature
3570  rawPatchAttraction,
3571  rawPatchConstraints
3572  );
3573  }
3574 
3575  // Print a bit about the attraction from patch point to feature
3576  if (debug)
3577  {
3578  Info<< "After baffle feature analysis : ";
3579  writeStats(pp, isPatchMasterPoint, rawPatchConstraints);
3580  }
3581 
3582 
3583  // Reverse lookup: Find nearest mesh point to feature edge
3584  // ~~~~~~~~~~~~~~~~----------------~~~~~~~~~~~~~~~~~~~~~~~
3585  // go through all edgeAttractors and find the nearest point on pp
3586 
3587  reverseAttractMeshPoints
3588  (
3589  iter,
3590 
3591  pp,
3592  snapDist,
3593 
3594  // Feature-point to pp point
3595  pointAttractor,
3596  pointConstraints,
3597  // Feature-edge to pp point
3598  edgeAttractors,
3599  edgeConstraints,
3600 
3601  // Estimated feature point
3602  rawPatchAttraction,
3603  rawPatchConstraints,
3604 
3605  // pp point to nearest feature
3606  patchAttraction,
3607  patchConstraints
3608  );
3609 
3610  // Print a bit about the attraction from patch point to feature
3611  if (debug)
3612  {
3613  Info<< "Reverse attract feature analysis : ";
3614  writeStats(pp, isPatchMasterPoint, patchConstraints);
3615  }
3616 
3617  // Dump
3619  {
3620  OBJstream featureEdgeStr
3621  (
3622  meshRefiner_.mesh().time().path()
3623  / "edgeAttractors_" + name(iter) + ".obj"
3624  );
3625  Info<< "Dumping feature-edge attraction to "
3626  << featureEdgeStr.name() << endl;
3627 
3628  OBJstream featurePointStr
3629  (
3630  meshRefiner_.mesh().time().path()
3631  / "pointAttractors_" + name(iter) + ".obj"
3632  );
3633  Info<< "Dumping feature-point attraction to "
3634  << featurePointStr.name() << endl;
3635 
3636  forAll(patchConstraints, pointi)
3637  {
3638  const point& pt = pp.localPoints()[pointi];
3639  const vector& attr = patchAttraction[pointi];
3640 
3641  if (patchConstraints[pointi].first() == 2)
3642  {
3643  featureEdgeStr.write(linePointRef(pt, pt+attr));
3644  }
3645  else if (patchConstraints[pointi].first() == 3)
3646  {
3647  featurePointStr.write(linePointRef(pt, pt+attr));
3648  }
3649  }
3650  }
3651 
3652 
3653  //MEJ: any faces that have multi-patch points only keep the
3654  // multi-patch
3655  // points. The other points on the face will be dragged along
3656  // (hopefully)
3657  if (releasePoints)
3658  {
3659  releasePointsNextToMultiPatch
3660  (
3661  iter,
3662  featureCos,
3663 
3664  pp,
3665  snapDist,
3666 
3667  pointFaceCentres,
3668  pointFacePatchID,
3669 
3670  rawPatchAttraction,
3671  rawPatchConstraints,
3672 
3673  patchAttraction,
3674  patchConstraints
3675  );
3676  }
3677 
3678 
3679  // Snap edges to feature edges
3680  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
3681  // Walk existing edges and snap remaining ones (that are marked as
3682  // feature edges in rawPatchConstraints)
3683  if (stringFeatures)
3684  {
3685  stringFeatureEdges
3686  (
3687  iter,
3688  featureCos,
3689 
3690  pp,
3691  snapDist,
3692 
3693  rawPatchAttraction,
3694  rawPatchConstraints,
3695 
3696  patchAttraction,
3697  patchConstraints
3698  );
3699  }
3700 
3701 
3702  // Avoid diagonal attraction
3703  // ~~~~~~~~~~~~~~~~~~~~~~~~~
3704  // Attract one of the non-diagonal points.
3705  if (avoidDiagonal)
3706  {
3707  avoidDiagonalAttraction
3708  (
3709  iter,
3710  featureCos,
3711  pp,
3712  patchAttraction,
3713  patchConstraints
3714  );
3715  }
3716 
3717 
3719  {
3720  dumpMove
3721  (
3722  meshRefiner_.mesh().time().path()
3723  / "patchAttraction_" + name(iter) + ".obj",
3724  pp.localPoints(),
3725  pp.localPoints() + patchAttraction
3726  );
3727  }
3728 }
3729 
3730 
3731 // Correct for squeezing of face
3732 void Foam::snappySnapDriver::preventFaceSqueeze
3733 (
3734  const label iter,
3735  const scalar featureCos,
3736 
3737  const indirectPrimitivePatch& pp,
3738  const scalarField& snapDist,
3739  const vectorField& nearestAttraction,
3740 
3741  vectorField& patchAttraction,
3742  List<pointConstraint>& patchConstraints
3743 ) const
3744 {
3745  autoPtr<OBJstream> strPtr;
3747  {
3748  strPtr.reset
3749  (
3750  new OBJstream
3751  (
3752  meshRefiner_.mesh().time().path()
3753  / "faceSqueeze_" + name(iter) + ".obj"
3754  )
3755  );
3756  Info<< "Dumping faceSqueeze corrections to "
3757  << strPtr().name() << endl;
3758  }
3759 
3761  face singleF;
3762  forAll(pp.localFaces(), facei)
3763  {
3764  const face& f = pp.localFaces()[facei];
3765 
3766  if (f.size() != points.size())
3767  {
3768  points.setSize(f.size());
3769  singleF.setSize(f.size());
3770  for (label i = 0; i < f.size(); i++)
3771  {
3772  singleF[i] = i;
3773  }
3774  }
3775  label nConstraints = 0;
3776  forAll(f, fp)
3777  {
3778  label pointi = f[fp];
3779  const point& pt = pp.localPoints()[pointi];
3780 
3781  if (patchConstraints[pointi].first() > 1)
3782  {
3783  points[fp] = pt + patchAttraction[pointi];
3784  nConstraints++;
3785  }
3786  else
3787  {
3788  points[fp] = pt;
3789  }
3790  }
3791 
3792  if (nConstraints == f.size())
3793  {
3794  if (f.size() == 3)
3795  {
3796  // Triangle: knock out attraction altogether
3797 
3798  // For now keep the points on the longest edge
3799  label maxFp = -1;
3800  scalar maxS = -1;
3801  forAll(f, fp)
3802  {
3803  const point& pt = pp.localPoints()[f[fp]];
3804  const point& nextPt = pp.localPoints()[f.nextLabel(fp)];
3805 
3806  scalar s = magSqr(pt-nextPt);
3807  if (s > maxS)
3808  {
3809  maxS = s;
3810  maxFp = fp;
3811  }
3812  }
3813  if (maxFp != -1)
3814  {
3815  label pointi = f.prevLabel(maxFp);
3816 
3817  // Reset attraction on pointi to nearest
3818 
3819  const point& pt = pp.localPoints()[pointi];
3820 
3821  //Pout<< "** on triangle " << pp.faceCentres()[facei]
3822  // << " knocking out attraction to " << pointi
3823  // << " at:" << pt
3824  // << endl;
3825 
3826  patchAttraction[pointi] = nearestAttraction[pointi];
3827 
3828  if (strPtr)
3829  {
3830  strPtr().write
3831  (
3832  linePointRef(pt, pt+patchAttraction[pointi])
3833  );
3834  }
3835  }
3836  }
3837  else
3838  {
3839  scalar oldArea = f.mag(pp.localPoints());
3840  scalar newArea = singleF.mag(points);
3841  if (newArea < 0.1*oldArea)
3842  {
3843  // For now remove the point with largest distance
3844  label maxFp = -1;
3845  scalar maxS = -1;
3846  forAll(f, fp)
3847  {
3848  scalar s = magSqr(patchAttraction[f[fp]]);
3849  if (s > maxS)
3850  {
3851  maxS = s;
3852  maxFp = fp;
3853  }
3854  }
3855  if (maxFp != -1)
3856  {
3857  label pointi = f[maxFp];
3858  // Lower attraction on pointi
3859  patchAttraction[pointi] *= 0.5;
3860  }
3861  }
3862  }
3863  }
3864  }
3865 }
3866 
3867 
3868 Foam::vectorField Foam::snappySnapDriver::calcNearestSurfaceFeature
3869 (
3870  const snapParameters& snapParams,
3871  const bool alignMeshEdges,
3872  const label iter,
3873  const scalar featureCos,
3874  const scalar featureAttract,
3875  const scalarField& snapDist,
3876  const vectorField& nearestDisp,
3877  const vectorField& nearestNormal,
3878  motionSmoother& meshMover,
3879  vectorField& patchAttraction,
3880  List<pointConstraint>& patchConstraints,
3881 
3882  DynamicList<label>& splitFaces,
3883  DynamicList<labelPair>& splits
3884 
3885 ) const
3886 {
3887  if (dryRun_)
3888  {
3889  return nearestDisp;
3890  }
3891 
3892 
3893  const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
3894  const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
3895  const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
3896 
3897  Info<< "Overriding displacement on features :" << nl
3898  << " implicit features : " << implicitFeatureAttraction << nl
3899  << " explicit features : " << explicitFeatureAttraction << nl
3900  << " multi-patch features : " << multiRegionFeatureSnap << nl
3901  << endl;
3902 
3903  const indirectPrimitivePatch& pp = meshMover.patch();
3904  const pointField& localPoints = pp.localPoints();
3905  const fvMesh& mesh = meshRefiner_.mesh();
3906 
3907 
3908  //const bitSet isMasterPoint(syncTools::getMasterPoints(mesh));
3909  const bitSet isPatchMasterPoint
3910  (
3912  (
3913  mesh,
3914  pp.meshPoints()
3915  )
3916  );
3917 
3918  // Per point, per surrounding face:
3919  // - faceSurfaceNormal
3920  // - faceDisp
3921  // - faceCentres
3922  List<List<point>> pointFaceSurfNormals;
3923  List<List<point>> pointFaceDisp;
3924  List<List<point>> pointFaceCentres;
3925  List<labelList> pointFacePatchID;
3926 
3927  {
3928  // Calculate attraction distance per face (from the attraction distance
3929  // per point)
3930  scalarField faceSnapDist(pp.size(), -GREAT);
3931  forAll(pp.localFaces(), facei)
3932  {
3933  const face& f = pp.localFaces()[facei];
3934  forAll(f, fp)
3935  {
3936  faceSnapDist[facei] = max(faceSnapDist[facei], snapDist[f[fp]]);
3937  }
3938  }
3939 
3940 
3941  // Displacement and orientation per pp face
3942  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3943 
3944  // vector from point on surface back to face centre
3945  vectorField faceDisp(pp.size(), Zero);
3946  // normal of surface at point on surface
3947  vectorField faceSurfaceNormal(pp.size(), Zero);
3948  labelList faceSurfaceGlobalRegion(pp.size(), -1);
3949  //vectorField faceRotation(pp.size(), Zero);
3950 
3951  calcNearestFace
3952  (
3953  iter,
3954  pp,
3955  faceSnapDist,
3956  faceDisp,
3957  faceSurfaceNormal,
3958  faceSurfaceGlobalRegion
3959  //faceRotation
3960  );
3961 
3962 
3963  // Collect (possibly remote) per point data of all surrounding faces
3964  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3965  // - faceSurfaceNormal
3966  // - faceDisp
3967  // - faceCentres
3968  calcNearestFacePointProperties
3969  (
3970  iter,
3971  pp,
3972 
3973  faceDisp,
3974  faceSurfaceNormal,
3975  faceSurfaceGlobalRegion,
3976 
3977  pointFaceSurfNormals,
3978  pointFaceDisp,
3979  pointFaceCentres,
3980  pointFacePatchID
3981  );
3982  }
3983 
3984 
3985  // Start off with nearest point on surface
3986  vectorField patchDisp = nearestDisp;
3987 
3988 
3989  // Main calculation
3990  // ~~~~~~~~~~~~~~~~
3991  // This is the main intelligence which calculates per point the vector to
3992  // attract it to the nearest surface. There are lots of possibilities
3993  // here.
3994 
3995  // Nearest feature
3996  patchAttraction.setSize(localPoints.size());
3997  patchAttraction = Zero;
3998  // Constraints at feature
3999  patchConstraints.setSize(localPoints.size());
4000  patchConstraints = pointConstraint();
4001 
4002  if (implicitFeatureAttraction)
4003  {
4004  // Sample faces around each point and see if nearest surface normal
4005  // differs. Reconstruct a feature edge/point if possible and snap to
4006  // it.
4007  featureAttractionUsingReconstruction
4008  (
4009  iter,
4010  featureCos,
4011 
4012  pp,
4013  snapDist,
4014  nearestDisp,
4015 
4016  pointFaceSurfNormals,
4017  pointFaceDisp,
4018  pointFaceCentres,
4019  pointFacePatchID,
4020 
4021  patchAttraction,
4022  patchConstraints
4023  );
4024  }
4025 
4026  if (explicitFeatureAttraction)
4027  {
4028  // Only do fancy stuff if alignMeshEdges
4029  bool releasePoints = false;
4030  bool stringFeatures = false;
4031  bool avoidDiagonal = false;
4032  if (alignMeshEdges)
4033  {
4034  releasePoints = snapParams.releasePoints();
4035  stringFeatures = snapParams.stringFeatures();
4036  avoidDiagonal = snapParams.avoidDiagonal();
4037  }
4038 
4039 
4040  // Sample faces around each point and see if nearest surface normal
4041  // differs. For those find the nearest real feature edge/point and
4042  // store the correspondence. Then loop over feature edge/point
4043  // and attract those nearest mesh point. (the first phase just is
4044  // a subsetting of candidate points, the second makes sure that only
4045  // one mesh point gets attracted per feature)
4046  featureAttractionUsingFeatureEdges
4047  (
4048  iter,
4049  multiRegionFeatureSnap,
4050 
4051  snapParams.detectBaffles(),
4052  snapParams.baffleFeaturePoints(), // all points on baffle edges
4053  // are attracted to feature pts
4054 
4055  releasePoints,
4056  stringFeatures,
4057  avoidDiagonal,
4058 
4059  featureCos,
4060 
4061  pp,
4062  snapDist,
4063  nearestDisp,
4064  nearestNormal,
4065 
4066  pointFaceSurfNormals,
4067  pointFaceDisp,
4068  pointFaceCentres,
4069  pointFacePatchID,
4070 
4071  patchAttraction,
4072  patchConstraints
4073  );
4074  }
4075 
4076  if (!alignMeshEdges)
4077  {
4078  const scalar concaveCos = Foam::cos
4079  (
4080  degToRad(snapParams.concaveAngle())
4081  );
4082  const scalar minAreaRatio = snapParams.minAreaRatio();
4083 
4084  Info<< "Experimental: introducing face splits to avoid rotating"
4085  << " mesh edges. Splitting faces when" << nl
4086  << indent << "- angle not concave by more than "
4087  << snapParams.concaveAngle() << " degrees" << nl
4088  << indent << "- resulting triangles of similar area "
4089  << " (ratio within " << minAreaRatio << ")" << nl
4090  << endl;
4091 
4092  splitDiagonals
4093  (
4094  featureCos,
4095  concaveCos,
4096  minAreaRatio,
4097  pp,
4098 
4099  nearestDisp,
4100  nearestNormal,
4101 
4102  patchAttraction,
4103  patchConstraints,
4104  splitFaces,
4105  splits
4106  );
4107 
4108  if (debug)
4109  {
4110  Info<< "Diagonal attraction feature correction : ";
4111  writeStats(pp, isPatchMasterPoint, patchConstraints);
4112  }
4113  }
4114 
4115 
4116  preventFaceSqueeze
4117  (
4118  iter,
4119  featureCos,
4120 
4121  pp,
4122  snapDist,
4123  nearestDisp,
4124 
4125  patchAttraction,
4126  patchConstraints
4127  );
4128 
4129  {
4130  vector avgPatchDisp = meshRefinement::gAverage
4131  (
4132  isPatchMasterPoint,
4133  patchDisp
4134  );
4135  vector avgPatchAttr = meshRefinement::gAverage
4136  (
4137  isPatchMasterPoint,
4138  patchAttraction
4139  );
4140 
4141  Info<< "Attraction:" << endl
4142  << " linear : max:" << gMaxMagSqr(patchDisp)
4143  << " avg:" << avgPatchDisp << endl
4144  << " feature : max:" << gMaxMagSqr(patchAttraction)
4145  << " avg:" << avgPatchAttr << endl;
4146  }
4147 
4148  // So now we have:
4149  // - patchDisp : point movement to go to nearest point on surface
4150  // (either direct or through interpolation of
4151  // face nearest)
4152  // - patchAttraction : direct attraction to features
4153  // - patchConstraints : type of features
4154 
4155  // Use any combination of patchDisp and direct feature attraction.
4156 
4157 
4158  // Mix with direct feature attraction
4159  forAll(patchConstraints, pointi)
4160  {
4161  if (patchConstraints[pointi].first() > 1)
4162  {
4163  patchDisp[pointi] =
4164  (1.0-featureAttract)*patchDisp[pointi]
4165  + featureAttract*patchAttraction[pointi];
4166  }
4167  }
4168 
4169 
4170 
4171  // Count
4172  {
4173  Info<< "Feature analysis : ";
4174  writeStats(pp, isPatchMasterPoint, patchConstraints);
4175  }
4176 
4177 
4178  // Now we have the displacement per patch point to move onto the surface
4179  // Split into tangential and normal direction.
4180  // - start off with all non-constrained points following the constrained
4181  // ones since point normals not relevant.
4182  // - finish with only tangential component smoothed.
4183  // Note: tangential is most
4184  // likely to come purely from face-centre snapping, not face rotation.
4185  // Note: could use the constraints here (constraintTransformation())
4186  // but this is not necessarily accurate and we're smoothing to
4187  // get out of problems.
4188 
4189  if (featureAttract < 1-0.001)
4190  {
4191  //const bitSet isMasterEdge(syncTools::getMasterEdges(mesh));
4192  const labelList meshEdges
4193  (
4194  pp.meshEdges(mesh.edges(), mesh.pointEdges())
4195  );
4196  const bitSet isPatchMasterEdge
4197  (
4199  (
4200  mesh,
4201  meshEdges
4202  )
4203  );
4204 
4205  const vectorField pointNormals
4206  (
4208  (
4209  mesh,
4210  pp
4211  )
4212  );
4213 
4214  // 1. Smoothed all displacement
4215  vectorField smoothedPatchDisp = patchDisp;
4216  smoothAndConstrain
4217  (
4218  isPatchMasterEdge,
4219  pp,
4220  meshEdges,
4221  patchConstraints,
4222  smoothedPatchDisp
4223  );
4224 
4225 
4226  // 2. Smoothed tangential component
4227  vectorField tangPatchDisp = patchDisp;
4228  tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
4229  smoothAndConstrain
4230  (
4231  isPatchMasterEdge,
4232  pp,
4233  meshEdges,
4234  patchConstraints,
4235  tangPatchDisp
4236  );
4237 
4238  // Re-add normal component
4239  tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
4240 
4242  {
4243  dumpMove
4244  (
4245  mesh.time().path()
4246  / "tangPatchDispConstrained_" + name(iter) + ".obj",
4247  pp.localPoints(),
4248  pp.localPoints() + tangPatchDisp
4249  );
4250  }
4251 
4252  patchDisp =
4253  (1.0-featureAttract)*smoothedPatchDisp
4254  + featureAttract*tangPatchDisp;
4255  }
4256 
4257 
4258  const scalar relax = featureAttract;
4259  patchDisp *= relax;
4260 
4261 
4262  // Points on zones in one domain but only present as point on other
4263  // will not do condition 2 on all. Sync explicitly.
4265  (
4266  mesh,
4267  pp.meshPoints(),
4268  patchDisp,
4269  minMagSqrEqOp<point>(), // combine op
4270  vector(GREAT, GREAT, GREAT) // null value (note: cannot use VGREAT)
4271  );
4272 
4273  return patchDisp;
4274 }
4275 
4276 
4277 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::syncTools::syncEdgeList
static void syncEdgeList(const polyMesh &mesh, List< T > &edgeValues, const CombineOp &cop, const T &nullValue, const TransformOp &top, const FlipOp &fop)
Synchronize values on all mesh edges.
Definition: syncToolsTemplates.C:800
Foam::polyMesh::points
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1069
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::PatchTools::pointNormals
static tmp< pointField > pointNormals(const polyMesh &, const PrimitivePatch< FaceList, PointField > &)
Return parallel consistent point normals for patches using mesh points.
Foam::meshRefinement::gAverage
static T gAverage(const bitSet &isMasterElem, const UList< T > &values)
Helper: calculate average.
Definition: meshRefinementTemplates.C:62
Foam::edgeList
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
snappySnapDriver.H
Foam::labelMax
constexpr label labelMax
Definition: label.H:61
Foam::VectorSpace< Vector< Cmpt >, Cmpt, 3 >::uniform
static Vector< Cmpt > uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Definition: VectorSpaceI.H:164
Foam::returnReduce
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Definition: PstreamReduceOps.H:94
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Definition: gmvOutputSpray.H:25
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
PatchTools.H
Foam::diag
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
Definition: pointPatchFieldFunctions.H:287
motionSmoother.H
polyTopoChange.H
indexedOctree.H
localPointRegion.H
Foam::meshRefinement::ATTRACTION
Definition: meshRefinement.H:97
unitConversion.H
Unit conversion functions.
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:65
Foam::surfaceZonesInfo::getNamedSurfaces
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
Definition: surfaceZonesInfo.C:263
Foam::polyMesh::boundaryMesh
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:444
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::primitiveMesh::edges
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
Definition: primitiveMeshEdges.C:505
Foam::syncTools::getMasterFaces
static bitSet getMasterFaces(const polyMesh &mesh)
Definition: syncTools.C:126
refinementFeatures.H
Foam::primitiveMesh::pointEdges
const labelListList & pointEdges() const
Definition: primitiveMeshEdges.C:516
Foam::transform
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
Definition: dimensionSet.C:521
syncTools.H
Foam::syncTools::syncPointList
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
Definition: syncToolsTemplates.C:721
Foam::primitiveMesh::nPoints
label nPoints() const noexcept
Number of mesh points.
Definition: primitiveMeshI.H:37
Foam::gMaxMagSqr
Type gMaxMagSqr(const UList< Type > &f, const label comm)
Definition: FieldFunctions.C:590
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:54
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::meshRefinement::getMasterEdges
static bitSet getMasterEdges(const polyMesh &mesh, const labelList &meshEdges)
Determine master edge for subset of edges. If coupled.
Definition: meshRefinement.C:3140
Foam::labelPair
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:54
Foam::wordList
List< word > wordList
A List of words.
Definition: fileName.H:62
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:51
Foam::PointIndexHit
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:52
Foam::constant::physicoChemical::c1
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
Foam::Field< vector >
treeDataPoint.H
plane.H
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
refinementSurfaces.H
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:222
patchID
label patchID
Definition: boundaryProcessorFaPatchPoints.H:5
Foam::polyMesh::faceZones
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
Definition: polyMesh.H:486
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::edgeMeshTools::writeStats
void writeStats(Ostream &os, const extendedFeatureEdgeMesh &emesh)
Write some information.
Definition: edgeMeshTools.C:36
Foam::FatalError
error FatalError
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:121
Foam::localPointRegion::findDuplicateFaces
static labelList findDuplicateFaces(const primitiveMesh &, const labelList &)
Helper routine to find baffles (two boundary faces using the.
Definition: localPointRegion.C:544
fvMesh.H
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::ZoneMesh::findZoneID
label findZoneID(const word &zoneName) const
Find zone index by name, return -1 if not found.
Definition: ZoneMesh.C:519
Foam::indent
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:339
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::degToRad
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Definition: unitConversion.H:48
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::surfaceZonesInfo::getUnnamedSurfaces
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
Definition: surfaceZonesInfo.C:242
Foam::indirectPrimitivePatch
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field.
Definition: indirectPrimitivePatch.H:49
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
points0
pointField points0(pointIOField(IOobject("points", mesh.time().constant(), polyMesh::meshSubDir, mesh, IOobject::MUST_READ, IOobject::NO_WRITE, false)))
Foam::linePointRef
line< point, const point & > linePointRef
A line using referred points.
Definition: linePointRef.H:47
Foam::polyMesh::faces
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1094
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::pointIndexHit
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::Pair< label >
Foam::Time::path
fileName path() const
Return path.
Definition: Time.H:358
f
labelList f(nPoints)
Foam::foamVersion::patch
const std::string patch
OpenFOAM patch number as a std::string.
Foam::primitiveMesh::nInternalFaces
label nInternalFaces() const noexcept
Number of internal faces.
Definition: primitiveMeshI.H:78
Foam::List
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: BitOps.H:63
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::primitiveMesh::faceCentres
const vectorField & faceCentres() const
Definition: primitiveMeshFaceCentresAndAreas.C:77
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::identity
labelList identity(const label len, label start=0)
Create identity map of the given length with (map[i] == i)
Definition: labelList.C:38
Foam::List::clear
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
x
x
Definition: LISASMDCalcMethod2.H:52
snapParameters.H
Foam::List::set
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:341
split
static bool split(const std::string &line, std::string &key, std::string &val)
Definition: cpuInfo.C:39
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
rndGen
Random rndGen
Definition: createFields.H:23
Foam::fvMesh::time
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:280
Foam::VectorSpace< Vector< scalar >, scalar, 3 >::zero
static const Vector< scalar > zero
Definition: VectorSpace.H:115
Foam::primitiveMesh::nFaces
label nFaces() const noexcept
Number of mesh faces.
Definition: primitiveMeshI.H:90
Foam::listPlusEqOp::operator()
void operator()(List< T > &x, const List< T > &y) const
Definition: snappySnapDriverFeature.C:56
Foam::sortedOrder
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
p0
const volScalarField & p0
Definition: EEqn.H:36
featureEdgeMesh.H
Foam::VectorSpace::size
static constexpr direction size() noexcept
The number of elements in the VectorSpace = Ncmpts.
Definition: VectorSpace.H:176
Foam::Tuple2
A 2-tuple for storing two objects of dissimilar types. The container is similar in purpose to std::pa...
Definition: stringOps.H:60
Foam::point
vector point
Point is a vector.
Definition: point.H:43
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
pyramidPointFaceRef.H
pFaces
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
Foam::listPlusEqOp
Definition: snappySnapDriverFeature.C:52
OBJstream.H
Foam::meshRefinement::getMasterPoints
static bitSet getMasterPoints(const polyMesh &mesh, const labelList &meshPoints)
Determine master point for subset of points. If coupled.
Definition: meshRefinement.C:3104
y
scalar y
Definition: LISASMDCalcMethod1.H:14
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:265
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: UIndirectList.H:58
relax
UEqn relax()
Foam::primitiveMesh::faceAreas
const vectorField & faceAreas() const
Definition: primitiveMeshFaceCentresAndAreas.C:89