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