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