motionSmootherAlgo.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-2017 OpenFOAM Foundation
9 Copyright (C) 2018-2022 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 "motionSmootherAlgo.H"
30#include "twoDPointCorrector.H"
31#include "faceSet.H"
32#include "pointSet.H"
34#include "pointConstraints.H"
35#include "syncTools.H"
36#include "meshTools.H"
37#include "OFstream.H"
38
39// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40
41namespace Foam
42{
44}
45
46
47// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48
49void Foam::motionSmootherAlgo::testSyncPositions
50(
51 const pointField& fld,
52 const scalar maxMag
53) const
54{
55 pointField syncedFld(fld);
56
58 (
59 mesh_,
60 syncedFld,
61 minEqOp<point>(), // combine op
62 point(GREAT,GREAT,GREAT) // null
63 );
64
65 forAll(syncedFld, i)
66 {
67 if (mag(syncedFld[i] - fld[i]) > maxMag)
68 {
70 << "On point " << i << " point:" << fld[i]
71 << " synchronised point:" << syncedFld[i]
72 << abort(FatalError);
73 }
74 }
75}
76
77
78void Foam::motionSmootherAlgo::checkFld(const pointScalarField& fld)
79{
80 forAll(fld, pointi)
81 {
82 const scalar val = fld[pointi];
83
84 if ((val > -GREAT) && (val < GREAT))
85 {}
86 else
87 {
89 << "Problem : point:" << pointi << " value:" << val
90 << abort(FatalError);
91 }
92 }
93}
94
95
96Foam::labelHashSet Foam::motionSmootherAlgo::getPoints
97(
98 const labelHashSet& faceLabels
99) const
100{
101 labelHashSet usedPoints(mesh_.nPoints()/100);
102
103 for (const label faceId : faceLabels)
104 {
105 usedPoints.insert(mesh_.faces()[faceId]);
106 }
107
108 return usedPoints;
109}
110
111
112Foam::tmp<Foam::scalarField> Foam::motionSmootherAlgo::calcEdgeWeights
113(
114 const pointField& points
115) const
116{
117 const edgeList& edges = mesh_.edges();
118
119 auto twght = tmp<scalarField>::New(edges.size());
120 auto& wght = twght.ref();
121
122 forAll(edges, edgei)
123 {
124 wght[edgei] = 1.0/(edges[edgei].mag(points)+SMALL);
125 }
126 return twght;
127}
128
129
130// Smooth on selected points (usually patch points)
131void Foam::motionSmootherAlgo::minSmooth
132(
133 const scalarField& edgeWeights,
134 const bitSet& isAffectedPoint,
135 const labelList& meshPoints,
136 const pointScalarField& fld,
137 pointScalarField& newFld
138) const
139{
140 tmp<pointScalarField> tavgFld = avg
141 (
142 fld,
143 edgeWeights //scalarField(mesh_.nEdges(), 1.0) // uniform weighting
144 );
145 const pointScalarField& avgFld = tavgFld();
146
147 for (const label pointi : meshPoints)
148 {
149 if (isAffectedPoint.test(pointi))
150 {
151 newFld[pointi] = min
152 (
153 fld[pointi],
154 0.5*fld[pointi] + 0.5*avgFld[pointi]
155 );
156 }
157 }
158
159 // Single and multi-patch constraints
160 pointConstraints::New(pMesh()).constrain(newFld, false);
161}
162
163
164// Smooth on all internal points
165void Foam::motionSmootherAlgo::minSmooth
166(
167 const scalarField& edgeWeights,
168 const bitSet& isAffectedPoint,
169 const pointScalarField& fld,
170 pointScalarField& newFld
171) const
172{
173 tmp<pointScalarField> tavgFld = avg
174 (
175 fld,
176 edgeWeights //scalarField(mesh_.nEdges(), 1.0) // uniform weighting
177 );
178 const pointScalarField& avgFld = tavgFld();
179
180 forAll(fld, pointi)
181 {
182 if (isAffectedPoint.test(pointi) && isInternalPoint_.test(pointi))
183 {
184 newFld[pointi] = min
185 (
186 fld[pointi],
187 0.5*fld[pointi] + 0.5*avgFld[pointi]
188 );
189 }
190 }
191
192 // Single and multi-patch constraints
193 pointConstraints::New(pMesh()).constrain(newFld, false);
194
195}
196
197
198// Scale on all internal points
199void Foam::motionSmootherAlgo::scaleField
200(
202 const scalar scale,
204) const
205{
206 for (const label pointi : pointLabels)
207 {
208 if (isInternalPoint_.test(pointi))
209 {
210 fld[pointi] *= scale;
211 }
212 }
213
214 // Single and multi-patch constraints
215 pointConstraints::New(pMesh()).constrain(fld, false);
216}
217
218
219// Scale on selected points (usually patch points)
220void Foam::motionSmootherAlgo::scaleField
221(
222 const labelList& meshPoints,
224 const scalar scale,
226) const
227{
228 for (const label pointi : meshPoints)
229 {
230 if (pointLabels.found(pointi))
231 {
232 fld[pointi] *= scale;
233 }
234 }
235}
236
237
238// Lower on internal points
239void Foam::motionSmootherAlgo::subtractField
240(
242 const scalar f,
244) const
245{
246 for (const label pointi : pointLabels)
247 {
248 if (isInternalPoint_.test(pointi))
249 {
250 fld[pointi] = max(0.0, fld[pointi]-f);
251 }
252 }
253
254 // Single and multi-patch constraints
255 pointConstraints::New(pMesh()).constrain(fld);
256}
257
258
259// Scale on selected points (usually patch points)
260void Foam::motionSmootherAlgo::subtractField
261(
262 const labelList& meshPoints,
264 const scalar f,
266) const
267{
268 for (const label pointi : meshPoints)
269 {
270 if (pointLabels.found(pointi))
271 {
272 fld[pointi] = max(0.0, fld[pointi]-f);
273 }
274 }
275}
276
277
278void Foam::motionSmootherAlgo::getAffectedFacesAndPoints
279(
280 const label nPointIter,
281 const faceSet& wrongFaces,
282
283 labelList& affectedFaces,
284 bitSet& isAffectedPoint
285) const
286{
287 isAffectedPoint.reset();
288 isAffectedPoint.resize(mesh_.nPoints());
289
290 faceSet nbrFaces(mesh_, "checkFaces", wrongFaces);
291
292 // Find possible points influenced by nPointIter iterations of
293 // scaling and smoothing by doing pointCellpoint walk.
294 // Also update faces-to-be-checked to extend one layer beyond the points
295 // that will get updated.
296
297 for (label i = 0; i < nPointIter; i++)
298 {
299 pointSet nbrPoints(mesh_, "grownPoints", getPoints(nbrFaces));
300
301 for (const label pointi : nbrPoints)
302 {
303 for (const label celli : mesh_.pointCells(pointi))
304 {
305 nbrFaces.set(mesh_.cells()[celli]); // set multiple
306 }
307 }
308 nbrFaces.sync(mesh_);
309
310 if (i == nPointIter - 2)
311 {
312 for (const label facei : nbrFaces)
313 {
314 isAffectedPoint.set(mesh_.faces()[facei]); // set multiple
315 }
316 }
317 }
318
319 affectedFaces = nbrFaces.toc();
320}
321
322
323// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
324
326(
327 polyMesh& mesh,
328 pointMesh& pMesh,
330 pointVectorField& displacement,
331 pointScalarField& scale,
332 pointField& oldPoints,
333 const labelList& adaptPatchIDs,
334 const dictionary& paramDict,
335 const bool dryRun
336)
337:
338 mesh_(mesh),
339 pMesh_(pMesh),
340 pp_(pp),
341 displacement_(displacement),
342 scale_(scale),
343 oldPoints_(oldPoints),
344 adaptPatchIDs_(adaptPatchIDs),
345 paramDict_(paramDict),
346 dryRun_(dryRun),
347 isInternalPoint_(mesh_.nPoints(), true)
348{
349 updateMesh();
350}
351
352
353// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
354
356{}
357
358
359// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
360
362{
363 return mesh_;
364}
365
366
368{
369 return pMesh_;
370}
371
372
374{
375 return pp_;
376}
377
378
380{
381 return adaptPatchIDs_;
382}
383
384
386{
387 return paramDict_;
388}
389
390
392{
393 oldPoints_ = mesh_.points();
394
395 scale_ = 1.0;
396
397 // No need to update twoDmotion corrector since only holds edge labels
398 // which will remain the same as before. So unless the mesh was distorted
399 // severely outside of motionSmootherAlgo there will be no need.
400}
401
402
404(
405 const labelList& patchIDs,
406 pointVectorField& displacement
407)
408{
409 pointVectorField::Boundary& displacementBf =
410 displacement.boundaryFieldRef();
411
412 // Adapt the fixedValue bc's (i.e. copy internal point data to
413 // boundaryField for all affected patches)
414 for (const label patchi : patchIDs)
415 {
416 displacementBf[patchi] ==
417 displacementBf[patchi].patchInternalField();
418 }
419
420 // Make consistent with non-adapted bc's by evaluating those now and
421 // resetting the displacement from the values.
422 // Note that we're just doing a correctBoundaryConditions with
423 // fixedValue bc's first.
424 labelHashSet adaptPatchSet(patchIDs);
425
426 const lduSchedule& patchSchedule =
427 displacement.mesh().globalData().patchSchedule();
428
429 for (const auto& schedEval : patchSchedule)
430 {
431 const label patchi = schedEval.patch;
432
433 if (!adaptPatchSet.found(patchi))
434 {
435 if (schedEval.init)
436 {
437 displacementBf[patchi]
438 .initEvaluate(Pstream::commsTypes::scheduled);
439 }
440 else
441 {
442 displacementBf[patchi]
444 }
445 }
446 }
447
448 // Multi-patch constraints
449 pointConstraints::New(displacement.mesh()).constrainCorners(displacement);
450
451 // Adapt the fixedValue bc's (i.e. copy internal point data to
452 // boundaryField for all affected patches) to take the changes caused
453 // by multi-corner constraints into account.
454 for (const label patchi : patchIDs)
455 {
456 displacementBf[patchi] == displacementBf[patchi].patchInternalField();
457 }
458}
459
460
462{
463 setDisplacementPatchFields(adaptPatchIDs_, displacement_);
464}
465
466
468(
469 const labelList& patchIDs,
470 const indirectPrimitivePatch& pp,
471 pointField& patchDisp,
472 pointVectorField& displacement
473)
474{
475 const polyMesh& mesh = displacement.mesh()();
476
477 // See comment in .H file about shared points.
478 // We want to disallow effect of loose coupled points - we only
479 // want to see effect of proper fixedValue boundary conditions
480
481 const labelList& cppMeshPoints =
483
484 const labelList& ppMeshPoints = pp.meshPoints();
485
486 // Knock out displacement on points which are not on pp but are coupled
487 // to them since we want 'proper' values from displacement to take
488 // precedence.
489 {
490 bitSet isPatchPoint(mesh.nPoints(), ppMeshPoints);
492 (
493 mesh,
494 isPatchPoint,
496 0u
497 );
498
499 for (const label pointi : cppMeshPoints)
500 {
501 if (isPatchPoint.test(pointi))
502 {
503 displacement[pointi] = Zero;
504 }
505 }
506 }
507
508
509 // Set internal point data from displacement on combined patch points.
510 forAll(ppMeshPoints, patchPointi)
511 {
512 displacement[ppMeshPoints[patchPointi]] = patchDisp[patchPointi];
513 }
514
515
516 // Combine any coupled points
518 (
519 mesh,
520 displacement,
521 maxMagEqOp(), // combine op
522 vector::zero // null value
523 );
524
525
526 // Adapt the fixedValue bc's (i.e. copy internal point data to
527 // boundaryField for all affected patches)
528 setDisplacementPatchFields(patchIDs, displacement);
529
530
531 if (debug)
532 {
533 OFstream str(mesh.db().path()/"changedPoints.obj");
534 label nVerts = 0;
535 forAll(ppMeshPoints, patchPointi)
536 {
537 const vector& newDisp = displacement[ppMeshPoints[patchPointi]];
538
539 if (mag(newDisp-patchDisp[patchPointi]) > SMALL)
540 {
541 const point& pt = mesh.points()[ppMeshPoints[patchPointi]];
542
543 meshTools::writeOBJ(str, pt);
544 nVerts++;
545 //Pout<< "Point:" << pt
546 // << " oldDisp:" << patchDisp[patchPointi]
547 // << " newDisp:" << newDisp << endl;
548 }
549 }
550 Pout<< "Written " << nVerts << " points that are changed to file "
551 << str.name() << endl;
552 }
553
554 // Now reset input displacement
555 forAll(ppMeshPoints, patchPointi)
556 {
557 patchDisp[patchPointi] = displacement[ppMeshPoints[patchPointi]];
558 }
559}
560
561
563{
564 setDisplacement(adaptPatchIDs_, pp_, patchDisp, displacement_);
565}
566
567
568// correctBoundaryConditions with fixedValue bc's first.
570(
571 pointVectorField& displacement
572) const
573{
574 labelHashSet adaptPatchSet(adaptPatchIDs_);
575
576 const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
577
578 auto& displacementBf = displacement.boundaryFieldRef();
579
580 // 1. evaluate on adaptPatches
581 for (const auto& schedEval : patchSchedule)
582 {
583 const label patchi = schedEval.patch;
584
585 if (adaptPatchSet.found(patchi))
586 {
587 if (schedEval.init)
588 {
589 displacementBf[patchi]
590 .initEvaluate(Pstream::commsTypes::blocking);
591 }
592 else
593 {
594 displacementBf[patchi]
596 }
597 }
598 }
599
600
601 // 2. evaluate on non-AdaptPatches
602 for (const auto& schedEval : patchSchedule)
603 {
604 const label patchi = schedEval.patch;
605
606 if (!adaptPatchSet.found(patchi))
607 {
608 if (schedEval.init)
609 {
610 displacementBf[patchi]
611 .initEvaluate(Pstream::commsTypes::blocking);
612 }
613 else
614 {
615 displacementBf[patchi]
617 }
618 }
619 }
620
621 // Multi-patch constraints
622 pointConstraints::New(displacement.mesh()).constrainCorners(displacement);
623
624 // Correct for problems introduced by corner constraints
626 (
627 mesh_,
628 displacement,
629 maxMagEqOp(), // combine op
630 vector::zero // null value
631 );
632}
633
634
635
637{
638 // Correct for 2-D motion
639 const twoDPointCorrector& twoDCorrector = twoDPointCorrector::New(mesh_);
640
641 if (twoDCorrector.required())
642 {
643 Info<< "Correcting 2-D mesh motion";
644
645 if (mesh_.globalData().parallel())
646 {
648 << "2D mesh-motion probably not correct in parallel" << endl;
649 }
650
651 // We do not want to move 3D planes so project all points onto those
652 const pointField& oldPoints = mesh_.points();
653 const edgeList& edges = mesh_.edges();
654
655 const labelList& neIndices = twoDCorrector.normalEdgeIndices();
656 const vector& pn = twoDCorrector.planeNormal();
657
658 for (const label edgei : neIndices)
659 forAll(neIndices, i)
660 {
661 const edge& e = edges[edgei];
662
663 point& pStart = newPoints[e.start()];
664
665 pStart += pn*(pn & (oldPoints[e.start()] - pStart));
666
667 point& pEnd = newPoints[e.end()];
668
669 pEnd += pn*(pn & (oldPoints[e.end()] - pEnd));
670 }
671
672 // Correct tangentially
673 twoDCorrector.correctPoints(newPoints);
674 Info<< " ...done" << endl;
675 }
676
677 if (debug)
678 {
679 Pout<< "motionSmootherAlgo::modifyMotionPoints :"
680 << " testing sync of newPoints."
681 << endl;
682 testSyncPositions(newPoints, 1e-6*mesh_.bounds().mag());
683 }
684}
685
686
688{
689 // Make sure to clear out tetPtIs since used in checks (sometimes, should
690 // really check)
691 mesh_.clearTetBasePtIs();
692 pp_.movePoints(mesh_.points());
693}
694
695
697(
698 const scalar errorReduction
699)
700{
701 scalar old = paramDict_.get<scalar>("errorReduction");
702
703 paramDict_.remove("errorReduction");
704 paramDict_.add("errorReduction", errorReduction);
705
706 return old;
707}
708
709
711(
712 labelList& checkFaces,
713 const bool smoothMesh,
714 const label nAllowableErrors
715)
716{
717 List<labelPair> emptyBaffles;
718 return scaleMesh
719 (
720 checkFaces,
721 emptyBaffles,
722 smoothMesh,
723 nAllowableErrors
724 );
725}
726
727
729(
730 labelList& checkFaces,
731 const List<labelPair>& baffles,
732 const bool smoothMesh,
733 const label nAllowableErrors
734)
735{
736 return scaleMesh
737 (
738 checkFaces,
739 baffles,
740 paramDict_,
741 paramDict_,
742 smoothMesh,
743 nAllowableErrors
744 );
745}
746
747
749{
750 // Set newPoints as old + scale*displacement
751
752 // Create overall displacement with same b.c.s as displacement_
753 wordList actualPatchTypes;
754 {
755 const pointBoundaryMesh& pbm = displacement_.mesh().boundary();
756 actualPatchTypes.setSize(pbm.size());
757 forAll(pbm, patchi)
758 {
759 actualPatchTypes[patchi] = pbm[patchi].type();
760 }
761 }
762
763 wordList actualPatchFieldTypes;
764 {
765 const pointVectorField::Boundary& pfld =
766 displacement_.boundaryField();
767 actualPatchFieldTypes.setSize(pfld.size());
768 forAll(pfld, patchi)
769 {
770 if (isA<fixedValuePointPatchField<vector>>(pfld[patchi]))
771 {
772 // Get rid of funny
773 actualPatchFieldTypes[patchi] =
775 }
776 else
777 {
778 actualPatchFieldTypes[patchi] = pfld[patchi].type();
779 }
780 }
781 }
782
783 pointVectorField totalDisplacement
784 (
786 (
787 "totalDisplacement",
788 mesh_.time().timeName(),
789 mesh_,
792 false
793 ),
794 scale_*displacement_,
795 actualPatchFieldTypes,
796 actualPatchTypes
797 );
798 correctBoundaryConditions(totalDisplacement);
799
800 if (debug)
801 {
802 Pout<< "scaleMesh : testing sync of totalDisplacement" << endl;
803 testSyncField
804 (
805 totalDisplacement,
806 maxMagEqOp(),
807 vector::zero, // null value
808 1e-6*mesh_.bounds().mag()
809 );
810 }
811
812 tmp<pointField> tnewPoints(oldPoints_ + totalDisplacement.primitiveField());
813
814 // Correct for 2-D motion
815 modifyMotionPoints(tnewPoints.ref());
816
817 return tnewPoints;
818}
819
820
822(
823 labelList& checkFaces,
824 const List<labelPair>& baffles,
825 const dictionary& paramDict,
826 const dictionary& meshQualityDict,
827 const bool smoothMesh,
828 const label nAllowableErrors
829)
830{
831 if (!smoothMesh && adaptPatchIDs_.empty())
832 {
834 << "You specified both no movement on the internal mesh points"
835 << " (smoothMesh = false)" << nl
836 << "and no movement on the patch (adaptPatchIDs is empty)" << nl
837 << "Hence nothing to adapt."
838 << exit(FatalError);
839 }
840
841 if (debug)
842 {
843 // Had a problem with patches moved non-synced. Check transformations.
844 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
845
846 Pout<< "Entering scaleMesh : coupled patches:" << endl;
847 forAll(patches, patchi)
848 {
849 if (patches[patchi].coupled())
850 {
851 const coupledPolyPatch& pp =
852 refCast<const coupledPolyPatch>(patches[patchi]);
853
854 Pout<< '\t' << patchi << '\t' << pp.name()
855 << " parallel:" << pp.parallel()
856 << " separated:" << pp.separated()
857 << " forwardT:" << pp.forwardT().size()
858 << endl;
859 }
860 }
861 }
862
863 const scalar errorReduction = get<scalar>
864 (
865 paramDict, "errorReduction", dryRun_, keyType::REGEX_RECURSIVE
866 );
867 const label nSmoothScale = get<label>
868 (
869 paramDict, "nSmoothScale", dryRun_, keyType::REGEX_RECURSIVE
870 );
871
872 // Note: displacement_ should already be synced already from setDisplacement
873 // but just to make sure.
875 (
876 mesh_,
877 displacement_,
878 maxMagEqOp(),
879 vector::zero // null value
880 );
881
882 Info<< "Moving mesh using displacement scaling :"
883 << " min:" << gMin(scale_.primitiveField())
884 << " max:" << gMax(scale_.primitiveField())
885 << endl;
886
887 // Get points using current displacement and scale. Optionally 2D corrected.
888 pointField newPoints(curPoints());
889
890 // Move. No need to do 2D correction since curPoints already done that.
891 mesh_.movePoints(newPoints);
892 movePoints();
893
894
895 // Check. Returns parallel number of incorrect faces.
896 faceSet wrongFaces(mesh_, "wrongFaces", mesh_.nFaces()/100+100);
897 checkMesh
898 (
899 false,
900 mesh_,
901 meshQualityDict,
902 checkFaces,
903 baffles,
904 wrongFaces,
905 dryRun_
906 );
907
908 if (returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
909 {
910 return true;
911 }
912 else
913 {
914 // Sync across coupled faces by extending the set.
915 wrongFaces.sync(mesh_);
916
917 // Special case:
918 // if errorReduction is set to zero, extend wrongFaces
919 // to face-Cell-faces to ensure quick return to previously valid mesh
920
921 if (mag(errorReduction) < SMALL)
922 {
923 labelHashSet newWrongFaces(wrongFaces);
924 for (const label facei : wrongFaces)
925 {
926 const label own = mesh_.faceOwner()[facei];
927 const cell& ownFaces = mesh_.cells()[own];
928
929 newWrongFaces.insert(ownFaces);
930
931 if (facei < mesh_.nInternalFaces())
932 {
933 const label nei = mesh_.faceNeighbour()[facei];
934 const cell& neiFaces = mesh_.cells()[nei];
935
936 newWrongFaces.insert(neiFaces);
937 }
938 }
939 wrongFaces.transfer(newWrongFaces);
940 wrongFaces.sync(mesh_);
941 }
942
943
944 // Find out points used by wrong faces and scale displacement.
945 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
946
947 pointSet usedPoints(mesh_, "usedPoints", getPoints(wrongFaces));
948 usedPoints.sync(mesh_);
949
950
951 // Grow a few layers to determine
952 // - points to be smoothed
953 // - faces to be checked in next iteration
954 bitSet isAffectedPoint(mesh_.nPoints());
955 getAffectedFacesAndPoints
956 (
957 nSmoothScale, // smoothing iterations
958 wrongFaces, // error faces
959 checkFaces,
960 isAffectedPoint
961 );
962
963 if (debug)
964 {
965 Pout<< "Faces in error:" << wrongFaces.size()
966 << " with points:" << usedPoints.size()
967 << endl;
968 }
969
970 if (adaptPatchIDs_.size())
971 {
972 // Scale conflicting patch points
973 scaleField(pp_.meshPoints(), usedPoints, errorReduction, scale_);
974 //subtractField(pp_.meshPoints(), usedPoints, 0.1, scale_);
975 }
976 if (smoothMesh)
977 {
978 // Scale conflicting internal points
979 scaleField(usedPoints, errorReduction, scale_);
980 //subtractField(usedPoints, 0.1, scale_);
981 }
982
983 scalarField eWeights(calcEdgeWeights(oldPoints_));
984
985 for (label i = 0; i < nSmoothScale; ++i)
986 {
987 if (adaptPatchIDs_.size())
988 {
989 // Smooth patch values
990 pointScalarField oldScale(scale_);
991 minSmooth
992 (
993 eWeights,
994 isAffectedPoint,
995 pp_.meshPoints(),
996 oldScale,
997 scale_
998 );
999 checkFld(scale_);
1000 }
1001 if (smoothMesh)
1002 {
1003 // Smooth internal values
1004 pointScalarField oldScale(scale_);
1005 minSmooth(eWeights, isAffectedPoint, oldScale, scale_);
1006 checkFld(scale_);
1007 }
1008 }
1009
1011 (
1012 mesh_,
1013 scale_,
1015 -GREAT // null value
1016 );
1017
1018
1019 if (debug)
1020 {
1021 Pout<< "scale_ after smoothing :"
1022 << " min:" << Foam::gMin(scale_)
1023 << " max:" << Foam::gMax(scale_)
1024 << endl;
1025 }
1026
1027 return false;
1028 }
1029}
1030
1031
1033{
1034 const pointBoundaryMesh& patches = pMesh_.boundary();
1035
1036 // Check whether displacement has fixed value b.c. on adaptPatchID
1037 for (const label patchi : adaptPatchIDs_)
1038 {
1039 if
1040 (
1041 !isA<fixedValuePointPatchVectorField>
1042 (
1043 displacement_.boundaryField()[patchi]
1044 )
1045 )
1046 {
1048 << "Patch " << patches[patchi].name()
1049 << " has wrong boundary condition "
1050 << displacement_.boundaryField()[patchi].type()
1051 << " on field " << displacement_.name() << nl
1052 << "Only type allowed is "
1053 << fixedValuePointPatchVectorField::typeName
1054 << exit(FatalError);
1055 }
1056 }
1057
1058
1059 // Determine internal points. Note that for twoD there are no internal
1060 // points so we use the points of adaptPatchIDs instead
1061
1062 isInternalPoint_.unset(pp_.meshPoints()); // unset multiple
1063
1064 // Calculate master edge addressing
1065 isMasterEdge_ = syncTools::getMasterEdges(mesh_);
1066}
1067
1068
1069// ************************************************************************* //
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
const Mesh & mesh() const
Return mesh.
void evaluate()
Evaluate boundary conditions.
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
void correctBoundaryConditions()
Correct boundary field.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
bool found(const Key &key) const
Return true if hashed entry is found in table.
Definition: HashTableI.H:100
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
const objectRegistry & db() const noexcept
Return the local objectRegistry.
Definition: IOobject.C:500
fileName path() const
The complete path.
Definition: IOobject.C:524
void setSize(const label n)
Alias for resize()
Definition: List.H:218
Output to file stream, using an OSstream.
Definition: OFstream.H:57
virtual const fileName & name() const
Read/write access to the name of the stream.
Definition: OSstream.H:107
A list of faces which address into the list of points.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
bool found(const T &val, label pos=0) const
True if the value if found in the list.
Definition: UListI.H:265
iterator end() noexcept
Return an iterator to end traversing the UList.
Definition: UListI.H:350
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
Definition: bitSetI.H:521
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
virtual bool separated() const
Are the planes separated.
virtual bool parallel() const
Are the cyclic planes parallel.
virtual const tensorField & forwardT() const
Return face transformation tensor.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
A list of face labels.
Definition: faceSet.H:54
virtual void sync(const polyMesh &mesh)
Sync faceSet across coupled patches.
Definition: faceSet.C:130
A FixedValue boundary condition for pointField.
const indirectPrimitivePatch & coupledPatch() const
Return patch of all coupled faces.
@ REGEX_RECURSIVE
Definition: keyType.H:87
Given a displacement moves the mesh by scaling the displacement back until there are no more mesh err...
void setDisplacementPatchFields()
Set patch fields on displacement to be consistent with.
tmp< pointField > curPoints() const
Get the current points (oldPoints+scale*displacement)
scalar setErrorReduction(const scalar)
Set the errorReduction (by how much to scale the displacement.
const labelList & adaptPatchIDs() const
Patch labels that are being adapted.
static void setDisplacement(const labelList &patchIDs, const indirectPrimitivePatch &pp, pointField &patchDisp, pointVectorField &displacement)
Set displacement field from displacement on patch points.
void correct()
Take over existing mesh position.
const polyMesh & mesh() const
Reference to mesh.
void movePoints()
Update for new mesh geometry.
const indirectPrimitivePatch & patch() const
Reference to patch.
const pointMesh & pMesh() const
Reference to pointMesh.
const dictionary & paramDict() const
bool scaleMesh(labelList &checkFaces, const bool smoothMesh=true, const label nAllow=0)
Move mesh with given scale. Return true if mesh ok or has.
void updateMesh()
Update for new mesh topology.
void modifyMotionPoints(pointField &newPoints) const
Apply optional point constraint (2d correction)
const word & name() const noexcept
The patch name.
Foam::pointBoundaryMesh.
const pointMesh & mesh() const noexcept
Return the mesh reference.
Mesh representing a set of points created from polyMesh.
Definition: pointMesh.H:55
const pointBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: pointMesh.H:114
A set of point labels.
Definition: pointSet.H:54
virtual void sync(const polyMesh &mesh)
Sync set across coupled patches. Adds coupled points to set.
Definition: pointSet.C:130
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1310
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
label nPoints() const noexcept
Number of mesh points.
static bitSet getMasterEdges(const polyMesh &mesh)
Definition: syncTools.C:97
static void syncPointPositions(const polyMesh &mesh, List< point > &positions, const CombineOp &cop, const point &nullValue)
Synchronize locations on all mesh points.
Definition: syncTools.H:201
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
A class for managing temporary objects.
Definition: tmp.H:65
T & ref() const
Definition: tmpI.H:227
Class applies a two-dimensional correction to mesh motion point field.
bool required() const
Is 2D correction required, i.e. is the mesh a wedge or slab.
void correctPoints(pointField &p) const
Correct motion points.
const labelList & normalEdgeIndices() const
Return indices of normal edges.
const vector & planeNormal() const
Return plane normal.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const polyBoundaryMesh & patches
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
label nPoints
label faceId(-1)
#define WarningInFunction
Report a warning using Foam::Warning.
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
Namespace for OpenFOAM.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
List< label > labelList
A List of labels.
Definition: List.H:66
const TargetType * isA(const Type &t)
Check if dynamic_cast to TargetType is possible.
Definition: typeInfo.H:197
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.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
errorManip< error > abort(error &err)
Definition: errorManip.H:144
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
error FatalError
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
Type gMin(const FieldField< Field, Type > &f)
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Type gMax(const FieldField< Field, Type > &f)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
labelList pointLabels(nPoints, -1)
cellMask correctBoundaryConditions()
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
A non-counting (dummy) refCount.
Definition: refCount.H:59