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 | Copyright (C) 2011-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "snappySnapDriver.H"
27 #include "polyTopoChange.H"
28 #include "syncTools.H"
29 #include "fvMesh.H"
30 #include "OBJstream.H"
31 #include "motionSmoother.H"
32 #include "refinementSurfaces.H"
33 #include "refinementFeatures.H"
34 #include "unitConversion.H"
35 #include "plane.H"
36 #include "featureEdgeMesh.H"
37 #include "treeDataPoint.H"
38 #include "indexedOctree.H"
39 #include "snapParameters.H"
40 #include "PatchTools.H"
41 
42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43 
44 namespace Foam
45 {
46  template<class T>
48  {
49  public:
50 
51  void operator()(List<T>& x, const List<T>& y) const
52  {
53  label sz = x.size();
54  x.setSize(sz+y.size());
55  forAll(y, i)
56  {
57  x[sz++] = y[i];
58  }
59  }
60  };
61 }
62 
63 
64 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
65 
66 bool Foam::snappySnapDriver::isFeaturePoint
67 (
68  const scalar featureCos,
69  const indirectPrimitivePatch& pp,
70  const PackedBoolList& isFeatureEdge,
71  const label pointi
72 ) const
73 {
74  const pointField& points = pp.localPoints();
75  const edgeList& edges = pp.edges();
76  const labelList& pEdges = pp.pointEdges()[pointi];
77 
78  label nFeatEdges = 0;
79 
80  forAll(pEdges, i)
81  {
82  if (isFeatureEdge[pEdges[i]])
83  {
84  nFeatEdges++;
85 
86  for (label j = i+1; j < pEdges.size(); j++)
87  {
88  if (isFeatureEdge[pEdges[j]])
89  {
90  const edge& eI = edges[pEdges[i]];
91  const edge& eJ = edges[pEdges[j]];
92 
93  const point& p = points[pointi];
94  const point& pI = points[eI.otherVertex(pointi)];
95  const point& pJ = points[eJ.otherVertex(pointi)];
96 
97  vector vI = p-pI;
98  scalar vIMag = mag(vI);
99 
100  vector vJ = pJ-p;
101  scalar vJMag = mag(vJ);
102 
103  if
104  (
105  vIMag > SMALL
106  && vJMag > SMALL
107  && ((vI/vIMag & vJ/vJMag) < featureCos)
108  )
109  {
110  return true;
111  }
112  }
113  }
114  }
115  }
116 
117  if (nFeatEdges == 1)
118  {
119  // End of feature-edge string
120  return true;
121  }
122 
123  return false;
124 }
125 
126 
127 void Foam::snappySnapDriver::smoothAndConstrain
128 (
129  const PackedBoolList& isPatchMasterEdge,
130  const indirectPrimitivePatch& pp,
131  const labelList& meshEdges,
132  const List<pointConstraint>& constraints,
133  vectorField& disp
134 ) const
135 {
136  const fvMesh& mesh = meshRefiner_.mesh();
137 
138  for (label avgIter = 0; avgIter < 20; avgIter++)
139  {
140  // Calculate average displacement of neighbours
141  // - unconstrained (i.e. surface) points use average of all
142  // neighbouring points
143  // - from testing it has been observed that it is not beneficial
144  // to have edge constrained points use average of all edge or point
145  // constrained neighbours since they're already attracted to
146  // the nearest point on the feature.
147  // Having them attract to point-constrained neighbours does not
148  // make sense either since there is usually just one of them so
149  // it severely distorts it.
150  // - same for feature points. They are already attracted to the
151  // nearest feature point.
152 
153  vectorField dispSum(pp.nPoints(), Zero);
154  labelList dispCount(pp.nPoints(), 0);
155 
156  const labelListList& pointEdges = pp.pointEdges();
157  const edgeList& edges = pp.edges();
158 
159  forAll(pointEdges, pointi)
160  {
161  const labelList& pEdges = pointEdges[pointi];
162 
163  label nConstraints = constraints[pointi].first();
164 
165  if (nConstraints <= 1)
166  {
167  forAll(pEdges, i)
168  {
169  label edgeI = pEdges[i];
170 
171  if (isPatchMasterEdge[edgeI])
172  {
173  label nbrPointi = edges[edgeI].otherVertex(pointi);
174  if (constraints[nbrPointi].first() >= nConstraints)
175  {
176  dispSum[pointi] += disp[nbrPointi];
177  dispCount[pointi]++;
178  }
179  }
180  }
181  }
182  }
183 
185  (
186  mesh,
187  pp.meshPoints(),
188  dispSum,
189  plusEqOp<point>(),
190  vector::zero,
192  );
194  (
195  mesh,
196  pp.meshPoints(),
197  dispCount,
198  plusEqOp<label>(),
199  label(0),
201  );
202 
203  // Constraints
204  forAll(constraints, pointi)
205  {
206  if (dispCount[pointi] > 0)
207  {
208  // Mix my displacement with neighbours' displacement
209  disp[pointi] =
210  0.5
211  *(disp[pointi] + dispSum[pointi]/dispCount[pointi]);
212  }
213  }
214  }
215 }
216 
217 
218 void Foam::snappySnapDriver::calcNearestFace
219 (
220  const label iter,
221  const indirectPrimitivePatch& pp,
222  const scalarField& faceSnapDist,
223  vectorField& faceDisp,
224  vectorField& faceSurfaceNormal,
225  labelList& faceSurfaceGlobalRegion,
226  vectorField& faceRotation
227 ) const
228 {
229  const fvMesh& mesh = meshRefiner_.mesh();
230  const refinementSurfaces& surfaces = meshRefiner_.surfaces();
231 
232  // Displacement and orientation per pp face.
233  faceDisp.setSize(pp.size());
234  faceDisp = Zero;
235  faceSurfaceNormal.setSize(pp.size());
236  faceSurfaceNormal = Zero;
237  faceSurfaceGlobalRegion.setSize(pp.size());
238  faceSurfaceGlobalRegion = -1;
239 
240  // Divide surfaces into zoned and unzoned
241  const labelList zonedSurfaces =
243  const labelList unzonedSurfaces =
245 
246  // Per pp face the current surface snapped to
247  labelList snapSurf(pp.size(), -1);
248 
249 
250  // Do zoned surfaces
251  // ~~~~~~~~~~~~~~~~~
252  // Zoned faces only attract to corresponding surface
253 
254  // Extract faces per zone
255  const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
256 
257  forAll(zonedSurfaces, i)
258  {
259  label zoneSurfI = zonedSurfaces[i];
260 
261  const word& faceZoneName = surfZones[zoneSurfI].faceZoneName();
262 
263  // Get indices of faces on pp that are also in zone
264  label zoneI = mesh.faceZones().findZoneID(faceZoneName);
265  if (zoneI == -1)
266  {
268  << "Problem. Cannot find zone " << faceZoneName
269  << exit(FatalError);
270  }
271  const faceZone& fZone = mesh.faceZones()[zoneI];
272  PackedBoolList isZonedFace(mesh.nFaces());
273  forAll(fZone, i)
274  {
275  isZonedFace[fZone[i]] = 1;
276  }
277 
278  DynamicList<label> ppFaces(fZone.size());
279  DynamicList<label> meshFaces(fZone.size());
280  forAll(pp.addressing(), i)
281  {
282  if (isZonedFace[pp.addressing()[i]])
283  {
284  snapSurf[i] = zoneSurfI;
285  ppFaces.append(i);
286  meshFaces.append(pp.addressing()[i]);
287  }
288  }
289 
290  //Pout<< "For faceZone " << fZone.name()
291  // << " found " << ppFaces.size() << " out of " << pp.size()
292  // << endl;
293 
294  pointField fc
295  (
297  (
298  IndirectList<face>(mesh.faces(), meshFaces),
299  mesh.points()
300  ).faceCentres()
301  );
302 
303  List<pointIndexHit> hitInfo;
304  labelList hitSurface;
305  labelList hitRegion;
306  vectorField hitNormal;
307  surfaces.findNearestRegion
308  (
309  labelList(1, zoneSurfI),
310  fc,
311  sqr(faceSnapDist),// sqr of attract dist
312  hitSurface,
313  hitInfo,
314  hitRegion,
315  hitNormal
316  );
317 
318  forAll(hitInfo, hitI)
319  {
320  if (hitInfo[hitI].hit())
321  {
322  label facei = ppFaces[hitI];
323  faceDisp[facei] = hitInfo[hitI].hitPoint() - fc[hitI];
324  faceSurfaceNormal[facei] = hitNormal[hitI];
325  faceSurfaceGlobalRegion[facei] = surfaces.globalRegion
326  (
327  hitSurface[hitI],
328  hitRegion[hitI]
329  );
330  }
331  else
332  {
334  << "Did not find surface near face centre " << fc[hitI]
335  << endl;
336  }
337  }
338  }
339 
340 
341  // Do unzoned surfaces
342  // ~~~~~~~~~~~~~~~~~~~
343  // Unzoned faces attract to any unzoned surface
344 
345  DynamicList<label> ppFaces(pp.size());
346  DynamicList<label> meshFaces(pp.size());
347  forAll(pp.addressing(), i)
348  {
349  if (snapSurf[i] == -1)
350  {
351  ppFaces.append(i);
352  meshFaces.append(pp.addressing()[i]);
353  }
354  }
355  //Pout<< "Found " << ppFaces.size() << " unzoned faces out of "
356  // << pp.size() << endl;
357 
358  pointField fc
359  (
361  (
362  IndirectList<face>(mesh.faces(), meshFaces),
363  mesh.points()
364  ).faceCentres()
365  );
366 
367  List<pointIndexHit> hitInfo;
368  labelList hitSurface;
369  labelList hitRegion;
370  vectorField hitNormal;
371  surfaces.findNearestRegion
372  (
373  unzonedSurfaces,
374  fc,
375  sqr(faceSnapDist),// sqr of attract dist
376  hitSurface,
377  hitInfo,
378  hitRegion,
379  hitNormal
380  );
381 
382  forAll(hitInfo, hitI)
383  {
384  if (hitInfo[hitI].hit())
385  {
386  label facei = ppFaces[hitI];
387  faceDisp[facei] = hitInfo[hitI].hitPoint() - fc[hitI];
388  faceSurfaceNormal[facei] = hitNormal[hitI];
389  faceSurfaceGlobalRegion[facei] = surfaces.globalRegion
390  (
391  hitSurface[hitI],
392  hitRegion[hitI]
393  );
394  }
395  else
396  {
398  << "Did not find surface near face centre " << fc[hitI]
399  << endl;
400  }
401  }
402 
403 
404  // Determine rotation
405  // ~~~~~~~~~~~~~~~~~~
406 
407  // Determine rotation axis
408  faceRotation.setSize(pp.size());
409  faceRotation = Zero;
410 
411  forAll(faceRotation, facei)
412  {
413  // Note: extend to >180 degrees checking
414  faceRotation[facei] =
415  pp.faceNormals()[facei]
416  ^ faceSurfaceNormal[facei];
417  }
418 
419  if (debug&meshRefinement::ATTRACTION)
420  {
421  dumpMove
422  (
423  mesh.time().path()
424  / "faceDisp_" + name(iter) + ".obj",
425  pp.faceCentres(),
426  pp.faceCentres() + faceDisp
427  );
428  dumpMove
429  (
430  mesh.time().path()
431  / "faceRotation_" + name(iter) + ".obj",
432  pp.faceCentres(),
433  pp.faceCentres() + faceRotation
434  );
435  }
436 }
437 
438 
439 void Foam::snappySnapDriver::calcNearestFacePointProperties
440 (
441  const label iter,
442  const indirectPrimitivePatch& pp,
443 
444  const vectorField& faceDisp,
445  const vectorField& faceSurfaceNormal,
446  const labelList& faceSurfaceGlobalRegion,
447 
448  List<List<point>>& pointFaceSurfNormals,
449  List<List<point>>& pointFaceDisp,
450  List<List<point>>& pointFaceCentres,
451  List<labelList>& pointFacePatchID
452 ) const
453 {
454  const fvMesh& mesh = meshRefiner_.mesh();
455 
456  const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
457 
458 
459  // For now just get all surrounding face data. Expensive - should just
460  // store and sync data on coupled points only
461  // (see e.g PatchToolsNormals.C)
462 
463  pointFaceSurfNormals.setSize(pp.nPoints());
464  pointFaceDisp.setSize(pp.nPoints());
465  pointFaceCentres.setSize(pp.nPoints());
466  pointFacePatchID.setSize(pp.nPoints());
467 
468  // Fill local data
469  forAll(pp.pointFaces(), pointi)
470  {
471  const labelList& pFaces = pp.pointFaces()[pointi];
472 
473  // Count valid face normals
474  label nFaces = 0;
475  forAll(pFaces, i)
476  {
477  label facei = pFaces[i];
478  if (isMasterFace[facei] && faceSurfaceGlobalRegion[facei] != -1)
479  {
480  nFaces++;
481  }
482  }
483 
484 
485  List<point>& pNormals = pointFaceSurfNormals[pointi];
486  pNormals.setSize(nFaces);
487  List<point>& pDisp = pointFaceDisp[pointi];
488  pDisp.setSize(nFaces);
489  List<point>& pFc = pointFaceCentres[pointi];
490  pFc.setSize(nFaces);
491  labelList& pFid = pointFacePatchID[pointi];
492  pFid.setSize(nFaces);
493 
494  nFaces = 0;
495  forAll(pFaces, i)
496  {
497  label facei = pFaces[i];
498  label globalRegionI = faceSurfaceGlobalRegion[facei];
499 
500  if (isMasterFace[facei] && globalRegionI != -1)
501  {
502  pNormals[nFaces] = faceSurfaceNormal[facei];
503  pDisp[nFaces] = faceDisp[facei];
504  pFc[nFaces] = pp.faceCentres()[facei];
505  pFid[nFaces] = globalToMasterPatch_[globalRegionI];
506  nFaces++;
507  }
508  }
509  }
510 
511 
512  // Collect additionally 'normal' boundary faces for boundaryPoints of pp
513  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
514  // points on the boundary of pp should pick up non-pp normals
515  // as well for the feature-reconstruction to behave correctly.
516  // (the movement is already constrained outside correctly so it
517  // is only that the unconstrained attraction vector is calculated
518  // correctly)
519  {
520  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
521  labelList patchID(pbm.patchID());
522 
523  // Unmark all non-coupled boundary faces
524  forAll(pbm, patchi)
525  {
526  const polyPatch& pp = pbm[patchi];
527 
528  if (pp.coupled() || isA<emptyPolyPatch>(pp))
529  {
530  forAll(pp, i)
531  {
532  label meshFacei = pp.start()+i;
533  patchID[meshFacei-mesh.nInternalFaces()] = -1;
534  }
535  }
536  }
537 
538  // Remove any meshed faces
539  forAll(pp.addressing(), i)
540  {
541  label meshFacei = pp.addressing()[i];
542  patchID[meshFacei-mesh.nInternalFaces()] = -1;
543  }
544 
545  // See if pp point uses any non-meshed boundary faces
546 
547  const labelList& boundaryPoints = pp.boundaryPoints();
548  forAll(boundaryPoints, i)
549  {
550  label pointi = boundaryPoints[i];
551  label meshPointi = pp.meshPoints()[pointi];
552  const point& pt = mesh.points()[meshPointi];
553  const labelList& pFaces = mesh.pointFaces()[meshPointi];
554 
555  List<point>& pNormals = pointFaceSurfNormals[pointi];
556  List<point>& pDisp = pointFaceDisp[pointi];
557  List<point>& pFc = pointFaceCentres[pointi];
558  labelList& pFid = pointFacePatchID[pointi];
559 
560  forAll(pFaces, i)
561  {
562  label meshFacei = pFaces[i];
563  if (!mesh.isInternalFace(meshFacei))
564  {
565  label patchi = patchID[meshFacei-mesh.nInternalFaces()];
566 
567  if (patchi != -1)
568  {
569  vector fn = mesh.faceAreas()[meshFacei];
570  pNormals.append(fn/mag(fn));
571  pDisp.append(mesh.faceCentres()[meshFacei]-pt);
572  pFc.append(mesh.faceCentres()[meshFacei]);
573  pFid.append(patchi);
574  }
575  }
576  }
577  }
578  }
579 
581  (
582  mesh,
583  pp.meshPoints(),
584  pointFaceSurfNormals,
586  List<point>(),
588  );
590  (
591  mesh,
592  pp.meshPoints(),
593  pointFaceDisp,
595  List<point>(),
597  );
599  (
600  mesh,
601  pp.meshPoints(),
602  pointFaceCentres,
604  List<point>(),
606  );
608  (
609  mesh,
610  pp.meshPoints(),
611  pointFacePatchID,
613  List<label>()
614  );
615 
616 
617  // Sort the data according to the face centres. This is only so we get
618  // consistent behaviour serial and parallel.
619  labelList visitOrder;
620  forAll(pointFaceDisp, pointi)
621  {
622  List<point>& pNormals = pointFaceSurfNormals[pointi];
623  List<point>& pDisp = pointFaceDisp[pointi];
624  List<point>& pFc = pointFaceCentres[pointi];
625  labelList& pFid = pointFacePatchID[pointi];
626 
627  sortedOrder(mag(pFc)(), visitOrder);
628 
629  pNormals = List<point>(pNormals, visitOrder);
630  pDisp = List<point>(pDisp, visitOrder);
631  pFc = List<point>(pFc, visitOrder);
632  pFid = UIndirectList<label>(pFid, visitOrder)();
633  }
634 }
635 
636 
637 void Foam::snappySnapDriver::correctAttraction
638 (
639  const DynamicList<point>& surfacePoints,
640  const DynamicList<label>& surfaceCounts,
641  const point& edgePt,
642  const vector& edgeNormal, // normalised normal
643  const point& pt,
644 
645  vector& edgeOffset // offset from pt to point on edge
646 ) const
647 {
648  // Tangential component along edge
649  scalar tang = ((pt-edgePt)&edgeNormal);
650 
651  labelList order;
652  Foam::sortedOrder(surfaceCounts, order);
653 
654  if (order[0] < order[1])
655  {
656  // There is a non-dominant plane. Use the point on the plane to
657  // attract to.
658  vector attractD = surfacePoints[order[0]]-edgePt;
659  // Tangential component along edge
660  scalar tang2 = (attractD&edgeNormal);
661  // Normal component
662  attractD -= tang2*edgeNormal;
663  // Calculate fraction of normal distances
664  scalar magAttractD = mag(attractD);
665  scalar fraction = magAttractD/(magAttractD+mag(edgeOffset));
666 
667  point linePt =
668  edgePt
669  + ((1.0-fraction)*tang2 + fraction*tang)*edgeNormal;
670  edgeOffset = linePt-pt;
671  }
672 }
673 
674 
675 Foam::pointIndexHit Foam::snappySnapDriver::findMultiPatchPoint
676 (
677  const point& pt,
678  const labelList& patchIDs,
679  const List<point>& faceCentres
680 ) const
681 {
682  // Determine if multiple patchIDs
683  if (patchIDs.size())
684  {
685  label patch0 = patchIDs[0];
686 
687  for (label i = 1; i < patchIDs.size(); i++)
688  {
689  if (patchIDs[i] != patch0)
690  {
691  return pointIndexHit(true, pt, labelMax);
692  }
693  }
694  }
695  return pointIndexHit(false, Zero, labelMax);
696 }
697 
698 
699 Foam::label Foam::snappySnapDriver::findNormal
700 (
701  const scalar featureCos,
702  const vector& n,
703  const DynamicList<vector>& surfaceNormals
704 ) const
705 {
706  label index = -1;
707 
708  forAll(surfaceNormals, j)
709  {
710  scalar cosAngle = (n&surfaceNormals[j]);
711 
712  if
713  (
714  (cosAngle >= featureCos)
715  || (cosAngle < (-1+0.001)) // triangle baffles
716  )
717  {
718  index = j;
719  break;
720  }
721  }
722  return index;
723 }
724 
725 
726 Foam::pointIndexHit Foam::snappySnapDriver::findMultiPatchPoint
727 (
728  const point& pt,
729  const labelList& patchIDs,
730  const DynamicList<vector>& surfaceNormals,
731  const labelList& faceToNormalBin
732 ) const
733 {
734  if (patchIDs.empty())
735  {
736  return pointIndexHit(false, pt, -1);
737  }
738 
739  // Detect single patch situation (to avoid allocation)
740  label patch0 = patchIDs[0];
741 
742  for (label i = 1; i < patchIDs.size(); i++)
743  {
744  if (patchIDs[i] != patch0)
745  {
746  patch0 = -1;
747  break;
748  }
749  }
750 
751  if (patch0 >= 0)
752  {
753  // Single patch
754  return pointIndexHit(false, pt, -1);
755  }
756  else
757  {
758  if (surfaceNormals.size() == 1)
759  {
760  // Same normals plane, flat region edge.
761  return pointIndexHit(true, pt, 1);
762  }
763  else
764  {
765  // Detect per normals bin
766  labelList normalToPatch(surfaceNormals.size(), -1);
767  forAll(faceToNormalBin, i)
768  {
769  if (faceToNormalBin[i] != -1)
770  {
771  label& patch = normalToPatch[faceToNormalBin[i]];
772  if (patch == -1)
773  {
774  // First occurence
775  patch = patchIDs[i];
776  }
777  else if (patch == -2)
778  {
779  // Already marked as being on multiple patches
780  }
781  else if (patch != patchIDs[i])
782  {
783  // Mark as being on multiple patches
784  patch = -2;
785  }
786  }
787  }
788 
789  forAll(normalToPatch, normalI)
790  {
791  if (normalToPatch[normalI] == -2)
792  {
793  // Multiple patches on same normals plane, flat region
794  // edge
795  return pointIndexHit(true, pt, 1);
796  }
797  }
798 
799  // All patches on either side of geometric feature anyway
800  return pointIndexHit(true, pt, 0);
801  }
802  }
803 }
804 
805 
806 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
807 (
808  const label iter,
809  const scalar featureCos,
810 
811  const indirectPrimitivePatch& pp,
812  const scalarField& snapDist,
813  const vectorField& nearestDisp,
814  const label pointi,
815 
816  const List<List<point>>& pointFaceSurfNormals,
817  const List<List<point>>& pointFaceDisp,
818  const List<List<point>>& pointFaceCentres,
819  const labelListList& pointFacePatchID,
820 
821  DynamicList<point>& surfacePoints,
822  DynamicList<vector>& surfaceNormals,
823  labelList& faceToNormalBin,
824 
825  vector& patchAttraction,
826  pointConstraint& patchConstraint
827 ) const
828 {
829  patchAttraction = Zero;
830  patchConstraint = pointConstraint();
831 
832  const List<point>& pfSurfNormals = pointFaceSurfNormals[pointi];
833  const List<point>& pfDisp = pointFaceDisp[pointi];
834  const List<point>& pfCentres = pointFaceCentres[pointi];
835 
836  // Bin according to surface normal
837  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
838 
839  //- Bins of differing normals:
840  // - one normal : flat(tish) surface
841  // - two normals : geometric feature edge
842  // - three normals: geometric feature point
843  // - four normals : too complex a feature
844  surfacePoints.clear();
845  surfaceNormals.clear();
846 
847  //- From face to above normals bin
848  faceToNormalBin.setSize(pfDisp.size());
849  faceToNormalBin = -1;
850 
851  forAll(pfSurfNormals, i)
852  {
853  const point& fc = pfCentres[i];
854  const vector& fSNormal = pfSurfNormals[i];
855  const vector& fDisp = pfDisp[i];
856 
857  // What to do with very far attraction? For now just ignore the face
858  if (magSqr(fDisp) < sqr(snapDist[pointi]) && mag(fSNormal) > VSMALL)
859  {
860  const point pt = fc + fDisp;
861 
862  // Do we already have surface normal?
863  faceToNormalBin[i] = findNormal
864  (
865  featureCos,
866  fSNormal,
867  surfaceNormals
868  );
869 
870  if (faceToNormalBin[i] != -1)
871  {
872  // Same normal
873  }
874  else
875  {
876  // Now check if the planes go through the same edge or point
877 
878  if (surfacePoints.size() <= 1)
879  {
880  surfacePoints.append(pt);
881  faceToNormalBin[i] = surfaceNormals.size();
882  surfaceNormals.append(fSNormal);
883  }
884  else if (surfacePoints.size() == 2)
885  {
886  plane pl0(surfacePoints[0], surfaceNormals[0]);
887  plane pl1(surfacePoints[1], surfaceNormals[1]);
888  plane::ray r(pl0.planeIntersect(pl1));
889  vector featureNormal = r.dir() / mag(r.dir());
890 
891  if (mag(fSNormal&featureNormal) >= 0.001)
892  {
893  // Definitely makes a feature point
894  surfacePoints.append(pt);
895  faceToNormalBin[i] = surfaceNormals.size();
896  surfaceNormals.append(fSNormal);
897  }
898  }
899  else if (surfacePoints.size() == 3)
900  {
901  // Have already feature point. See if this new plane is
902  // the same point or not.
903  plane pl0(surfacePoints[0], surfaceNormals[0]);
904  plane pl1(surfacePoints[1], surfaceNormals[1]);
905  plane pl2(surfacePoints[2], surfaceNormals[2]);
906  point p012(pl0.planePlaneIntersect(pl1, pl2));
907 
908  plane::ray r(pl0.planeIntersect(pl1));
909  vector featureNormal = r.dir() / mag(r.dir());
910  if (mag(fSNormal&featureNormal) >= 0.001)
911  {
912  plane pl3(pt, fSNormal);
913  point p013(pl0.planePlaneIntersect(pl1, pl3));
914 
915  if (mag(p012-p013) > snapDist[pointi])
916  {
917  // Different feature point
918  surfacePoints.append(pt);
919  faceToNormalBin[i] = surfaceNormals.size();
920  surfaceNormals.append(fSNormal);
921  }
922  }
923  }
924  }
925  }
926  }
927 
928 
929  const point& pt = pp.localPoints()[pointi];
930 
931  // Check the number of directions
932  if (surfaceNormals.size() == 1)
933  {
934  // Normal distance to plane
935  vector d =
936  ((surfacePoints[0]-pt) & surfaceNormals[0])
937  *surfaceNormals[0];
938 
939  // Trim to snap distance
940  if (magSqr(d) > sqr(snapDist[pointi]))
941  {
942  d *= Foam::sqrt(sqr(snapDist[pointi])/magSqr(d));
943  }
944 
945  patchAttraction = d;
946 
947  // Store constraints
948  patchConstraint.applyConstraint(surfaceNormals[0]);
949  }
950  else if (surfaceNormals.size() == 2)
951  {
952  plane pl0(surfacePoints[0], surfaceNormals[0]);
953  plane pl1(surfacePoints[1], surfaceNormals[1]);
954  plane::ray r(pl0.planeIntersect(pl1));
955  vector n = r.dir() / mag(r.dir());
956 
957  // Get nearest point on infinite ray
958  vector d = r.refPoint()-pt;
959  d -= (d&n)*n;
960 
961  // Trim to snap distance
962  if (magSqr(d) > sqr(snapDist[pointi]))
963  {
964  d *= Foam::sqrt(sqr(snapDist[pointi])/magSqr(d));
965  }
966 
967  patchAttraction = d;
968 
969  // Store constraints
970  patchConstraint.applyConstraint(surfaceNormals[0]);
971  patchConstraint.applyConstraint(surfaceNormals[1]);
972  }
973  else if (surfaceNormals.size() == 3)
974  {
975  // Calculate point from the faces.
976  plane pl0(surfacePoints[0], surfaceNormals[0]);
977  plane pl1(surfacePoints[1], surfaceNormals[1]);
978  plane pl2(surfacePoints[2], surfaceNormals[2]);
979  point cornerPt(pl0.planePlaneIntersect(pl1, pl2));
980  vector d = cornerPt - pt;
981 
982  // Trim to snap distance
983  if (magSqr(d) > sqr(snapDist[pointi]))
984  {
985  d *= Foam::sqrt(sqr(snapDist[pointi])/magSqr(d));
986  }
987 
988  patchAttraction = d;
989 
990  // Store constraints
991  patchConstraint.applyConstraint(surfaceNormals[0]);
992  patchConstraint.applyConstraint(surfaceNormals[1]);
993  patchConstraint.applyConstraint(surfaceNormals[2]);
994  }
995 }
996 
997 
998 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
999 (
1000  const label iter,
1001  const bool avoidSnapProblems,
1002  const scalar featureCos,
1003 
1004  const indirectPrimitivePatch& pp,
1005  const scalarField& snapDist,
1006  const vectorField& nearestDisp,
1007 
1008  const List<List<point>>& pointFaceSurfNormals,
1009  const List<List<point>>& pointFaceDisp,
1010  const List<List<point>>& pointFaceCentres,
1011  const labelListList& pointFacePatchID,
1012 
1013  vectorField& patchAttraction,
1014  List<pointConstraint>& patchConstraints
1015 ) const
1016 {
1017  autoPtr<OBJstream> feStr;
1018  autoPtr<OBJstream> fpStr;
1019  if (debug&meshRefinement::ATTRACTION)
1020  {
1021  feStr.reset
1022  (
1023  new OBJstream
1024  (
1025  meshRefiner_.mesh().time().path()
1026  / "implicitFeatureEdge_" + name(iter) + ".obj"
1027  )
1028  );
1029  Info<< "Dumping implicit feature-edge direction to "
1030  << feStr().name() << endl;
1031 
1032  fpStr.reset
1033  (
1034  new OBJstream
1035  (
1036  meshRefiner_.mesh().time().path()
1037  / "implicitFeaturePoint_" + name(iter) + ".obj"
1038  )
1039  );
1040  Info<< "Dumping implicit feature-point direction to "
1041  << fpStr().name() << endl;
1042  }
1043 
1044 
1045  DynamicList<point> surfacePoints(4);
1046  DynamicList<vector> surfaceNormals(4);
1047  labelList faceToNormalBin;
1048 
1049  forAll(pp.localPoints(), pointi)
1050  {
1051  vector attraction = Zero;
1052  pointConstraint constraint;
1053 
1054  featureAttractionUsingReconstruction
1055  (
1056  iter,
1057  featureCos,
1058 
1059  pp,
1060  snapDist,
1061  nearestDisp,
1062 
1063  pointi,
1064 
1065  pointFaceSurfNormals,
1066  pointFaceDisp,
1067  pointFaceCentres,
1068  pointFacePatchID,
1069 
1070  surfacePoints,
1071  surfaceNormals,
1072  faceToNormalBin,
1073 
1074  attraction,
1075  constraint
1076  );
1077 
1078  if
1079  (
1080  (constraint.first() > patchConstraints[pointi].first())
1081  || (
1082  (constraint.first() == patchConstraints[pointi].first())
1083  && (magSqr(attraction) < magSqr(patchAttraction[pointi]))
1084  )
1085  )
1086  {
1087  patchAttraction[pointi] = attraction;
1088  patchConstraints[pointi] = constraint;
1089 
1090  const point& pt = pp.localPoints()[pointi];
1091 
1092  if (patchConstraints[pointi].first() == 2 && feStr.valid())
1093  {
1094  feStr().write(linePointRef(pt, pt+patchAttraction[pointi]));
1095  }
1096  else if (patchConstraints[pointi].first() == 3 && fpStr.valid())
1097  {
1098  fpStr().write(linePointRef(pt, pt+patchAttraction[pointi]));
1099  }
1100  }
1101  }
1102 }
1103 
1104 
1105 void Foam::snappySnapDriver::stringFeatureEdges
1106 (
1107  const label iter,
1108  const scalar featureCos,
1109 
1110  const indirectPrimitivePatch& pp,
1111  const scalarField& snapDist,
1112 
1113  const vectorField& rawPatchAttraction,
1114  const List<pointConstraint>& rawPatchConstraints,
1115 
1116  vectorField& patchAttraction,
1117  List<pointConstraint>& patchConstraints
1118 ) const
1119 {
1120  // Snap edges to feature edges
1121  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1122  // Walk existing edges and snap remaining ones (that are marked as
1123  // feature edges in rawPatchConstraints)
1124 
1125  // What this does is fill in any faces where not all points
1126  // on the face are being attracted:
1127  /*
1128  +
1129  / \
1130  / \
1131  ---+ +---
1132  \ /
1133  \ /
1134  +
1135  */
1136  // so the top and bottom will never get attracted since the nearest
1137  // back from the feature edge will always be one of the left or right
1138  // points since the face is diamond like. So here we walk the feature edges
1139  // and add any non-attracted points.
1140 
1141 
1142  while (true)
1143  {
1144  label nChanged = 0;
1145 
1146  const labelListList& pointEdges = pp.pointEdges();
1147  forAll(pointEdges, pointi)
1148  {
1149  if (patchConstraints[pointi].first() == 2)
1150  {
1151  const point& pt = pp.localPoints()[pointi];
1152  const labelList& pEdges = pointEdges[pointi];
1153  const vector& featVec = patchConstraints[pointi].second();
1154 
1155  // Detect whether there are edges in both directions.
1156  // (direction along the feature edge that is)
1157  bool hasPos = false;
1158  bool hasNeg = false;
1159 
1160  forAll(pEdges, pEdgeI)
1161  {
1162  const edge& e = pp.edges()[pEdges[pEdgeI]];
1163  label nbrPointi = e.otherVertex(pointi);
1164 
1165  if (patchConstraints[nbrPointi].first() > 1)
1166  {
1167  const point& nbrPt = pp.localPoints()[nbrPointi];
1168  const point featPt =
1169  nbrPt + patchAttraction[nbrPointi];
1170  const scalar cosAngle = (featVec & (featPt-pt));
1171 
1172  if (cosAngle > 0)
1173  {
1174  hasPos = true;
1175  }
1176  else
1177  {
1178  hasNeg = true;
1179  }
1180  }
1181  }
1182 
1183  if (!hasPos || !hasNeg)
1184  {
1185  //Pout<< "**Detected feature string end at "
1186  // << pp.localPoints()[pointi] << endl;
1187 
1188  // No string. Assign best choice on either side
1189  label bestPosPointi = -1;
1190  scalar minPosDistSqr = GREAT;
1191  label bestNegPointi = -1;
1192  scalar minNegDistSqr = GREAT;
1193 
1194  forAll(pEdges, pEdgeI)
1195  {
1196  const edge& e = pp.edges()[pEdges[pEdgeI]];
1197  label nbrPointi = e.otherVertex(pointi);
1198 
1199  if
1200  (
1201  patchConstraints[nbrPointi].first() <= 1
1202  && rawPatchConstraints[nbrPointi].first() > 1
1203  )
1204  {
1205  const vector& nbrFeatVec =
1206  rawPatchConstraints[pointi].second();
1207 
1208  if (mag(featVec&nbrFeatVec) > featureCos)
1209  {
1210  // nbrPointi attracted to sameish feature
1211  // Note: also check on position.
1212 
1213  scalar d2 = magSqr
1214  (
1215  rawPatchAttraction[nbrPointi]
1216  );
1217 
1218  const point featPt =
1219  pp.localPoints()[nbrPointi]
1220  + rawPatchAttraction[nbrPointi];
1221  const scalar cosAngle =
1222  (featVec & (featPt-pt));
1223 
1224  if (cosAngle > 0)
1225  {
1226  if (!hasPos && d2 < minPosDistSqr)
1227  {
1228  minPosDistSqr = d2;
1229  bestPosPointi = nbrPointi;
1230  }
1231  }
1232  else
1233  {
1234  if (!hasNeg && d2 < minNegDistSqr)
1235  {
1236  minNegDistSqr = d2;
1237  bestNegPointi = nbrPointi;
1238  }
1239  }
1240  }
1241  }
1242  }
1243 
1244  if (bestPosPointi != -1)
1245  {
1246  // Use reconstructed-feature attraction. Use only
1247  // part of it since not sure...
1248  //const point& bestPt =
1249  // pp.localPoints()[bestPosPointi];
1250  //Pout<< "**Overriding point " << bestPt
1251  // << " on reconstructed feature edge at "
1252  // << rawPatchAttraction[bestPosPointi]+bestPt
1253  // << " to attracted-to-feature-edge." << endl;
1254  patchAttraction[bestPosPointi] =
1255  0.5*rawPatchAttraction[bestPosPointi];
1256  patchConstraints[bestPosPointi] =
1257  rawPatchConstraints[bestPosPointi];
1258 
1259  nChanged++;
1260  }
1261  if (bestNegPointi != -1)
1262  {
1263  // Use reconstructed-feature attraction. Use only
1264  // part of it since not sure...
1265  //const point& bestPt =
1266  // pp.localPoints()[bestNegPointi];
1267  //Pout<< "**Overriding point " << bestPt
1268  // << " on reconstructed feature edge at "
1269  // << rawPatchAttraction[bestNegPointi]+bestPt
1270  // << " to attracted-to-feature-edge." << endl;
1271  patchAttraction[bestNegPointi] =
1272  0.5*rawPatchAttraction[bestNegPointi];
1273  patchConstraints[bestNegPointi] =
1274  rawPatchConstraints[bestNegPointi];
1275 
1276  nChanged++;
1277  }
1278  }
1279  }
1280  }
1281 
1282 
1283  reduce(nChanged, sumOp<label>());
1284  Info<< "Stringing feature edges : changed " << nChanged << " points"
1285  << endl;
1286  if (nChanged == 0)
1287  {
1288  break;
1289  }
1290  }
1291 }
1292 
1293 
1294 void Foam::snappySnapDriver::releasePointsNextToMultiPatch
1295 (
1296  const label iter,
1297  const scalar featureCos,
1298 
1299  const indirectPrimitivePatch& pp,
1300  const scalarField& snapDist,
1301 
1302  const List<List<point>>& pointFaceCentres,
1303  const labelListList& pointFacePatchID,
1304 
1305  const vectorField& rawPatchAttraction,
1306  const List<pointConstraint>& rawPatchConstraints,
1307 
1308  vectorField& patchAttraction,
1309  List<pointConstraint>& patchConstraints
1310 ) const
1311 {
1312  autoPtr<OBJstream> multiPatchStr;
1313  if (debug&meshRefinement::ATTRACTION)
1314  {
1315  multiPatchStr.reset
1316  (
1317  new OBJstream
1318  (
1319  meshRefiner_.mesh().time().path()
1320  / "multiPatch_" + name(iter) + ".obj"
1321  )
1322  );
1323  Info<< "Dumping removed constraints due to same-face"
1324  << " multi-patch points to "
1325  << multiPatchStr().name() << endl;
1326  }
1327 
1328 
1329  // 1. Mark points on multiple patches
1330  PackedBoolList isMultiPatchPoint(pp.size());
1331 
1332  forAll(pointFacePatchID, pointi)
1333  {
1334  pointIndexHit multiPatchPt = findMultiPatchPoint
1335  (
1336  pp.localPoints()[pointi],
1337  pointFacePatchID[pointi],
1338  pointFaceCentres[pointi]
1339  );
1340  isMultiPatchPoint[pointi] = multiPatchPt.hit();
1341  }
1342 
1343  // 2. Make sure multi-patch points are also attracted
1344  forAll(isMultiPatchPoint, pointi)
1345  {
1346  if (isMultiPatchPoint[pointi])
1347  {
1348  if
1349  (
1350  patchConstraints[pointi].first() <= 1
1351  && rawPatchConstraints[pointi].first() > 1
1352  )
1353  {
1354  patchAttraction[pointi] = rawPatchAttraction[pointi];
1355  patchConstraints[pointi] = rawPatchConstraints[pointi];
1356 
1357  //if (multiPatchStr.valid())
1358  //{
1359  // Pout<< "Adding constraint on multiPatchPoint:"
1360  // << pp.localPoints()[pointi]
1361  // << " constraint:" << patchConstraints[pointi]
1362  // << " attraction:" << patchAttraction[pointi]
1363  // << endl;
1364  //}
1365  }
1366  }
1367  }
1368 
1369  // Up to here it is all parallel ok.
1370 
1371 
1372  // 3. Knock out any attraction on faces with multi-patch points
1373  label nChanged = 0;
1374  forAll(pp.localFaces(), facei)
1375  {
1376  const face& f = pp.localFaces()[facei];
1377 
1378  label nMultiPatchPoints = 0;
1379  forAll(f, fp)
1380  {
1381  label pointi = f[fp];
1382  if
1383  (
1384  isMultiPatchPoint[pointi]
1385  && patchConstraints[pointi].first() > 1
1386  )
1387  {
1388  nMultiPatchPoints++;
1389  }
1390  }
1391 
1392  if (nMultiPatchPoints > 0)
1393  {
1394  forAll(f, fp)
1395  {
1396  label pointi = f[fp];
1397  if
1398  (
1399  !isMultiPatchPoint[pointi]
1400  && patchConstraints[pointi].first() > 1
1401  )
1402  {
1403  //Pout<< "Knocking out constraint"
1404  // << " on non-multiPatchPoint:"
1405  // << pp.localPoints()[pointi] << endl;
1406  patchAttraction[pointi] = Zero;
1407  patchConstraints[pointi] = pointConstraint();
1408  nChanged++;
1409 
1410  if (multiPatchStr.valid())
1411  {
1412  multiPatchStr().write(pp.localPoints()[pointi]);
1413  }
1414  }
1415  }
1416  }
1417  }
1418 
1419  reduce(nChanged, sumOp<label>());
1420  Info<< "Removing constraints near multi-patch points : changed "
1421  << nChanged << " points" << endl;
1422 }
1423 
1424 
1425 Foam::labelPair Foam::snappySnapDriver::findDiagonalAttraction
1426 (
1427  const indirectPrimitivePatch& pp,
1428  const vectorField& patchAttraction,
1429  const List<pointConstraint>& patchConstraints,
1430  const label facei
1431 ) const
1432 {
1433  const face& f = pp.localFaces()[facei];
1434 
1435  // For now just detect any attraction. Improve this to look at
1436  // actual attraction position and orientation
1437 
1438  labelPair attractIndices(-1, -1);
1439 
1440  if (f.size() >= 4)
1441  {
1442  forAll(f, fp)
1443  {
1444  label pointi = f[fp];
1445  if (patchConstraints[pointi].first() >= 2)
1446  {
1447  // Attract to feature edge or point
1448  if (attractIndices[0] == -1)
1449  {
1450  // First attraction. Store
1451  attractIndices[0] = fp;
1452  }
1453  else if (attractIndices[1] == -1)
1454  {
1455  // Second attraction. Check if not consecutive to first
1456  // attraction
1457  label fp0 = attractIndices[0];
1458  if (f.fcIndex(fp0) == fp || f.fcIndex(fp) == fp0)
1459  {
1460  // Consecutive. Skip.
1461  attractIndices = labelPair(-1, -1);
1462  break;
1463  }
1464  else
1465  {
1466  attractIndices[1] = fp;
1467  }
1468  }
1469  else
1470  {
1471  // More than two attractions. Skip.
1472  attractIndices = labelPair(-1, -1);
1473  break;
1474  }
1475  }
1476  }
1477 
1478 
1479  if (attractIndices[1] == -1)
1480  {
1481  // Found only one attraction. Skip.
1482  attractIndices = labelPair(-1, -1);
1483  }
1484  }
1485  return attractIndices;
1486 }
1487 
1488 
1489 void Foam::snappySnapDriver::avoidDiagonalAttraction
1490 (
1491  const label iter,
1492  const scalar featureCos,
1493 
1494  const indirectPrimitivePatch& pp,
1495 
1496  vectorField& patchAttraction,
1497  List<pointConstraint>& patchConstraints
1498 ) const
1499 {
1500  forAll(pp.localFaces(), facei)
1501  {
1502  const face& f = pp.localFaces()[facei];
1503 
1504  labelPair diag = findDiagonalAttraction
1505  (
1506  pp,
1507  patchAttraction,
1508  patchConstraints,
1509  facei
1510  );
1511 
1512  if (diag[0] != -1 && diag[1] != -1)
1513  {
1514  // Found two diagonal points that being attracted.
1515  // For now just attract my one to the average of those.
1516  const label i0 = f[diag[0]];
1517  const point pt0 =
1518  pp.localPoints()[i0]+patchAttraction[i0];
1519  const label i1 = f[diag[1]];
1520  const point pt1 =
1521  pp.localPoints()[i1]+patchAttraction[i1];
1522  const point mid = 0.5*(pt0+pt1);
1523 
1524  const scalar cosAngle = mag
1525  (
1526  patchConstraints[i0].second()
1527  & patchConstraints[i1].second()
1528  );
1529 
1530  //Pout<< "Found diagonal attraction at indices:"
1531  // << diag[0]
1532  // << " and " << diag[1]
1533  // << " with cosAngle:" << cosAngle
1534  // << " mid:" << mid << endl;
1535 
1536  if (cosAngle > featureCos)
1537  {
1538  // Add the nearest of the other two points as
1539  // attractor
1540  label minFp = -1;
1541  scalar minDistSqr = GREAT;
1542  forAll(f, fp)
1543  {
1544  label pointi = f[fp];
1545  if (patchConstraints[pointi].first() <= 1)
1546  {
1547  const point& pt = pp.localPoints()[pointi];
1548  scalar distSqr = magSqr(mid-pt);
1549  if (distSqr < minDistSqr)
1550  {
1551  distSqr = minDistSqr;
1552  minFp = fp;
1553  }
1554  }
1555  }
1556  if (minFp != -1)
1557  {
1558  label minPointi = f[minFp];
1559  patchAttraction[minPointi] =
1560  mid-pp.localPoints()[minPointi];
1561  patchConstraints[minPointi] =
1562  patchConstraints[f[diag[0]]];
1563  }
1564  }
1565  else
1566  {
1567  //Pout<< "Diagonal attractors at" << nl
1568  // << " pt0:" << pt0
1569  // << " constraint:"
1570  // << patchConstraints[i0].second() << nl
1571  // << " pt1:" << pt1
1572  // << " constraint:"
1573  // << patchConstraints[i1].second() << nl
1574  // << " make too large an angle:"
1575  // << mag
1576  // (
1577  // patchConstraints[i0].second()
1578  // & patchConstraints[i1].second()
1579  // )
1580  // << endl;
1581  }
1582  }
1583  }
1584 }
1585 
1586 
1588 Foam::snappySnapDriver::findNearFeatureEdge
1589 (
1590  const bool isRegionEdge,
1591 
1592  const indirectPrimitivePatch& pp,
1593  const scalarField& snapDist,
1594  const label pointi,
1595  const point& estimatedPt,
1596 
1597  List<List<DynamicList<point>>>& edgeAttractors,
1598  List<List<DynamicList<pointConstraint>>>& edgeConstraints,
1599  vectorField& patchAttraction,
1600  List<pointConstraint>& patchConstraints
1601 ) const
1602 {
1603  const refinementFeatures& features = meshRefiner_.features();
1604 
1605  labelList nearEdgeFeat;
1606  List<pointIndexHit> nearEdgeInfo;
1607  vectorField nearNormal;
1608 
1609  if (isRegionEdge)
1610  {
1611  features.findNearestRegionEdge
1612  (
1613  pointField(1, estimatedPt),
1614  scalarField(1, sqr(snapDist[pointi])),
1615  nearEdgeFeat,
1616  nearEdgeInfo,
1617  nearNormal
1618  );
1619  }
1620  else
1621  {
1622  features.findNearestEdge
1623  (
1624  pointField(1, estimatedPt),
1625  scalarField(1, sqr(snapDist[pointi])),
1626  nearEdgeFeat,
1627  nearEdgeInfo,
1628  nearNormal
1629  );
1630  }
1631 
1632  const pointIndexHit& nearInfo = nearEdgeInfo[0];
1633  label featI = nearEdgeFeat[0];
1634 
1635  if (nearInfo.hit())
1636  {
1637  // So we have a point on the feature edge. Use this
1638  // instead of our estimate from planes.
1639  edgeAttractors[featI][nearInfo.index()].append
1640  (
1641  nearInfo.hitPoint()
1642  );
1643  pointConstraint c(Tuple2<label, vector>(2, nearNormal[0]));
1644  edgeConstraints[featI][nearInfo.index()].append(c);
1645 
1646  // Store for later use
1647  patchAttraction[pointi] =
1648  nearInfo.hitPoint()-pp.localPoints()[pointi];
1649  patchConstraints[pointi] = c;
1650  }
1651  return Tuple2<label, pointIndexHit>(featI, nearInfo);
1652 }
1653 
1654 
1656 Foam::snappySnapDriver::findNearFeaturePoint
1657 (
1658  const bool isRegionPoint,
1659 
1660  const indirectPrimitivePatch& pp,
1661  const scalarField& snapDist,
1662  const label pointi,
1663  const point& estimatedPt,
1664 
1665  // Feature-point to pp point
1666  List<labelList>& pointAttractor,
1668  // Feature-edge to pp point
1669  List<List<DynamicList<point>>>& edgeAttractors,
1670  List<List<DynamicList<pointConstraint>>>& edgeConstraints,
1671  // pp point to nearest feature
1672  vectorField& patchAttraction,
1673  List<pointConstraint>& patchConstraints
1674 ) const
1675 {
1676  const refinementFeatures& features = meshRefiner_.features();
1677 
1678  labelList nearFeat;
1680  features.findNearestPoint
1681  (
1682  pointField(1, estimatedPt),
1683  scalarField(1, sqr(snapDist[pointi])),
1684  nearFeat,
1685  nearInfo
1686  );
1687 
1688  label featI = nearFeat[0];
1689 
1690  if (featI != -1)
1691  {
1692  const point& pt = pp.localPoints()[pointi];
1693 
1694  label featPointi = nearInfo[0].index();
1695  const point& featPt = nearInfo[0].hitPoint();
1696  scalar distSqr = magSqr(featPt-pt);
1697 
1698  // Check if already attracted
1699  label oldPointi = pointAttractor[featI][featPointi];
1700 
1701  if (oldPointi != -1)
1702  {
1703  // Check distance
1704  if (distSqr >= magSqr(featPt-pp.localPoints()[oldPointi]))
1705  {
1706  // oldPointi nearest. Keep.
1707  featI = -1;
1708  featPointi = -1;
1709  }
1710  else
1711  {
1712  // Current pointi nearer.
1713  pointAttractor[featI][featPointi] = pointi;
1714  pointConstraints[featI][featPointi].first() = 3;
1715  pointConstraints[featI][featPointi].second() = Zero;
1716 
1717  // Store for later use
1718  patchAttraction[pointi] = featPt-pt;
1719  patchConstraints[pointi] =
1720  pointConstraints[featI][featPointi];
1721 
1722  // Reset oldPointi to nearest on feature edge
1723  patchAttraction[oldPointi] = Zero;
1724  patchConstraints[oldPointi] = pointConstraint();
1725 
1726  findNearFeatureEdge
1727  (
1728  isRegionPoint, // search region edges only
1729 
1730  pp,
1731  snapDist,
1732  oldPointi,
1733  pp.localPoints()[oldPointi],
1734 
1735  edgeAttractors,
1736  edgeConstraints,
1737  patchAttraction,
1738  patchConstraints
1739  );
1740  }
1741  }
1742  else
1743  {
1744  // Current pointi nearer.
1745  pointAttractor[featI][featPointi] = pointi;
1746  pointConstraints[featI][featPointi].first() = 3;
1747  pointConstraints[featI][featPointi].second() = Zero;
1748 
1749  // Store for later use
1750  patchAttraction[pointi] = featPt-pt;
1751  patchConstraints[pointi] = pointConstraints[featI][featPointi];
1752  }
1753  }
1754 
1755  return Tuple2<label, pointIndexHit>(featI, nearInfo[0]);
1756 }
1757 
1758 
1759 void Foam::snappySnapDriver::determineFeatures
1760 (
1761  const label iter,
1762  const scalar featureCos,
1763  const bool multiRegionFeatureSnap,
1764 
1765  const indirectPrimitivePatch& pp,
1766  const scalarField& snapDist,
1767  const vectorField& nearestDisp,
1768 
1769  const List<List<point>>& pointFaceSurfNormals,
1770  const List<List<point>>& pointFaceDisp,
1771  const List<List<point>>& pointFaceCentres,
1772  const labelListList& pointFacePatchID,
1773 
1774  // Feature-point to pp point
1775  List<labelList>& pointAttractor,
1777  // Feature-edge to pp point
1778  List<List<DynamicList<point>>>& edgeAttractors,
1779  List<List<DynamicList<pointConstraint>>>& edgeConstraints,
1780  // pp point to nearest feature
1781  vectorField& patchAttraction,
1782  List<pointConstraint>& patchConstraints
1783 ) const
1784 {
1785  autoPtr<OBJstream> featureEdgeStr;
1786  autoPtr<OBJstream> missedEdgeStr;
1787  autoPtr<OBJstream> featurePointStr;
1788 
1789  if (debug&meshRefinement::ATTRACTION)
1790  {
1791  featureEdgeStr.reset
1792  (
1793  new OBJstream
1794  (
1795  meshRefiner_.mesh().time().path()
1796  / "featureEdge_" + name(iter) + ".obj"
1797  )
1798  );
1799  Info<< "Dumping feature-edge sampling to "
1800  << featureEdgeStr().name() << endl;
1801 
1802  missedEdgeStr.reset
1803  (
1804  new OBJstream
1805  (
1806  meshRefiner_.mesh().time().path()
1807  / "missedFeatureEdge_" + name(iter) + ".obj"
1808  )
1809  );
1810  Info<< "Dumping feature-edges that are too far away to "
1811  << missedEdgeStr().name() << endl;
1812 
1813  featurePointStr.reset
1814  (
1815  new OBJstream
1816  (
1817  meshRefiner_.mesh().time().path()
1818  / "featurePoint_" + name(iter) + ".obj"
1819  )
1820  );
1821  Info<< "Dumping feature-point sampling to "
1822  << featurePointStr().name() << endl;
1823  }
1824 
1825 
1826  DynamicList<point> surfacePoints(4);
1827  DynamicList<vector> surfaceNormals(4);
1828  labelList faceToNormalBin;
1829 
1830  forAll(pp.localPoints(), pointi)
1831  {
1832  const point& pt = pp.localPoints()[pointi];
1833 
1834  vector attraction = Zero;
1835  pointConstraint constraint;
1836 
1837  featureAttractionUsingReconstruction
1838  (
1839  iter,
1840  featureCos,
1841 
1842  pp,
1843  snapDist,
1844  nearestDisp,
1845 
1846  pointi,
1847 
1848  pointFaceSurfNormals,
1849  pointFaceDisp,
1850  pointFaceCentres,
1851  pointFacePatchID,
1852 
1853  surfacePoints,
1854  surfaceNormals,
1855  faceToNormalBin,
1856 
1857  attraction,
1858  constraint
1859  );
1860 
1861 
1862  if
1863  (
1864  (constraint.first() > patchConstraints[pointi].first())
1865  || (
1866  (constraint.first() == patchConstraints[pointi].first())
1867  && (magSqr(attraction) < magSqr(patchAttraction[pointi]))
1868  )
1869  )
1870  {
1871  patchAttraction[pointi] = attraction;
1872  patchConstraints[pointi] = constraint;
1873 
1874  // Check the number of directions
1875  if (patchConstraints[pointi].first() == 1)
1876  {
1877  // Flat surface. Check for different patchIDs
1878  if (multiRegionFeatureSnap)
1879  {
1880  const point estimatedPt(pt + nearestDisp[pointi]);
1881  pointIndexHit multiPatchPt
1882  (
1883  findMultiPatchPoint
1884  (
1885  estimatedPt,
1886  pointFacePatchID[pointi],
1887  surfaceNormals,
1888  faceToNormalBin
1889  )
1890  );
1891 
1892  if (multiPatchPt.hit())
1893  {
1894  // Behave like when having two surface normals so
1895  // attract to nearest feature edge (with a guess for
1896  // the multipatch point as starting point)
1898  findNearFeatureEdge
1899  (
1900  true, // isRegionPoint
1901  pp,
1902  snapDist,
1903  pointi,
1904  multiPatchPt.hitPoint(), // estimatedPt
1905 
1906  edgeAttractors,
1907  edgeConstraints,
1908 
1909  patchAttraction,
1910  patchConstraints
1911  );
1912 
1913  const pointIndexHit& info = nearInfo.second();
1914  if (info.hit())
1915  {
1916  // Dump
1917  if (featureEdgeStr.valid())
1918  {
1919  featureEdgeStr().write
1920  (
1921  linePointRef(pt, info.hitPoint())
1922  );
1923  }
1924  }
1925  else
1926  {
1927  if (missedEdgeStr.valid())
1928  {
1929  missedEdgeStr().write
1930  (
1931  linePointRef(pt, info.missPoint())
1932  );
1933  }
1934  }
1935  }
1936  }
1937  }
1938  else if (patchConstraints[pointi].first() == 2)
1939  {
1940  // Mark point on the nearest feature edge. Note that we
1941  // only search within the surrounding since the plane
1942  // reconstruction might find a feature where there isn't one.
1943  const point estimatedPt(pt + patchAttraction[pointi]);
1944 
1946 
1947  // Geometric feature edge. Check for different patchIDs
1948  bool hasSnapped = false;
1949  if (multiRegionFeatureSnap)
1950  {
1951  pointIndexHit multiPatchPt
1952  (
1953  findMultiPatchPoint
1954  (
1955  estimatedPt,
1956  pointFacePatchID[pointi],
1957  surfaceNormals,
1958  faceToNormalBin
1959  )
1960  );
1961  if (multiPatchPt.hit())
1962  {
1963  if (multiPatchPt.index() == 0)
1964  {
1965  // Geometric feature edge is also region edge
1966  nearInfo = findNearFeatureEdge
1967  (
1968  true, // isRegionPoint
1969  pp,
1970  snapDist,
1971  pointi,
1972  estimatedPt,
1973 
1974  edgeAttractors,
1975  edgeConstraints,
1976 
1977  patchAttraction,
1978  patchConstraints
1979  );
1980  hasSnapped = true;
1981  }
1982  else
1983  {
1984  // One of planes of feature contains multiple
1985  // regions
1986  nearInfo = findNearFeaturePoint
1987  (
1988  true, // isRegionPoint
1989  pp,
1990  snapDist,
1991  pointi,
1992  estimatedPt,
1993 
1994  // Feature-point to pp point
1995  pointAttractor,
1997  // Feature-edge to pp point
1998  edgeAttractors,
1999  edgeConstraints,
2000  // pp point to nearest feature
2001  patchAttraction,
2002  patchConstraints
2003  );
2004  hasSnapped = true;
2005  }
2006  }
2007  }
2008 
2009  if (!hasSnapped)
2010  {
2011  // Determine nearest point on feature edge. Store
2012  // constraint
2013  // (calculated from feature edge, alternative would be to
2014  // use constraint calculated from both surfaceNormals)
2015  nearInfo = findNearFeatureEdge
2016  (
2017  false, // isRegionPoint
2018  pp,
2019  snapDist,
2020  pointi,
2021  estimatedPt,
2022 
2023  edgeAttractors,
2024  edgeConstraints,
2025 
2026  patchAttraction,
2027  patchConstraints
2028  );
2029  hasSnapped = true;
2030  }
2031 
2032  // Dump to obj
2033  const pointIndexHit& info = nearInfo.second();
2034  if (info.hit())
2035  {
2036  if
2037  (
2038  patchConstraints[pointi].first() == 3
2039  && featurePointStr.valid()
2040  )
2041  {
2042  featurePointStr().write
2043  (
2044  linePointRef(pt, info.hitPoint())
2045  );
2046  }
2047  else if
2048  (
2049  patchConstraints[pointi].first() == 2
2050  && featureEdgeStr.valid()
2051  )
2052  {
2053  featureEdgeStr().write
2054  (
2055  linePointRef(pt, info.hitPoint())
2056  );
2057  }
2058  }
2059  else
2060  {
2061  if (missedEdgeStr.valid())
2062  {
2063  missedEdgeStr().write
2064  (
2065  linePointRef(pt, info.missPoint())
2066  );
2067  }
2068  }
2069  }
2070  else if (patchConstraints[pointi].first() == 3)
2071  {
2072  // Mark point on the nearest feature point.
2073  const point estimatedPt(pt + patchAttraction[pointi]);
2074 
2076 
2077  if (multiRegionFeatureSnap)
2078  {
2079  pointIndexHit multiPatchPt
2080  (
2081  findMultiPatchPoint
2082  (
2083  estimatedPt,
2084  pointFacePatchID[pointi],
2085  surfaceNormals,
2086  faceToNormalBin
2087  )
2088  );
2089  if (multiPatchPt.hit())
2090  {
2091  // Multiple regions
2092  nearInfo = findNearFeaturePoint
2093  (
2094  true, // isRegionPoint
2095  pp,
2096  snapDist,
2097  pointi,
2098  estimatedPt,
2099 
2100  // Feature-point to pp point
2101  pointAttractor,
2103  // Feature-edge to pp point
2104  edgeAttractors,
2105  edgeConstraints,
2106  // pp point to nearest feature
2107  patchAttraction,
2108  patchConstraints
2109  );
2110  }
2111  else
2112  {
2113  nearInfo = findNearFeaturePoint
2114  (
2115  false, // isRegionPoint
2116  pp,
2117  snapDist,
2118  pointi,
2119  estimatedPt,
2120 
2121  // Feature-point to pp point
2122  pointAttractor,
2124  // Feature-edge to pp point
2125  edgeAttractors,
2126  edgeConstraints,
2127  // pp point to nearest feature
2128  patchAttraction,
2129  patchConstraints
2130  );
2131  }
2132  }
2133  else
2134  {
2135  // No multi-patch snapping
2136  nearInfo = findNearFeaturePoint
2137  (
2138  false, // isRegionPoint
2139  pp,
2140  snapDist,
2141  pointi,
2142  estimatedPt,
2143 
2144  // Feature-point to pp point
2145  pointAttractor,
2147  // Feature-edge to pp point
2148  edgeAttractors,
2149  edgeConstraints,
2150  // pp point to nearest feature
2151  patchAttraction,
2152  patchConstraints
2153  );
2154  }
2155 
2156  const pointIndexHit& info = nearInfo.second();
2157  if (info.hit() && featurePointStr.valid())
2158  {
2159  featurePointStr().write
2160  (
2161  linePointRef(pt, info.hitPoint())
2162  );
2163  }
2164  }
2165  }
2166  }
2167 }
2168 
2169 
2170 void Foam::snappySnapDriver::determineBaffleFeatures
2171 (
2172  const label iter,
2173  const scalar featureCos,
2174 
2175  const indirectPrimitivePatch& pp,
2176  const scalarField& snapDist,
2177 
2178  // Feature-point to pp point
2179  List<labelList>& pointAttractor,
2181  // Feature-edge to pp point
2182  List<List<DynamicList<point>>>& edgeAttractors,
2183  List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2184  // pp point to nearest feature
2185  vectorField& patchAttraction,
2186  List<pointConstraint>& patchConstraints
2187 ) const
2188 {
2189  // Override pointAttractor, edgeAttractor, patchAttration etc. to
2190  // implement 'baffle' handling.
2191  // Baffle: the mesh pp point originates from a loose standing
2192  // baffle.
2193  // Sampling the surface with the surrounding face-centres only picks up
2194  // a single triangle normal so above determineFeatures will not have
2195  // detected anything. So explicitly pick up feature edges on the pp
2196  // (after duplicating points & smoothing so will already have been
2197  // expanded) and match these to the features.
2198 
2199  const fvMesh& mesh = meshRefiner_.mesh();
2200  const refinementFeatures& features = meshRefiner_.features();
2201 
2202  // Calculate edge-faces
2203  List<List<point>> edgeFaceNormals(pp.nEdges());
2204 
2205  // Fill local data
2206  forAll(pp.edgeFaces(), edgeI)
2207  {
2208  const labelList& eFaces = pp.edgeFaces()[edgeI];
2209  List<point>& eFc = edgeFaceNormals[edgeI];
2210  eFc.setSize(eFaces.size());
2211  forAll(eFaces, i)
2212  {
2213  label facei = eFaces[i];
2214  eFc[i] = pp.faceNormals()[facei];
2215  }
2216  }
2217 
2218  {
2219  // Precalculate mesh edges for pp.edges.
2220  const labelList meshEdges
2221  (
2222  pp.meshEdges(mesh.edges(), mesh.pointEdges())
2223  );
2225  (
2226  mesh,
2227  meshEdges,
2228  edgeFaceNormals,
2230  List<point>(),
2232  );
2233  }
2234 
2235  // Detect baffle edges. Assume initial mesh will have 0,90 or 180
2236  // (baffle) degree angles so smoothing should make 0,90
2237  // to be less than 90.
2238  const scalar baffleFeatureCos = Foam::cos(degToRad(91));
2239 
2240 
2241  autoPtr<OBJstream> baffleEdgeStr;
2242  if (debug&meshRefinement::ATTRACTION)
2243  {
2244  baffleEdgeStr.reset
2245  (
2246  new OBJstream
2247  (
2248  meshRefiner_.mesh().time().path()
2249  / "baffleEdge_" + name(iter) + ".obj"
2250  )
2251  );
2252  Info<< nl << "Dumping baffle-edges to "
2253  << baffleEdgeStr().name() << endl;
2254  }
2255 
2256 
2257  // Is edge on baffle
2258  PackedBoolList isBaffleEdge(pp.nEdges());
2259  label nBaffleEdges = 0;
2260  // Is point on
2261  // 0 : baffle-edge (0)
2262  // 1 : baffle-feature-point (1)
2263  // -1 : rest
2264  labelList pointStatus(pp.nPoints(), -1);
2265 
2266  forAll(edgeFaceNormals, edgeI)
2267  {
2268  const List<point>& efn = edgeFaceNormals[edgeI];
2269 
2270  if (efn.size() == 2 && (efn[0]&efn[1]) < baffleFeatureCos)
2271  {
2272  isBaffleEdge[edgeI] = true;
2273  nBaffleEdges++;
2274  const edge& e = pp.edges()[edgeI];
2275  pointStatus[e[0]] = 0;
2276  pointStatus[e[1]] = 0;
2277 
2278  if (baffleEdgeStr.valid())
2279  {
2280  const point& p0 = pp.localPoints()[e[0]];
2281  const point& p1 = pp.localPoints()[e[1]];
2282  baffleEdgeStr().write(linePointRef(p0, p1));
2283  }
2284  }
2285  }
2286 
2287  Info<< "Detected "
2288  << returnReduce(nBaffleEdges, sumOp<label>())
2289  << " baffle edges out of "
2290  << returnReduce(pp.nEdges(), sumOp<label>())
2291  << " edges." << endl;
2292 
2293 
2294  forAll(pp.pointEdges(), pointi)
2295  {
2296  if
2297  (
2298  isFeaturePoint
2299  (
2300  featureCos,
2301  pp,
2302  isBaffleEdge,
2303  pointi
2304  )
2305  )
2306  {
2307  //Pout<< "Detected feature point:" << pp.localPoints()[pointi]
2308  // << endl;
2309  //-TEMPORARILY DISABLED:
2310  //pointStatus[pointi] = 1;
2311  }
2312  }
2313 
2314 
2315  forAll(pointStatus, pointi)
2316  {
2317  const point& pt = pp.localPoints()[pointi];
2318 
2319  if (pointStatus[pointi] == 0) // baffle edge
2320  {
2321  Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
2322  (
2323  false, // isRegionPoint?
2324  pp,
2325  snapDist,
2326  pointi,
2327  pt,
2328 
2329  edgeAttractors,
2330  edgeConstraints,
2331  patchAttraction,
2332  patchConstraints
2333  );
2334 
2335  if (!nearInfo.second().hit())
2336  {
2337  //Pout<< "*** Failed to find close edge to point " << pt
2338  // << endl;
2339  }
2340  }
2341  else if (pointStatus[pointi] == 1) // baffle point
2342  {
2343  labelList nearFeat;
2345  features.findNearestPoint
2346  (
2347  pointField(1, pt),
2348  scalarField(1, sqr(snapDist[pointi])),
2349  nearFeat,
2350  nearInfo
2351  );
2352 
2353  label featI = nearFeat[0];
2354 
2355  if (featI != -1)
2356  {
2357  label featPointi = nearInfo[0].index();
2358  const point& featPt = nearInfo[0].hitPoint();
2359  scalar distSqr = magSqr(featPt-pt);
2360 
2361  // Check if already attracted
2362  label oldPointi = pointAttractor[featI][featPointi];
2363 
2364  if
2365  (
2366  oldPointi == -1
2367  || (
2368  distSqr
2369  < magSqr(featPt-pp.localPoints()[oldPointi])
2370  )
2371  )
2372  {
2373  pointAttractor[featI][featPointi] = pointi;
2374  pointConstraints[featI][featPointi].first() = 3;
2375  pointConstraints[featI][featPointi].second() =
2376  Zero;
2377 
2378  // Store for later use
2379  patchAttraction[pointi] = featPt-pt;
2380  patchConstraints[pointi] =
2381  pointConstraints[featI][featPointi];
2382 
2383  if (oldPointi != -1)
2384  {
2385  // The current point is closer so wins. Reset
2386  // the old point to attract to nearest edge
2387  // instead.
2388  findNearFeatureEdge
2389  (
2390  false, // isRegionPoint
2391  pp,
2392  snapDist,
2393  oldPointi,
2394  pp.localPoints()[oldPointi],
2395 
2396  edgeAttractors,
2397  edgeConstraints,
2398  patchAttraction,
2399  patchConstraints
2400  );
2401  }
2402  }
2403  else
2404  {
2405  // Make it fall through to check below
2406  featI = -1;
2407  }
2408  }
2409 
2410  // Not found a feature point or another point is already
2411  // closer to that feature
2412  if (featI == -1)
2413  {
2414  //Pout<< "*** Falling back to finding nearest feature"
2415  // << " edge"
2416  // << " for baffle-feature-point " << pt
2417  // << endl;
2418 
2419  findNearFeatureEdge
2420  (
2421  false, // isRegionPoint
2422  pp,
2423  snapDist,
2424  pointi,
2425  pt, // starting point
2426 
2427  edgeAttractors,
2428  edgeConstraints,
2429  patchAttraction,
2430  patchConstraints
2431  );
2432  }
2433  }
2434  }
2435 }
2436 
2437 
2438 void Foam::snappySnapDriver::reverseAttractMeshPoints
2439 (
2440  const label iter,
2441 
2442  const indirectPrimitivePatch& pp,
2443  const scalarField& snapDist,
2444 
2445  // Feature-point to pp point
2446  const List<labelList>& pointAttractor,
2448  // Feature-edge to pp point
2449  const List<List<DynamicList<point>>>& edgeAttractors,
2450  const List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2451 
2452  const vectorField& rawPatchAttraction,
2453  const List<pointConstraint>& rawPatchConstraints,
2454 
2455  // pp point to nearest feature
2456  vectorField& patchAttraction,
2457  List<pointConstraint>& patchConstraints
2458 ) const
2459 {
2460  const refinementFeatures& features = meshRefiner_.features();
2461 
2462  // Find nearest mesh point to feature edge
2463  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2464  // Reverse lookup : go through all edgeAttractors and find the
2465  // nearest point on pp
2466 
2467  // Get search domain and extend it a bit
2468  treeBoundBox bb(pp.localPoints());
2469  {
2470  // Random number generator. Bit dodgy since not exactly random ;-)
2471  Random rndGen(65431);
2472 
2473  // Slightly extended bb. Slightly off-centred just so on symmetric
2474  // geometry there are less face/edge aligned items.
2475  bb = bb.extend(rndGen, 1e-4);
2476  bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
2477  bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
2478  }
2479 
2480  // Collect candidate points for attraction
2481  DynamicList<label> attractPoints(pp.nPoints());
2482  {
2483  const fvMesh& mesh = meshRefiner_.mesh();
2484 
2485  boolList isFeatureEdgeOrPoint(pp.nPoints(), false);
2486  label nFeats = 0;
2487  forAll(rawPatchConstraints, pointi)
2488  {
2489  if (rawPatchConstraints[pointi].first() >= 2)
2490  {
2491  isFeatureEdgeOrPoint[pointi] = true;
2492  nFeats++;
2493  }
2494  }
2495 
2496  Info<< "Initially selected " << returnReduce(nFeats, sumOp<label>())
2497  << " points out of " << returnReduce(pp.nPoints(), sumOp<label>())
2498  << " for reverse attraction." << endl;
2499 
2500  // Make sure is synchronised (note: check if constraint is already
2501  // synced in which case this is not needed here)
2503  (
2504  mesh,
2505  pp.meshPoints(),
2506  isFeatureEdgeOrPoint,
2507  orEqOp<bool>(), // combine op
2508  false
2509  );
2510 
2511  for (label nGrow = 0; nGrow < 1; nGrow++)
2512  {
2513  boolList newIsFeatureEdgeOrPoint(isFeatureEdgeOrPoint);
2514 
2515  forAll(pp.localFaces(), facei)
2516  {
2517  const face& f = pp.localFaces()[facei];
2518 
2519  forAll(f, fp)
2520  {
2521  if (isFeatureEdgeOrPoint[f[fp]])
2522  {
2523  // Mark all points on face
2524  forAll(f, fp)
2525  {
2526  newIsFeatureEdgeOrPoint[f[fp]] = true;
2527  }
2528  break;
2529  }
2530  }
2531  }
2532 
2533  isFeatureEdgeOrPoint = newIsFeatureEdgeOrPoint;
2534 
2536  (
2537  mesh,
2538  pp.meshPoints(),
2539  isFeatureEdgeOrPoint,
2540  orEqOp<bool>(), // combine op
2541  false
2542  );
2543  }
2544 
2545 
2546  // Collect attractPoints
2547  forAll(isFeatureEdgeOrPoint, pointi)
2548  {
2549  if (isFeatureEdgeOrPoint[pointi])
2550  {
2551  attractPoints.append(pointi);
2552  }
2553  }
2554 
2555  Info<< "Selected "
2556  << returnReduce(attractPoints.size(), sumOp<label>())
2557  << " points out of " << returnReduce(pp.nPoints(), sumOp<label>())
2558  << " for reverse attraction." << endl;
2559  }
2560 
2561 
2563  (
2564  treeDataPoint(pp.localPoints(), attractPoints),
2565  bb, // overall search domain
2566  8, // maxLevel
2567  10, // leafsize
2568  3.0 // duplicity
2569  );
2570 
2571  // Per mesh point the point on nearest feature edge.
2572  patchAttraction.setSize(pp.nPoints());
2573  patchAttraction = Zero;
2574  patchConstraints.setSize(pp.nPoints());
2575  patchConstraints = pointConstraint();
2576 
2577  forAll(edgeAttractors, featI)
2578  {
2579  const List<DynamicList<point>>& edgeAttr = edgeAttractors[featI];
2580  const List<DynamicList<pointConstraint>>& edgeConstr =
2581  edgeConstraints[featI];
2582 
2583  forAll(edgeAttr, featEdgeI)
2584  {
2585  const DynamicList<point>& attr = edgeAttr[featEdgeI];
2586  forAll(attr, i)
2587  {
2588  // Find nearest pp point
2589  const point& featPt = attr[i];
2590  pointIndexHit nearInfo = ppTree.findNearest
2591  (
2592  featPt,
2593  sqr(GREAT)
2594  );
2595 
2596  if (nearInfo.hit())
2597  {
2598  label pointi =
2599  ppTree.shapes().pointLabels()[nearInfo.index()];
2600  const point attraction = featPt-pp.localPoints()[pointi];
2601 
2602  // Check if this point is already being attracted. If so
2603  // override it only if nearer.
2604  if
2605  (
2606  patchConstraints[pointi].first() <= 1
2607  || magSqr(attraction) < magSqr(patchAttraction[pointi])
2608  )
2609  {
2610  patchAttraction[pointi] = attraction;
2611  patchConstraints[pointi] = edgeConstr[featEdgeI][i];
2612  }
2613  }
2614  else
2615  {
2617  << "Did not find pp point near " << featPt
2618  << endl;
2619  }
2620  }
2621  }
2622  }
2623 
2624 
2625  // Different procs might have different patchAttraction,patchConstraints
2626  // however these only contain geometric information, no topology
2627  // so as long as we synchronise after overriding with feature points
2628  // there is no problem, just possibly a small error.
2629 
2630 
2631  // Find nearest mesh point to feature point
2632  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2633  // (overrides attraction to feature edge)
2634  forAll(pointAttractor, featI)
2635  {
2636  const labelList& pointAttr = pointAttractor[featI];
2637  const List<pointConstraint>& pointConstr = pointConstraints[featI];
2638 
2639  forAll(pointAttr, featPointi)
2640  {
2641  if (pointAttr[featPointi] != -1)
2642  {
2643  const point& featPt = features[featI].points()
2644  [
2645  featPointi
2646  ];
2647 
2648  // Find nearest pp point
2649  pointIndexHit nearInfo = ppTree.findNearest
2650  (
2651  featPt,
2652  sqr(GREAT)
2653  );
2654 
2655  if (nearInfo.hit())
2656  {
2657  label pointi =
2658  ppTree.shapes().pointLabels()[nearInfo.index()];
2659 
2660  const point& pt = pp.localPoints()[pointi];
2661  const point attraction = featPt-pt;
2662 
2663  // - already attracted to feature edge : point always wins
2664  // - already attracted to feature point: nearest wins
2665 
2666  if (patchConstraints[pointi].first() <= 1)
2667  {
2668  patchAttraction[pointi] = attraction;
2669  patchConstraints[pointi] = pointConstr[featPointi];
2670  }
2671  else if (patchConstraints[pointi].first() == 2)
2672  {
2673  patchAttraction[pointi] = attraction;
2674  patchConstraints[pointi] = pointConstr[featPointi];
2675  }
2676  else if (patchConstraints[pointi].first() == 3)
2677  {
2678  // Only if nearer
2679  if
2680  (
2681  magSqr(attraction)
2682  < magSqr(patchAttraction[pointi])
2683  )
2684  {
2685  patchAttraction[pointi] = attraction;
2686  patchConstraints[pointi] =
2687  pointConstr[featPointi];
2688  }
2689  }
2690  }
2691  }
2692  }
2693  }
2694 }
2695 
2696 
2697 void Foam::snappySnapDriver::featureAttractionUsingFeatureEdges
2698 (
2699  const label iter,
2700  const bool avoidSnapProblems,
2701  const scalar featureCos,
2702  const bool multiRegionFeatureSnap,
2703 
2704  const indirectPrimitivePatch& pp,
2705  const scalarField& snapDist,
2706  const vectorField& nearestDisp,
2707 
2708  const List<List<point>>& pointFaceSurfNormals,
2709  const List<List<point>>& pointFaceDisp,
2710  const List<List<point>>& pointFaceCentres,
2711  const labelListList& pointFacePatchID,
2712 
2713  vectorField& patchAttraction,
2714  List<pointConstraint>& patchConstraints
2715 ) const
2716 {
2717  const refinementFeatures& features = meshRefiner_.features();
2718 
2719  // Collect ordered attractions on feature edges
2720  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2721 
2722  // Per feature, per feature-edge a list of attraction points and their
2723  // originating vertex.
2724  List<List<DynamicList<point>>> edgeAttractors(features.size());
2725  List<List<DynamicList<pointConstraint>>> edgeConstraints
2726  (
2727  features.size()
2728  );
2729  forAll(features, featI)
2730  {
2731  label nFeatEdges = features[featI].edges().size();
2732  edgeAttractors[featI].setSize(nFeatEdges);
2733  edgeConstraints[featI].setSize(nFeatEdges);
2734  }
2735 
2736  // Per feature, per feature-point the pp point that is attracted to it.
2737  // This list is only used to subset the feature-points that are actually
2738  // used.
2739  List<labelList> pointAttractor(features.size());
2741  forAll(features, featI)
2742  {
2743  label nFeatPoints = features[featI].points().size();
2744  pointAttractor[featI].setSize(nFeatPoints, -1);
2745  pointConstraints[featI].setSize(nFeatPoints);
2746  }
2747 
2748  // Reverse: from pp point to nearest feature
2749  vectorField rawPatchAttraction(pp.nPoints(), Zero);
2750  List<pointConstraint> rawPatchConstraints(pp.nPoints());
2751 
2752  determineFeatures
2753  (
2754  iter,
2755  featureCos,
2756  multiRegionFeatureSnap,
2757 
2758  pp,
2759  snapDist, // per point max distance and nearest surface
2760  nearestDisp,
2761 
2762  pointFaceSurfNormals, // per face nearest surface
2763  pointFaceDisp,
2764  pointFaceCentres,
2765  pointFacePatchID,
2766 
2767  // Feature-point to pp point
2768  pointAttractor,
2770  // Feature-edge to pp point
2771  edgeAttractors,
2772  edgeConstraints,
2773  // pp point to nearest feature
2774  rawPatchAttraction,
2775  rawPatchConstraints
2776  );
2777 
2778 
2779 
2780  // Baffle handling
2781  // ~~~~~~~~~~~~~~~
2782  // Override pointAttractor, edgeAttractor, rawPatchAttration etc. to
2783  // implement 'baffle' handling.
2784  // Baffle: the mesh pp point originates from a loose standing
2785  // baffle.
2786  // Sampling the surface with the surrounding face-centres only picks up
2787  // a single triangle normal so above determineFeatures will not have
2788  // detected anything. So explicitly pick up feature edges on the pp
2789  // (after duplicating points & smoothing so will already have been
2790  // expanded) and match these to the features.
2791  determineBaffleFeatures
2792  (
2793  iter,
2794  featureCos,
2795 
2796  pp,
2797  snapDist,
2798 
2799  // Feature-point to pp point
2800  pointAttractor,
2802  // Feature-edge to pp point
2803  edgeAttractors,
2804  edgeConstraints,
2805  // pp point to nearest feature
2806  rawPatchAttraction,
2807  rawPatchConstraints
2808  );
2809 
2810 
2811 
2812  // Reverse lookup: Find nearest mesh point to feature edge
2813  // ~~~~~~~~~~~~~~~~----------------~~~~~~~~~~~~~~~~~~~~~~~
2814  // go through all edgeAttractors and find the nearest point on pp
2815 
2816  reverseAttractMeshPoints
2817  (
2818  iter,
2819 
2820  pp,
2821  snapDist,
2822 
2823  // Feature-point to pp point
2824  pointAttractor,
2826  // Feature-edge to pp point
2827  edgeAttractors,
2828  edgeConstraints,
2829 
2830  // Estimated feature point
2831  rawPatchAttraction,
2832  rawPatchConstraints,
2833 
2834  // pp point to nearest feature
2835  patchAttraction,
2836  patchConstraints
2837  );
2838 
2839 
2840  // Dump
2841  if (debug&meshRefinement::ATTRACTION)
2842  {
2843  OBJstream featureEdgeStr
2844  (
2845  meshRefiner_.mesh().time().path()
2846  / "edgeAttractors_" + name(iter) + ".obj"
2847  );
2848  Info<< "Dumping feature-edge attraction to "
2849  << featureEdgeStr.name() << endl;
2850 
2851  OBJstream featurePointStr
2852  (
2853  meshRefiner_.mesh().time().path()
2854  / "pointAttractors_" + name(iter) + ".obj"
2855  );
2856  Info<< "Dumping feature-point attraction to "
2857  << featurePointStr.name() << endl;
2858 
2859  forAll(patchConstraints, pointi)
2860  {
2861  const point& pt = pp.localPoints()[pointi];
2862  const vector& attr = patchAttraction[pointi];
2863 
2864  if (patchConstraints[pointi].first() == 2)
2865  {
2866  featureEdgeStr.write(linePointRef(pt, pt+attr));
2867  }
2868  else if (patchConstraints[pointi].first() == 3)
2869  {
2870  featurePointStr.write(linePointRef(pt, pt+attr));
2871  }
2872  }
2873  }
2874 
2875 
2876  if (avoidSnapProblems)
2877  {
2878  //MEJ: any faces that have multi-patch points only keep the multi-patch
2879  // points. The other points on the face will be dragged along
2880  // (hopefully)
2881  if (multiRegionFeatureSnap)
2882  {
2883  releasePointsNextToMultiPatch
2884  (
2885  iter,
2886  featureCos,
2887 
2888  pp,
2889  snapDist,
2890 
2891  pointFaceCentres,
2892  pointFacePatchID,
2893 
2894  rawPatchAttraction,
2895  rawPatchConstraints,
2896 
2897  patchAttraction,
2898  patchConstraints
2899  );
2900  }
2901 
2902 
2903  // Snap edges to feature edges
2904  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
2905  // Walk existing edges and snap remaining ones (that are marked as
2906  // feature edges in rawPatchConstraints)
2907 
2908  stringFeatureEdges
2909  (
2910  iter,
2911  featureCos,
2912 
2913  pp,
2914  snapDist,
2915 
2916  rawPatchAttraction,
2917  rawPatchConstraints,
2918 
2919  patchAttraction,
2920  patchConstraints
2921  );
2922 
2923 
2924  // Avoid diagonal attraction
2925  // ~~~~~~~~~~~~~~~~~~~~~~~~~
2926  // Attract one of the non-diagonal points.
2927 
2928  avoidDiagonalAttraction
2929  (
2930  iter,
2931  featureCos,
2932  pp,
2933  patchAttraction,
2934  patchConstraints
2935  );
2936  }
2937 
2938  if (debug&meshRefinement::ATTRACTION)
2939  {
2940  dumpMove
2941  (
2942  meshRefiner_.mesh().time().path()
2943  / "patchAttraction_" + name(iter) + ".obj",
2944  pp.localPoints(),
2945  pp.localPoints() + patchAttraction
2946  );
2947  }
2948 }
2949 
2950 
2951 void Foam::snappySnapDriver::preventFaceSqueeze
2952 (
2953  const label iter,
2954  const scalar featureCos,
2955 
2956  const indirectPrimitivePatch& pp,
2957  const scalarField& snapDist,
2958 
2959  vectorField& patchAttraction,
2960  List<pointConstraint>& patchConstraints
2961 ) const
2962 {
2964  face singleF;
2965  forAll(pp.localFaces(), facei)
2966  {
2967  const face& f = pp.localFaces()[facei];
2968 
2969  if (f.size() != points.size())
2970  {
2971  points.setSize(f.size());
2972  singleF.setSize(f.size());
2973  for (label i = 0; i < f.size(); i++)
2974  {
2975  singleF[i] = i;
2976  }
2977  }
2978  label nConstraints = 0;
2979  forAll(f, fp)
2980  {
2981  label pointi = f[fp];
2982  const point& pt = pp.localPoints()[pointi];
2983 
2984  if (patchConstraints[pointi].first() > 1)
2985  {
2986  points[fp] = pt + patchAttraction[pointi];
2987  nConstraints++;
2988  }
2989  else
2990  {
2991  points[fp] = pt;
2992  }
2993  }
2994 
2995  if (nConstraints == f.size())
2996  {
2997  scalar oldArea = f.mag(pp.localPoints());
2998  scalar newArea = singleF.mag(points);
2999  if (newArea < 0.1*oldArea)
3000  {
3001  // For now remove the point with largest distance
3002  label maxFp = -1;
3003  scalar maxS = -1;
3004  forAll(f, fp)
3005  {
3006  scalar s = magSqr(patchAttraction[f[fp]]);
3007  if (s > maxS)
3008  {
3009  maxS = s;
3010  maxFp = fp;
3011  }
3012  }
3013  if (maxFp != -1)
3014  {
3015  label pointi = f[maxFp];
3016  // Lower attraction on pointi
3017  patchAttraction[pointi] *= 0.5;
3018  }
3019  }
3020  }
3021  }
3022 }
3023 
3024 
3025 Foam::vectorField Foam::snappySnapDriver::calcNearestSurfaceFeature
3026 (
3027  const snapParameters& snapParams,
3028  const bool avoidSnapProblems,
3029  const label iter,
3030  const scalar featureCos,
3031  const scalar featureAttract,
3032  const scalarField& snapDist,
3033  const vectorField& nearestDisp,
3034  motionSmoother& meshMover,
3035  vectorField& patchAttraction,
3036  List<pointConstraint>& patchConstraints
3037 ) const
3038 {
3039  const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
3040  const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
3041  const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
3042 
3043  Info<< "Overriding displacement on features :" << nl
3044  << " implicit features : " << implicitFeatureAttraction << nl
3045  << " explicit features : " << explicitFeatureAttraction << nl
3046  << " multi-patch features : " << multiRegionFeatureSnap << nl
3047  << endl;
3048 
3049 
3050  const indirectPrimitivePatch& pp = meshMover.patch();
3051  const pointField& localPoints = pp.localPoints();
3052  const fvMesh& mesh = meshRefiner_.mesh();
3053 
3054 
3055 
3056  // Per point, per surrounding face:
3057  // - faceSurfaceNormal
3058  // - faceDisp
3059  // - faceCentres
3060  List<List<point>> pointFaceSurfNormals;
3061  List<List<point>> pointFaceDisp;
3062  List<List<point>> pointFaceCentres;
3063  List<labelList> pointFacePatchID;
3064 
3065  {
3066  // Calculate attraction distance per face (from the attraction distance
3067  // per point)
3068  scalarField faceSnapDist(pp.size(), -GREAT);
3069  forAll(pp.localFaces(), facei)
3070  {
3071  const face& f = pp.localFaces()[facei];
3072  forAll(f, fp)
3073  {
3074  faceSnapDist[facei] = max(faceSnapDist[facei], snapDist[f[fp]]);
3075  }
3076  }
3077 
3078 
3079  // Displacement and orientation per pp face
3080  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3081 
3082  // vector from point on surface back to face centre
3083  vectorField faceDisp(pp.size(), Zero);
3084  // normal of surface at point on surface
3085  vectorField faceSurfaceNormal(pp.size(), Zero);
3086  labelList faceSurfaceGlobalRegion(pp.size(), -1);
3087  vectorField faceRotation(pp.size(), Zero);
3088 
3089  calcNearestFace
3090  (
3091  iter,
3092  pp,
3093  faceSnapDist,
3094  faceDisp,
3095  faceSurfaceNormal,
3096  faceSurfaceGlobalRegion,
3097  faceRotation
3098  );
3099 
3100 
3101  // Collect (possibly remote) per point data of all surrounding faces
3102  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3103  // - faceSurfaceNormal
3104  // - faceDisp
3105  // - faceCentres
3106  calcNearestFacePointProperties
3107  (
3108  iter,
3109  pp,
3110 
3111  faceDisp,
3112  faceSurfaceNormal,
3113  faceSurfaceGlobalRegion,
3114 
3115  pointFaceSurfNormals,
3116  pointFaceDisp,
3117  pointFaceCentres,
3118  pointFacePatchID
3119  );
3120  }
3121 
3122 
3123  // Start off with nearest point on surface
3124  vectorField patchDisp = nearestDisp;
3125 
3126 
3127  // Main calculation
3128  // ~~~~~~~~~~~~~~~~
3129  // This is the main intelligence which calculates per point the vector to
3130  // attract it to the nearest surface. There are lots of possibilities
3131  // here.
3132 
3133  // Nearest feature
3134  patchAttraction.setSize(localPoints.size());
3135  patchAttraction = Zero;
3136  // Constraints at feature
3137  patchConstraints.setSize(localPoints.size());
3138  patchConstraints = pointConstraint();
3139 
3140  if (implicitFeatureAttraction)
3141  {
3142  // Sample faces around each point and see if nearest surface normal
3143  // differs. Reconstruct a feature edge/point if possible and snap to
3144  // it.
3145  featureAttractionUsingReconstruction
3146  (
3147  iter,
3148  avoidSnapProblems,
3149  featureCos,
3150 
3151  pp,
3152  snapDist,
3153  nearestDisp,
3154 
3155  pointFaceSurfNormals,
3156  pointFaceDisp,
3157  pointFaceCentres,
3158  pointFacePatchID,
3159 
3160  patchAttraction,
3161  patchConstraints
3162  );
3163  }
3164 
3165  if (explicitFeatureAttraction)
3166  {
3167  // Sample faces around each point and see if nearest surface normal
3168  // differs. For those find the nearest real feature edge/point and
3169  // store the correspondence. Then loop over feature edge/point
3170  // and attract those nearest mesh point. (the first phase just is
3171  // a subsetting of candidate points, the second makes sure that only
3172  // one mesh point gets attracted per feature)
3173  featureAttractionUsingFeatureEdges
3174  (
3175  iter,
3176  avoidSnapProblems,
3177  featureCos,
3178  multiRegionFeatureSnap,
3179 
3180  pp,
3181  snapDist,
3182  nearestDisp,
3183 
3184  pointFaceSurfNormals,
3185  pointFaceDisp,
3186  pointFaceCentres,
3187  pointFacePatchID,
3188 
3189  patchAttraction,
3190  patchConstraints
3191  );
3192  }
3193 
3194  preventFaceSqueeze
3195  (
3196  iter,
3197  featureCos,
3198 
3199  pp,
3200  snapDist,
3201 
3202  patchAttraction,
3203  patchConstraints
3204  );
3205 
3206  //const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
3207  const PackedBoolList isPatchMasterPoint
3208  (
3210  (
3211  mesh,
3212  pp.meshPoints()
3213  )
3214  );
3215  {
3216  vector avgPatchDisp = meshRefinement::gAverage
3217  (
3218  isPatchMasterPoint,
3219  patchDisp
3220  );
3221  vector avgPatchAttr = meshRefinement::gAverage
3222  (
3223  isPatchMasterPoint,
3224  patchAttraction
3225  );
3226 
3227  Info<< "Attraction:" << endl
3228  << " linear : max:" << gMaxMagSqr(patchDisp)
3229  << " avg:" << avgPatchDisp << endl
3230  << " feature : max:" << gMaxMagSqr(patchAttraction)
3231  << " avg:" << avgPatchAttr << endl;
3232  }
3233 
3234  // So now we have:
3235  // - patchDisp : point movement to go to nearest point on surface
3236  // (either direct or through interpolation of
3237  // face nearest)
3238  // - patchAttraction : direct attraction to features
3239  // - patchConstraints : type of features
3240 
3241  // Use any combination of patchDisp and direct feature attraction.
3242 
3243 
3244  // Mix with direct feature attraction
3245  forAll(patchConstraints, pointi)
3246  {
3247  if (patchConstraints[pointi].first() > 1)
3248  {
3249  patchDisp[pointi] =
3250  (1.0-featureAttract)*patchDisp[pointi]
3251  + featureAttract*patchAttraction[pointi];
3252  }
3253  }
3254 
3255 
3256 
3257  // Count
3258  {
3259  label nMasterPoints = 0;
3260  label nPlanar = 0;
3261  label nEdge = 0;
3262  label nPoint = 0;
3263 
3264  forAll(patchConstraints, pointi)
3265  {
3266  if (isPatchMasterPoint[pointi])
3267  {
3268  nMasterPoints++;
3269 
3270  if (patchConstraints[pointi].first() == 1)
3271  {
3272  nPlanar++;
3273  }
3274  else if (patchConstraints[pointi].first() == 2)
3275  {
3276  nEdge++;
3277  }
3278  else if (patchConstraints[pointi].first() == 3)
3279  {
3280  nPoint++;
3281  }
3282  }
3283  }
3284 
3285  reduce(nMasterPoints, sumOp<label>());
3286  reduce(nPlanar, sumOp<label>());
3287  reduce(nEdge, sumOp<label>());
3288  reduce(nPoint, sumOp<label>());
3289  Info<< "Feature analysis : total master points:"
3290  << nMasterPoints
3291  << " attraction to :" << nl
3292  << " feature point : " << nPoint << nl
3293  << " feature edge : " << nEdge << nl
3294  << " nearest surface : " << nPlanar << nl
3295  << " rest : " << nMasterPoints-nPoint-nEdge-nPlanar
3296  << nl
3297  << endl;
3298  }
3299 
3300 
3301  // Now we have the displacement per patch point to move onto the surface
3302  // Split into tangential and normal direction.
3303  // - start off with all non-constrained points following the constrained
3304  // ones since point normals not relevant.
3305  // - finish with only tangential component smoothed.
3306  // Note: tangential is most
3307  // likely to come purely from face-centre snapping, not face rotation.
3308  // Note: could use the constraints here (constraintTransformation())
3309  // but this is not necessarily accurate and we're smoothing to
3310  // get out of problems.
3311 
3312  if (featureAttract < 1-0.001)
3313  {
3314  //const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
3315  const labelList meshEdges
3316  (
3317  pp.meshEdges(mesh.edges(), mesh.pointEdges())
3318  );
3319  const PackedBoolList isPatchMasterEdge
3320  (
3322  (
3323  mesh,
3324  meshEdges
3325  )
3326  );
3327 
3328  const vectorField pointNormals
3329  (
3331  (
3332  mesh,
3333  pp
3334  )
3335  );
3336 
3337 
3338  // 1. Smoothed all displacement
3339  vectorField smoothedPatchDisp = patchDisp;
3340  smoothAndConstrain
3341  (
3342  isPatchMasterEdge,
3343  pp,
3344  meshEdges,
3345  patchConstraints,
3346  smoothedPatchDisp
3347  );
3348 
3349 
3350  // 2. Smoothed tangential component
3351  vectorField tangPatchDisp = patchDisp;
3352  tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
3353  smoothAndConstrain
3354  (
3355  isPatchMasterEdge,
3356  pp,
3357  meshEdges,
3358  patchConstraints,
3359  tangPatchDisp
3360  );
3361 
3362  // Re-add normal component
3363  tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
3364 
3365  if (debug&meshRefinement::ATTRACTION)
3366  {
3367  dumpMove
3368  (
3369  mesh.time().path()
3370  / "tangPatchDispConstrained_" + name(iter) + ".obj",
3371  pp.localPoints(),
3372  pp.localPoints() + tangPatchDisp
3373  );
3374  }
3375 
3376  patchDisp =
3377  (1.0-featureAttract)*smoothedPatchDisp
3378  + featureAttract*tangPatchDisp;
3379  }
3380 
3381 
3382  const scalar relax = featureAttract;
3383  patchDisp *= relax;
3384 
3385 
3386  // Points on zones in one domain but only present as point on other
3387  // will not do condition 2 on all. Sync explicitly.
3389  (
3390  mesh,
3391  pp.meshPoints(),
3392  patchDisp,
3393  minMagSqrEqOp<point>(), // combine op
3394  vector(GREAT, GREAT, GREAT) // null value (note: cant use VGREAT)
3395  );
3396 
3397  return patchDisp;
3398 }
3399 
3400 
3401 // ************************************************************************* //
const labelListList & pointFaces() const
Given a displacement moves the mesh by scaling the displacement back until there are no more mesh err...
cachedRandom rndGen(label(0),-1)
const labelList & patchID() const
Per boundary face label the patch index.
label nPoints() const
Return number of points supporting patch faces.
Switch multiRegionFeatureSnap() const
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
Holds (reference to) pointField. Encapsulation of data needed for octree searches. Used for searching for nearest point. No bounding boxes around points. Only overlaps and calcNearest are implemented, rest makes little sense.
Definition: treeDataPoint.H:59
A direction and a reference point.
Definition: plane.H:73
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const Type1 & first() const
Return first.
Definition: Tuple2.H:94
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
UEqn relax()
const double e
Elementary charge.
Definition: doubleFloat.H:78
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
A 2-tuple for storing two objects of different types.
Definition: Tuple2.H:47
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:313
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
Switch explicitFeatureSnap() const
bool hit() const
Is there a hit.
const vectorField & faceAreas() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
const labelList & boundaryPoints() const
Return list of boundary points,.
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
Definition: Switch.H:60
label otherVertex(const label a) const
Given one vertex, return the other.
Definition: edgeI.H:103
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const indirectPrimitivePatch & patch() const
Reference to patch.
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
const vectorField & faceCentres() const
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
T & first()
Return the first element of the list.
Definition: UListI.H:114
void findNearestRegion(const labelList &surfacesToTest, const pointField &samples, const scalarField &nearestDistSqr, labelList &hitSurface, labelList &hitRegion) const
Find nearest point on surfaces. Return surface and region on.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
static PackedBoolList getMasterFaces(const polyMesh &)
Get per face whether it is uncoupled or a master of a.
Definition: syncTools.C:153
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
const labelListList & pointEdges() const
Return point-edge addressing.
Container for data on surfaces used for surface-driven refinement. Contains all the data about the le...
void operator()(List< T > &x, const List< T > &y) const
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
Application of (multi-)patch point contraints.
Type gMaxMagSqr(const UList< Type > &f, const label comm)
Encapsulates queries for features.
void applyConstraint(const vector &cd)
Apply and accumulate the effect of the given constraint direction.
scalar y
A list of faces which address into the list of points.
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
label findZoneID(const word &zoneName) const
Find zone index given a name.
Definition: ZoneMesh.C:341
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().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))
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
ray planeIntersect(const plane &) const
Return the cutting line between this plane and another.
Definition: plane.C:323
void reset(T *=0)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
dynamicFvMesh & mesh
Default transformation behaviour.
dimensionedScalar cos(const dimensionedScalar &ds)
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
const pointField & points
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:310
A class for handling words, derived from string.
Definition: word.H:59
static PackedBoolList getMasterEdges(const polyMesh &mesh, const labelList &meshEdges)
Determine master edge for subset of edges. If coupled.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Simple container to keep together snap specific information.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:97
void findNearestPoint(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo) const
Find nearest feature point. Sets:
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints,-1);const cellModel &hex=*(cellModeller::lookup("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]=cellShape(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
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
static const label labelMax
Definition: label.H:62
List< label > labelList
A List of labels.
Definition: labelList.H:56
const labelListList & pointEdges() const
static const zero Zero
Definition: zero.H:91
Accumulates point constraints through successive applications of the applyConstraint function...
point planePlaneIntersect(const plane &, const plane &) const
Return the cutting point between this plane and two other planes.
Definition: plane.C:391
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
Default transformation behaviour for position.
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set).
Definition: autoPtrI.H:83
Simple random number generator.
Definition: Random.H:49
static tmp< pointField > pointNormals(const polyMesh &, const PrimitivePatch< Face, FaceList, PointField, PointType > &)
Return parallel consistent point normals for patches using mesh points.
Foam::polyBoundaryMesh.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
OFstream which keeps track of vertices.
Definition: OBJstream.H:53
static const char nl
Definition: Ostream.H:262
const Field< PointType > & faceCentres() const
Return face centres for patch.
const Type2 & second() const
Return second.
Definition: Tuple2.H:106
labelList f(nPoints)
const Point & missPoint() const
Return miss point.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
label nEdges() const
Return number of edges in patch.
label globalRegion(const label surfI, const label regionI) const
From surface and region on surface to global region.
label nFaces() const
void setSize(const label)
Reset size of List.
Definition: List.C:295
A bit-packed bool list.
const PtrList< surfaceZonesInfo > & surfZones() const
static T gAverage(const PackedBoolList &isMasterElem, const UList< T > &values)
Helper: calculate average.
label patchi
const Field< PointType > & faceNormals() const
Return face normals for patch.
vector point
Point is a vector.
Definition: point.H:41
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
#define WarningInFunction
Report a warning using Foam::Warning.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:62
const Type & shapes() const
Reference to shape.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
label index() const
Return index.
A List with indirect addressing.
Definition: fvMatrix.H:106
const dimensionedScalar c
Speed of light in a vacuum.
const Point & hitPoint() const
Return hit point.
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
const labelList & surfaces() const
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:463
const labelListList & edgeFaces() const
Return edge-face addressing.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:53
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
const labelListList & pointFaces() const
Return point-face addressing.
volScalarField & p
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
static PackedBoolList getMasterPoints(const polyMesh &mesh, const labelList &meshPoints)
Determine master point for subset of points. If coupled.
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:224
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
fileName path() const
Return path.
Definition: Time.H:269
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
Return labels of patch edges in the global edge list using.
label nInternalFaces() const
A List with indirect addressing.
Definition: IndirectList.H:102
scalar mag(const pointField &) const
Magnitude of face area.
Definition: faceI.H:97
const Field< PointType > & localPoints() const
Return pointField of points in patch.
void findNearestEdge(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo, vectorField &nearNormal) const
Find nearest point on nearest feature edge. Sets:
void findNearestRegionEdge(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo, vectorField &nearNormal) const
Find nearest point on nearest region edge. Sets:
Tuple2< scalar, label > nearInfo
Private class for finding nearest.
Switch implicitFeatureSnap() const
Namespace for OpenFOAM.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:465
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29