snappySnapDriverFeature.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 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 occurrence
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.0));
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  bb = bb.extend(1e-4);
2470 
2471  // Collect candidate points for attraction
2472  DynamicList<label> attractPoints(pp.nPoints());
2473  {
2474  const fvMesh& mesh = meshRefiner_.mesh();
2475 
2476  boolList isFeatureEdgeOrPoint(pp.nPoints(), false);
2477  label nFeats = 0;
2478  forAll(rawPatchConstraints, pointi)
2479  {
2480  if (rawPatchConstraints[pointi].first() >= 2)
2481  {
2482  isFeatureEdgeOrPoint[pointi] = true;
2483  nFeats++;
2484  }
2485  }
2486 
2487  Info<< "Initially selected " << returnReduce(nFeats, sumOp<label>())
2488  << " points out of " << returnReduce(pp.nPoints(), sumOp<label>())
2489  << " for reverse attraction." << endl;
2490 
2491  // Make sure is synchronised (note: check if constraint is already
2492  // synced in which case this is not needed here)
2494  (
2495  mesh,
2496  pp.meshPoints(),
2497  isFeatureEdgeOrPoint,
2498  orEqOp<bool>(), // combine op
2499  false
2500  );
2501 
2502  for (label nGrow = 0; nGrow < 1; nGrow++)
2503  {
2504  boolList newIsFeatureEdgeOrPoint(isFeatureEdgeOrPoint);
2505 
2506  forAll(pp.localFaces(), facei)
2507  {
2508  const face& f = pp.localFaces()[facei];
2509 
2510  forAll(f, fp)
2511  {
2512  if (isFeatureEdgeOrPoint[f[fp]])
2513  {
2514  // Mark all points on face
2515  forAll(f, fp)
2516  {
2517  newIsFeatureEdgeOrPoint[f[fp]] = true;
2518  }
2519  break;
2520  }
2521  }
2522  }
2523 
2524  isFeatureEdgeOrPoint = newIsFeatureEdgeOrPoint;
2525 
2527  (
2528  mesh,
2529  pp.meshPoints(),
2530  isFeatureEdgeOrPoint,
2531  orEqOp<bool>(), // combine op
2532  false
2533  );
2534  }
2535 
2536 
2537  // Collect attractPoints
2538  forAll(isFeatureEdgeOrPoint, pointi)
2539  {
2540  if (isFeatureEdgeOrPoint[pointi])
2541  {
2542  attractPoints.append(pointi);
2543  }
2544  }
2545 
2546  Info<< "Selected "
2547  << returnReduce(attractPoints.size(), sumOp<label>())
2548  << " points out of " << returnReduce(pp.nPoints(), sumOp<label>())
2549  << " for reverse attraction." << endl;
2550  }
2551 
2552 
2554  (
2555  treeDataPoint(pp.localPoints(), attractPoints),
2556  bb, // overall search domain
2557  8, // maxLevel
2558  10, // leafsize
2559  3.0 // duplicity
2560  );
2561 
2562  // Per mesh point the point on nearest feature edge.
2563  patchAttraction.setSize(pp.nPoints());
2564  patchAttraction = Zero;
2565  patchConstraints.setSize(pp.nPoints());
2566  patchConstraints = pointConstraint();
2567 
2568  forAll(edgeAttractors, feati)
2569  {
2570  const List<DynamicList<point>>& edgeAttr = edgeAttractors[feati];
2571  const List<DynamicList<pointConstraint>>& edgeConstr =
2572  edgeConstraints[feati];
2573 
2574  forAll(edgeAttr, featEdgei)
2575  {
2576  const DynamicList<point>& attr = edgeAttr[featEdgei];
2577  forAll(attr, i)
2578  {
2579  // Find nearest pp point
2580  const point& featPt = attr[i];
2581  pointIndexHit nearInfo = ppTree.findNearest
2582  (
2583  featPt,
2584  sqr(great)
2585  );
2586 
2587  if (nearInfo.hit())
2588  {
2589  label pointi =
2590  ppTree.shapes().pointLabels()[nearInfo.index()];
2591  const point attraction = featPt-pp.localPoints()[pointi];
2592 
2593  // Check if this point is already being attracted. If so
2594  // override it only if nearer.
2595  if
2596  (
2597  patchConstraints[pointi].first() <= 1
2598  || magSqr(attraction) < magSqr(patchAttraction[pointi])
2599  )
2600  {
2601  patchAttraction[pointi] = attraction;
2602  patchConstraints[pointi] = edgeConstr[featEdgei][i];
2603  }
2604  }
2605  else
2606  {
2608  << "Did not find pp point near " << featPt
2609  << endl;
2610  }
2611  }
2612  }
2613  }
2614 
2615 
2616  // Different procs might have different patchAttraction,patchConstraints
2617  // however these only contain geometric information, no topology
2618  // so as long as we synchronise after overriding with feature points
2619  // there is no problem, just possibly a small error.
2620 
2621 
2622  // Find nearest mesh point to feature point
2623  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2624  // (overrides attraction to feature edge)
2625  forAll(pointAttractor, feati)
2626  {
2627  const labelList& pointAttr = pointAttractor[feati];
2628  const List<pointConstraint>& pointConstr = pointConstraints[feati];
2629 
2630  forAll(pointAttr, featPointi)
2631  {
2632  if (pointAttr[featPointi] != -1)
2633  {
2634  const point& featPt = features[feati].points()
2635  [
2636  featPointi
2637  ];
2638 
2639  // Find nearest pp point
2640  pointIndexHit nearInfo = ppTree.findNearest
2641  (
2642  featPt,
2643  sqr(great)
2644  );
2645 
2646  if (nearInfo.hit())
2647  {
2648  label pointi =
2649  ppTree.shapes().pointLabels()[nearInfo.index()];
2650 
2651  const point& pt = pp.localPoints()[pointi];
2652  const point attraction = featPt-pt;
2653 
2654  // - already attracted to feature edge : point always wins
2655  // - already attracted to feature point: nearest wins
2656 
2657  if (patchConstraints[pointi].first() <= 1)
2658  {
2659  patchAttraction[pointi] = attraction;
2660  patchConstraints[pointi] = pointConstr[featPointi];
2661  }
2662  else if (patchConstraints[pointi].first() == 2)
2663  {
2664  patchAttraction[pointi] = attraction;
2665  patchConstraints[pointi] = pointConstr[featPointi];
2666  }
2667  else if (patchConstraints[pointi].first() == 3)
2668  {
2669  // Only if nearer
2670  if
2671  (
2672  magSqr(attraction)
2673  < magSqr(patchAttraction[pointi])
2674  )
2675  {
2676  patchAttraction[pointi] = attraction;
2677  patchConstraints[pointi] =
2678  pointConstr[featPointi];
2679  }
2680  }
2681  }
2682  }
2683  }
2684  }
2685 }
2686 
2687 
2688 void Foam::snappySnapDriver::featureAttractionUsingFeatureEdges
2689 (
2690  const label iter,
2691  const bool avoidSnapProblems,
2692  const scalar featureCos,
2693  const bool multiRegionFeatureSnap,
2694 
2695  const indirectPrimitivePatch& pp,
2696  const scalarField& snapDist,
2697  const vectorField& nearestDisp,
2698 
2699  const List<List<point>>& pointFaceSurfNormals,
2700  const List<List<point>>& pointFaceDisp,
2701  const List<List<point>>& pointFaceCentres,
2702  const labelListList& pointFacePatchID,
2703 
2704  vectorField& patchAttraction,
2705  List<pointConstraint>& patchConstraints
2706 ) const
2707 {
2708  const refinementFeatures& features = meshRefiner_.features();
2709 
2710  // Collect ordered attractions on feature edges
2711  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2712 
2713  // Per feature, per feature-edge a list of attraction points and their
2714  // originating vertex.
2715  List<List<DynamicList<point>>> edgeAttractors(features.size());
2716  List<List<DynamicList<pointConstraint>>> edgeConstraints
2717  (
2718  features.size()
2719  );
2720  forAll(features, feati)
2721  {
2722  label nFeatEdges = features[feati].edges().size();
2723  edgeAttractors[feati].setSize(nFeatEdges);
2724  edgeConstraints[feati].setSize(nFeatEdges);
2725  }
2726 
2727  // Per feature, per feature-point the pp point that is attracted to it.
2728  // This list is only used to subset the feature-points that are actually
2729  // used.
2730  List<labelList> pointAttractor(features.size());
2732  forAll(features, feati)
2733  {
2734  label nFeatPoints = features[feati].points().size();
2735  pointAttractor[feati].setSize(nFeatPoints, -1);
2736  pointConstraints[feati].setSize(nFeatPoints);
2737  }
2738 
2739  // Reverse: from pp point to nearest feature
2740  vectorField rawPatchAttraction(pp.nPoints(), Zero);
2741  List<pointConstraint> rawPatchConstraints(pp.nPoints());
2742 
2743  determineFeatures
2744  (
2745  iter,
2746  featureCos,
2747  multiRegionFeatureSnap,
2748 
2749  pp,
2750  snapDist, // per point max distance and nearest surface
2751  nearestDisp,
2752 
2753  pointFaceSurfNormals, // per face nearest surface
2754  pointFaceDisp,
2755  pointFaceCentres,
2756  pointFacePatchID,
2757 
2758  // Feature-point to pp point
2759  pointAttractor,
2761  // Feature-edge to pp point
2762  edgeAttractors,
2763  edgeConstraints,
2764  // pp point to nearest feature
2765  rawPatchAttraction,
2766  rawPatchConstraints
2767  );
2768 
2769 
2770 
2771  // Baffle handling
2772  // ~~~~~~~~~~~~~~~
2773  // Override pointAttractor, edgeAttractor, rawPatchAttration etc. to
2774  // implement 'baffle' handling.
2775  // Baffle: the mesh pp point originates from a loose standing
2776  // baffle.
2777  // Sampling the surface with the surrounding face-centres only picks up
2778  // a single triangle normal so above determineFeatures will not have
2779  // detected anything. So explicitly pick up feature edges on the pp
2780  // (after duplicating points & smoothing so will already have been
2781  // expanded) and match these to the features.
2782  determineBaffleFeatures
2783  (
2784  iter,
2785  featureCos,
2786 
2787  pp,
2788  snapDist,
2789 
2790  // Feature-point to pp point
2791  pointAttractor,
2793  // Feature-edge to pp point
2794  edgeAttractors,
2795  edgeConstraints,
2796  // pp point to nearest feature
2797  rawPatchAttraction,
2798  rawPatchConstraints
2799  );
2800 
2801 
2802 
2803  // Reverse lookup: Find nearest mesh point to feature edge
2804  // ~~~~~~~~~~~~~~~~----------------~~~~~~~~~~~~~~~~~~~~~~~
2805  // go through all edgeAttractors and find the nearest point on pp
2806 
2807  reverseAttractMeshPoints
2808  (
2809  iter,
2810 
2811  pp,
2812  snapDist,
2813 
2814  // Feature-point to pp point
2815  pointAttractor,
2817  // Feature-edge to pp point
2818  edgeAttractors,
2819  edgeConstraints,
2820 
2821  // Estimated feature point
2822  rawPatchAttraction,
2823  rawPatchConstraints,
2824 
2825  // pp point to nearest feature
2826  patchAttraction,
2827  patchConstraints
2828  );
2829 
2830 
2831  // Dump
2832  if (debug&meshRefinement::ATTRACTION)
2833  {
2834  OBJstream featureEdgeStr
2835  (
2836  meshRefiner_.mesh().time().path()
2837  / "edgeAttractors_" + name(iter) + ".obj"
2838  );
2839  Info<< "Dumping feature-edge attraction to "
2840  << featureEdgeStr.name() << endl;
2841 
2842  OBJstream featurePointStr
2843  (
2844  meshRefiner_.mesh().time().path()
2845  / "pointAttractors_" + name(iter) + ".obj"
2846  );
2847  Info<< "Dumping feature-point attraction to "
2848  << featurePointStr.name() << endl;
2849 
2850  forAll(patchConstraints, pointi)
2851  {
2852  const point& pt = pp.localPoints()[pointi];
2853  const vector& attr = patchAttraction[pointi];
2854 
2855  if (patchConstraints[pointi].first() == 2)
2856  {
2857  featureEdgeStr.write(linePointRef(pt, pt+attr));
2858  }
2859  else if (patchConstraints[pointi].first() == 3)
2860  {
2861  featurePointStr.write(linePointRef(pt, pt+attr));
2862  }
2863  }
2864  }
2865 
2866 
2867  if (avoidSnapProblems)
2868  {
2869  // MEJ: any faces that have multi-patch points only keep the multi-patch
2870  // points. The other points on the face will be dragged along
2871  // (hopefully)
2872  if (multiRegionFeatureSnap)
2873  {
2874  releasePointsNextToMultiPatch
2875  (
2876  iter,
2877  featureCos,
2878 
2879  pp,
2880  snapDist,
2881 
2882  pointFaceCentres,
2883  pointFacePatchID,
2884 
2885  rawPatchAttraction,
2886  rawPatchConstraints,
2887 
2888  patchAttraction,
2889  patchConstraints
2890  );
2891  }
2892 
2893 
2894  // Snap edges to feature edges
2895  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
2896  // Walk existing edges and snap remaining ones (that are marked as
2897  // feature edges in rawPatchConstraints)
2898 
2899  stringFeatureEdges
2900  (
2901  iter,
2902  featureCos,
2903 
2904  pp,
2905  snapDist,
2906 
2907  rawPatchAttraction,
2908  rawPatchConstraints,
2909 
2910  patchAttraction,
2911  patchConstraints
2912  );
2913 
2914 
2915  // Avoid diagonal attraction
2916  // ~~~~~~~~~~~~~~~~~~~~~~~~~
2917  // Attract one of the non-diagonal points.
2918 
2919  avoidDiagonalAttraction
2920  (
2921  iter,
2922  featureCos,
2923  pp,
2924  patchAttraction,
2925  patchConstraints
2926  );
2927  }
2928 
2929  if (debug&meshRefinement::ATTRACTION)
2930  {
2931  dumpMove
2932  (
2933  meshRefiner_.mesh().time().path()
2934  / "patchAttraction_" + name(iter) + ".obj",
2935  pp.localPoints(),
2936  pp.localPoints() + patchAttraction
2937  );
2938  }
2939 }
2940 
2941 
2942 void Foam::snappySnapDriver::preventFaceSqueeze
2943 (
2944  const label iter,
2945  const scalar featureCos,
2946 
2947  const indirectPrimitivePatch& pp,
2948  const scalarField& snapDist,
2949 
2950  vectorField& patchAttraction,
2951  List<pointConstraint>& patchConstraints
2952 ) const
2953 {
2955  face singleF;
2956  forAll(pp.localFaces(), facei)
2957  {
2958  const face& f = pp.localFaces()[facei];
2959 
2960  if (f.size() != points.size())
2961  {
2962  points.setSize(f.size());
2963  singleF.setSize(f.size());
2964  for (label i = 0; i < f.size(); i++)
2965  {
2966  singleF[i] = i;
2967  }
2968  }
2969  label nConstraints = 0;
2970  forAll(f, fp)
2971  {
2972  label pointi = f[fp];
2973  const point& pt = pp.localPoints()[pointi];
2974 
2975  if (patchConstraints[pointi].first() > 1)
2976  {
2977  points[fp] = pt + patchAttraction[pointi];
2978  nConstraints++;
2979  }
2980  else
2981  {
2982  points[fp] = pt;
2983  }
2984  }
2985 
2986  if (nConstraints == f.size())
2987  {
2988  scalar oldArea = f.mag(pp.localPoints());
2989  scalar newArea = singleF.mag(points);
2990  if (newArea < 0.1*oldArea)
2991  {
2992  // For now remove the point with largest distance
2993  label maxFp = -1;
2994  scalar maxS = -1;
2995  forAll(f, fp)
2996  {
2997  scalar s = magSqr(patchAttraction[f[fp]]);
2998  if (s > maxS)
2999  {
3000  maxS = s;
3001  maxFp = fp;
3002  }
3003  }
3004  if (maxFp != -1)
3005  {
3006  label pointi = f[maxFp];
3007  // Lower attraction on pointi
3008  patchAttraction[pointi] *= 0.5;
3009  }
3010  }
3011  }
3012  }
3013 }
3014 
3015 
3016 Foam::vectorField Foam::snappySnapDriver::calcNearestSurfaceFeature
3017 (
3018  const snapParameters& snapParams,
3019  const bool avoidSnapProblems,
3020  const label iter,
3021  const scalar featureCos,
3022  const scalar featureAttract,
3023  const scalarField& snapDist,
3024  const vectorField& nearestDisp,
3025  motionSmoother& meshMover,
3026  vectorField& patchAttraction,
3027  List<pointConstraint>& patchConstraints
3028 ) const
3029 {
3030  const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
3031  const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
3032  const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
3033 
3034  Info<< "Overriding displacement on features :" << nl
3035  << " implicit features : " << implicitFeatureAttraction << nl
3036  << " explicit features : " << explicitFeatureAttraction << nl
3037  << " multi-patch features : " << multiRegionFeatureSnap << nl
3038  << endl;
3039 
3040 
3041  const indirectPrimitivePatch& pp = meshMover.patch();
3042  const pointField& localPoints = pp.localPoints();
3043  const fvMesh& mesh = meshRefiner_.mesh();
3044 
3045 
3046 
3047  // Per point, per surrounding face:
3048  // - faceSurfaceNormal
3049  // - faceDisp
3050  // - faceCentres
3051  List<List<point>> pointFaceSurfNormals;
3052  List<List<point>> pointFaceDisp;
3053  List<List<point>> pointFaceCentres;
3054  List<labelList> pointFacePatchID;
3055 
3056  {
3057  // Calculate attraction distance per face (from the attraction distance
3058  // per point)
3059  scalarField faceSnapDist(pp.size(), -great);
3060  forAll(pp.localFaces(), facei)
3061  {
3062  const face& f = pp.localFaces()[facei];
3063  forAll(f, fp)
3064  {
3065  faceSnapDist[facei] = max(faceSnapDist[facei], snapDist[f[fp]]);
3066  }
3067  }
3068 
3069 
3070  // Displacement and orientation per pp face
3071  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3072 
3073  // vector from point on surface back to face centre
3074  vectorField faceDisp(pp.size(), Zero);
3075  // normal of surface at point on surface
3076  vectorField faceSurfaceNormal(pp.size(), Zero);
3077  labelList faceSurfaceGlobalRegion(pp.size(), -1);
3078  vectorField faceRotation(pp.size(), Zero);
3079 
3080  calcNearestFace
3081  (
3082  iter,
3083  pp,
3084  faceSnapDist,
3085  faceDisp,
3086  faceSurfaceNormal,
3087  faceSurfaceGlobalRegion,
3088  faceRotation
3089  );
3090 
3091 
3092  // Collect (possibly remote) per point data of all surrounding faces
3093  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3094  // - faceSurfaceNormal
3095  // - faceDisp
3096  // - faceCentres
3097  calcNearestFacePointProperties
3098  (
3099  iter,
3100  pp,
3101 
3102  faceDisp,
3103  faceSurfaceNormal,
3104  faceSurfaceGlobalRegion,
3105 
3106  pointFaceSurfNormals,
3107  pointFaceDisp,
3108  pointFaceCentres,
3109  pointFacePatchID
3110  );
3111  }
3112 
3113 
3114  // Start off with nearest point on surface
3115  vectorField patchDisp = nearestDisp;
3116 
3117 
3118  // Main calculation
3119  // ~~~~~~~~~~~~~~~~
3120  // This is the main intelligence which calculates per point the vector to
3121  // attract it to the nearest surface. There are lots of possibilities
3122  // here.
3123 
3124  // Nearest feature
3125  patchAttraction.setSize(localPoints.size());
3126  patchAttraction = Zero;
3127  // Constraints at feature
3128  patchConstraints.setSize(localPoints.size());
3129  patchConstraints = pointConstraint();
3130 
3131  if (implicitFeatureAttraction)
3132  {
3133  // Sample faces around each point and see if nearest surface normal
3134  // differs. Reconstruct a feature edge/point if possible and snap to
3135  // it.
3136  featureAttractionUsingReconstruction
3137  (
3138  iter,
3139  avoidSnapProblems,
3140  featureCos,
3141 
3142  pp,
3143  snapDist,
3144  nearestDisp,
3145 
3146  pointFaceSurfNormals,
3147  pointFaceDisp,
3148  pointFaceCentres,
3149  pointFacePatchID,
3150 
3151  patchAttraction,
3152  patchConstraints
3153  );
3154  }
3155 
3156  if (explicitFeatureAttraction)
3157  {
3158  // Sample faces around each point and see if nearest surface normal
3159  // differs. For those find the nearest real feature edge/point and
3160  // store the correspondence. Then loop over feature edge/point
3161  // and attract those nearest mesh point. (the first phase just is
3162  // a subsetting of candidate points, the second makes sure that only
3163  // one mesh point gets attracted per feature)
3164  featureAttractionUsingFeatureEdges
3165  (
3166  iter,
3167  avoidSnapProblems,
3168  featureCos,
3169  multiRegionFeatureSnap,
3170 
3171  pp,
3172  snapDist,
3173  nearestDisp,
3174 
3175  pointFaceSurfNormals,
3176  pointFaceDisp,
3177  pointFaceCentres,
3178  pointFacePatchID,
3179 
3180  patchAttraction,
3181  patchConstraints
3182  );
3183  }
3184 
3185  preventFaceSqueeze
3186  (
3187  iter,
3188  featureCos,
3189 
3190  pp,
3191  snapDist,
3192 
3193  patchAttraction,
3194  patchConstraints
3195  );
3196 
3197  // const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
3198  const PackedBoolList isPatchMasterPoint
3199  (
3201  (
3202  mesh,
3203  pp.meshPoints()
3204  )
3205  );
3206  {
3207  vector avgPatchDisp = meshRefinement::gAverage
3208  (
3209  isPatchMasterPoint,
3210  patchDisp
3211  );
3212  vector avgPatchAttr = meshRefinement::gAverage
3213  (
3214  isPatchMasterPoint,
3215  patchAttraction
3216  );
3217 
3218  Info<< "Attraction:" << endl
3219  << " linear : max:" << gMaxMagSqr(patchDisp)
3220  << " avg:" << avgPatchDisp << endl
3221  << " feature : max:" << gMaxMagSqr(patchAttraction)
3222  << " avg:" << avgPatchAttr << endl;
3223  }
3224 
3225  // So now we have:
3226  // - patchDisp : point movement to go to nearest point on surface
3227  // (either direct or through interpolation of
3228  // face nearest)
3229  // - patchAttraction : direct attraction to features
3230  // - patchConstraints : type of features
3231 
3232  // Use any combination of patchDisp and direct feature attraction.
3233 
3234 
3235  // Mix with direct feature attraction
3236  forAll(patchConstraints, pointi)
3237  {
3238  if (patchConstraints[pointi].first() > 1)
3239  {
3240  patchDisp[pointi] =
3241  (1.0-featureAttract)*patchDisp[pointi]
3242  + featureAttract*patchAttraction[pointi];
3243  }
3244  }
3245 
3246 
3247 
3248  // Count
3249  {
3250  label nMasterPoints = 0;
3251  label nPlanar = 0;
3252  label nEdge = 0;
3253  label nPoint = 0;
3254 
3255  forAll(patchConstraints, pointi)
3256  {
3257  if (isPatchMasterPoint[pointi])
3258  {
3259  nMasterPoints++;
3260 
3261  if (patchConstraints[pointi].first() == 1)
3262  {
3263  nPlanar++;
3264  }
3265  else if (patchConstraints[pointi].first() == 2)
3266  {
3267  nEdge++;
3268  }
3269  else if (patchConstraints[pointi].first() == 3)
3270  {
3271  nPoint++;
3272  }
3273  }
3274  }
3275 
3276  reduce(nMasterPoints, sumOp<label>());
3277  reduce(nPlanar, sumOp<label>());
3278  reduce(nEdge, sumOp<label>());
3279  reduce(nPoint, sumOp<label>());
3280  Info<< "Feature analysis : total master points:"
3281  << nMasterPoints
3282  << " attraction to :" << nl
3283  << " feature point : " << nPoint << nl
3284  << " feature edge : " << nEdge << nl
3285  << " nearest surface : " << nPlanar << nl
3286  << " rest : " << nMasterPoints-nPoint-nEdge-nPlanar
3287  << nl
3288  << endl;
3289  }
3290 
3291 
3292  // Now we have the displacement per patch point to move onto the surface
3293  // Split into tangential and normal direction.
3294  // - start off with all non-constrained points following the constrained
3295  // ones since point normals not relevant.
3296  // - finish with only tangential component smoothed.
3297  // Note: tangential is most
3298  // likely to come purely from face-centre snapping, not face rotation.
3299  // Note: could use the constraints here (constraintTransformation())
3300  // but this is not necessarily accurate and we're smoothing to
3301  // get out of problems.
3302 
3303  if (featureAttract < 1-0.001)
3304  {
3305  // const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
3306  const labelList meshEdges
3307  (
3308  pp.meshEdges(mesh.edges(), mesh.pointEdges())
3309  );
3310  const PackedBoolList isPatchMasterEdge
3311  (
3313  (
3314  mesh,
3315  meshEdges
3316  )
3317  );
3318 
3319  const vectorField pointNormals
3320  (
3322  (
3323  mesh,
3324  pp
3325  )
3326  );
3327 
3328 
3329  // 1. Smoothed all displacement
3330  vectorField smoothedPatchDisp = patchDisp;
3331  smoothAndConstrain
3332  (
3333  isPatchMasterEdge,
3334  pp,
3335  meshEdges,
3336  patchConstraints,
3337  smoothedPatchDisp
3338  );
3339 
3340 
3341  // 2. Smoothed tangential component
3342  vectorField tangPatchDisp = patchDisp;
3343  tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
3344  smoothAndConstrain
3345  (
3346  isPatchMasterEdge,
3347  pp,
3348  meshEdges,
3349  patchConstraints,
3350  tangPatchDisp
3351  );
3352 
3353  // Re-add normal component
3354  tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
3355 
3356  if (debug&meshRefinement::ATTRACTION)
3357  {
3358  dumpMove
3359  (
3360  mesh.time().path()
3361  / "tangPatchDispConstrained_" + name(iter) + ".obj",
3362  pp.localPoints(),
3363  pp.localPoints() + tangPatchDisp
3364  );
3365  }
3366 
3367  patchDisp =
3368  (1.0-featureAttract)*smoothedPatchDisp
3369  + featureAttract*tangPatchDisp;
3370  }
3371 
3372 
3373  const scalar relax = featureAttract;
3374  patchDisp *= relax;
3375 
3376 
3377  // Points on zones in one domain but only present as point on other
3378  // will not do condition 2 on all. Sync explicitly.
3380  (
3381  mesh,
3382  pp.meshPoints(),
3383  patchDisp,
3384  minMagSqrEqOp<point>(), // combine op
3385  vector(great, great, great) // null value (note: can't use vGreat)
3386  );
3387 
3388  return patchDisp;
3389 }
3390 
3391 
3392 // ************************************************************************* //
Given a displacement moves the mesh by scaling the displacement back until there are no more mesh err...
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
const labelListList & pointEdges() const
Return point-edge addressing.
labelList first(const UList< labelPair > &p)
Definition: patchToPatch.C:38
label nPoints() const
Return number of points supporting patch faces.
label findZoneID(const word &zoneName) const
Find zone index given a name.
Definition: MeshZones.C:341
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:325
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
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
const labelList & surfaces() const
A direction and a reference point.
Definition: plane.H:72
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
void reset(T *=nullptr)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
void findNearestPoint(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo) const
Find nearest feature point. Sets:
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
const Field< PointType > & faceCentres() const
Return face centres for patch.
A 2-tuple for storing two objects of different types.
Definition: HashTable.H:65
label nInternalFaces() const
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
const labelListList & pointEdges() const
const Type1 & first() const
Return first.
Definition: Tuple2.H:116
label nFaces() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
const Point & missPoint() const
Return miss point.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Default transformation behaviour.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
const labelList & boundaryPoints() const
Return list of boundary points,.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const labelList & patchID() const
Per boundary face label the patch index.
ray planeIntersect(const plane &) const
Return the cutting line between this plane and another.
Definition: plane.C:386
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/any.
Definition: Switch.H:60
Switch multiRegionFeatureSnap() const
Tuple2< scalar, label > nearInfo
Private class for finding nearest.
point planePlaneIntersect(const plane &, const plane &) const
Return the cutting point between this plane and two other planes.
Definition: plane.C:454
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
label otherVertex(const label a) const
Given one vertex, return the other.
Definition: edgeI.H:140
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
const Field< PointType > & localPoints() const
Return pointField of points in patch.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:372
fvMesh & mesh
T & first()
Return the first element of the list.
Definition: UListI.H:114
const dimensionedScalar c
Speed of light in a vacuum.
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.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
Switch explicitFeatureSnap() const
Container for data on surfaces used for surface-driven refinement. Contains all the data about the le...
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1211
Application of (multi-)patch point constraints.
Type gMaxMagSqr(const UList< Type > &f, const label comm)
Switch implicitFeatureSnap() const
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.
const Point & hitPoint() const
Return hit point.
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
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
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:316
dimensionedScalar cos(const dimensionedScalar &ds)
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
const pointField & points
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
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.
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)
Synchronise values on all mesh edges.
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:178
Default transformation behaviour for position.
scalar mag(const pointField &) const
Return scalar magnitude.
Definition: faceI.H:81
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
bool hit() const
Is there a hit.
const Type & shapes() const
Reference to shape.
const labelListList & edgeFaces() const
Return edge-face addressing.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
static const label labelMax
Definition: label.H:62
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1224
const Field< PointType > & faceNormals() const
Return face normals for patch.
Accumulates point constraints through successive applications of the applyConstraint function...
const PtrList< surfaceZonesInfo > & surfZones() const
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set)
Definition: autoPtrI.H:83
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
Foam::polyBoundaryMesh.
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, &mergedCyclicPolyPatch::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
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
const labelListList & pointFaces() const
Return point-face addressing.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
label nEdges() const
Return number of edges in patch.
OFstream which keeps track of vertices.
Definition: OBJstream.H:53
static const char nl
Definition: Ostream.H:260
label globalRegion(const label surfi, const label regioni) const
From surface and region on surface to global region.
const indirectPrimitivePatch & patch() const
Reference to patch.
labelList f(nPoints)
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
labelList second(const UList< labelPair > &p)
Definition: patchToPatch.C:48
UEqn relax()
static tmp< pointField > pointNormals(const polyMesh &, const PrimitivePatch< FaceList, PointField > &)
Return parallel consistent point normals for patches using mesh points.
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const vectorField & faceCentres() const
void setSize(const label)
Reset size of List.
Definition: List.C:281
A bit-packed bool list.
static T gAverage(const PackedBoolList &isMasterElem, const UList< T > &values)
Helper: calculate average.
label patchi
void findNearestRegionEdge(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo, vectorField &nearNormal) const
Find nearest point on nearest region edge. Sets:
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:70
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
Return labels of patch edges in the global edge list using.
const meshFaceZones & faceZones() const
Return face zones.
Definition: polyMesh.H:495
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:95
A List with indirect addressing.
Definition: fvMatrix.H:106
const vectorField & faceAreas() const
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:306
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh points.
void findNearestEdge(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo, vectorField &nearNormal) const
Find nearest point on nearest feature edge. Sets:
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
const labelListList & pointFaces() const
const Type2 & second() const
Return second.
Definition: Tuple2.H:128
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
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:52
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:65
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:236
fileName path() const
Return path.
Definition: TimePaths.H:139
void operator()(List< T > &x, const List< T > &y) const
label index() const
Return index.
treeBoundBox extend(const scalar s) const
Return asymetrically extended bounding box, with guaranteed.
Namespace for OpenFOAM.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:483