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