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