FreeFOAM The Cross-Platform CFD Toolkit
motionSmoother.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) 1991-2010 OpenCFD Ltd.
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 "motionSmoother.H"
28 #include <meshTools/faceSet.H>
29 #include <meshTools/pointSet.H>
32 #include <OpenFOAM/syncTools.H>
33 #include <meshTools/meshTools.H>
34 #include <OpenFOAM/OFstream.H>
35 
36 namespace Foam
37 {
38 
39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
40 
41 defineTypeNameAndDebug(motionSmoother, 0);
42 
43 }
44 
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 // From pointPatchInterpolation
49 void Foam::motionSmoother::makePatchPatchAddressing()
50 {
51  if (debug)
52  {
53  Pout<< "motionSmoother::makePatchPatchAddressing() : "
54  << "constructing boundary addressing"
55  << endl;
56  }
57 
58  const polyBoundaryMesh& bm = mesh_.boundaryMesh();
59  const pointBoundaryMesh& pbm = pMesh_.boundary();
60 
61  // first count the total number of patch-patch points
62 
63  label nPatchPatchPoints = 0;
64 
65  forAll(bm, patchi)
66  {
67  if(!isA<emptyPolyPatch>(bm[patchi]))
68  {
69  nPatchPatchPoints += bm[patchi].boundaryPoints().size();
70  }
71  }
72 
73 
74  // Go through all patches and mark up the external edge points
75  Map<label> patchPatchPointSet(2*nPatchPatchPoints);
76 
77  labelList patchPatchPoints(nPatchPatchPoints);
78 
79  List<pointConstraint> patchPatchPointConstraints(nPatchPatchPoints);
80 
81  label pppi = 0;
82 
83  forAll(bm, patchi)
84  {
85  if(!isA<emptyPolyPatch>(bm[patchi]))
86  {
87  const labelList& bp = bm[patchi].boundaryPoints();
88  const labelList& meshPoints = bm[patchi].meshPoints();
89 
90  forAll(bp, pointi)
91  {
92  label ppp = meshPoints[bp[pointi]];
93 
94  Map<label>::iterator iter = patchPatchPointSet.find(ppp);
95 
96  if (iter == patchPatchPointSet.end())
97  {
98  patchPatchPointSet.insert(ppp, pppi);
99  patchPatchPoints[pppi] = ppp;
100  pbm[patchi].applyConstraint
101  (
102  bp[pointi],
103  patchPatchPointConstraints[pppi]
104  );
105  pppi++;
106  }
107  else
108  {
109  pbm[patchi].applyConstraint
110  (
111  bp[pointi],
112  patchPatchPointConstraints[iter()]
113  );
114  }
115  }
116  }
117  }
118 
119  nPatchPatchPoints = pppi;
120  patchPatchPoints.setSize(nPatchPatchPoints);
121  patchPatchPointConstraints.setSize(nPatchPatchPoints);
122 
123  patchPatchPointConstraintPoints_.setSize(nPatchPatchPoints);
124  patchPatchPointConstraintTensors_.setSize(nPatchPatchPoints);
125 
126  label nConstraints = 0;
127 
128  forAll(patchPatchPointConstraints, i)
129  {
130  if (patchPatchPointConstraints[i].first() != 0)
131  {
132  patchPatchPointConstraintPoints_[nConstraints] =
133  patchPatchPoints[i];
134 
135  patchPatchPointConstraintTensors_[nConstraints] =
136  patchPatchPointConstraints[i].constraintTransformation();
137 
138  nConstraints++;
139  }
140  }
141 
142  patchPatchPointConstraintPoints_.setSize(nConstraints);
143  patchPatchPointConstraintTensors_.setSize(nConstraints);
144 
145 
146  if (debug)
147  {
148  OFstream str(mesh_.db().path()/"constraintPoints.obj");
149 
150  Pout<< "Dumping " << patchPatchPointConstraintPoints_.size()
151  << " constraintPoints to " << str.name() << endl;
152  forAll(patchPatchPointConstraintPoints_, i)
153  {
154  label pointI = patchPatchPointConstraintPoints_[i];
155 
156  meshTools::writeOBJ(str, mesh_.points()[pointI]);
157  }
158 
159  Pout<< "motionSmoother::makePatchPatchAddressing() : "
160  << "finished constructing boundary addressing"
161  << endl;
162  }
163 }
164 
165 
166 void Foam::motionSmoother::checkFld(const pointScalarField& fld)
167 {
168  forAll(fld, pointI)
169  {
170  const scalar val = fld[pointI];
171 
172  if ((val > -GREAT) && (val < GREAT))
173  {}
174  else
175  {
176  FatalErrorIn("motionSmoother::checkFld")
177  << "Problem : point:" << pointI << " value:" << val
178  << abort(FatalError);
179  }
180  }
181 }
182 
183 
184 Foam::labelHashSet Foam::motionSmoother::getPoints
185 (
186  const labelHashSet& faceLabels
187 ) const
188 {
189  labelHashSet usedPoints(mesh_.nPoints()/100);
190 
191  forAllConstIter(labelHashSet, faceLabels, iter)
192  {
193  const face& f = mesh_.faces()[iter.key()];
194 
195  forAll(f, fp)
196  {
197  usedPoints.insert(f[fp]);
198  }
199  }
200 
201  return usedPoints;
202 }
203 
204 
205 // Smooth on selected points (usually patch points)
206 void Foam::motionSmoother::minSmooth
207 (
208  const PackedBoolList& isAffectedPoint,
209  const labelList& meshPoints,
210  const pointScalarField& fld,
211  pointScalarField& newFld
212 ) const
213 {
214  tmp<pointScalarField> tavgFld = avg
215  (
216  fld,
217  scalarField(mesh_.nEdges(), 1.0), // uniform weighting
218  false // fld is not position.
219  );
220  const pointScalarField& avgFld = tavgFld();
221 
222  forAll(meshPoints, i)
223  {
224  label pointI = meshPoints[i];
225  if (isAffectedPoint.get(pointI) == 1)
226  {
227  newFld[pointI] = min
228  (
229  fld[pointI],
230  0.5*fld[pointI] + 0.5*avgFld[pointI]
231  );
232  }
233  }
234 
235  newFld.correctBoundaryConditions();
236  applyCornerConstraints(newFld);
237 }
238 
239 
240 // Smooth on all internal points
241 void Foam::motionSmoother::minSmooth
242 (
243  const PackedBoolList& isAffectedPoint,
244  const pointScalarField& fld,
245  pointScalarField& newFld
246 ) const
247 {
248  tmp<pointScalarField> tavgFld = avg
249  (
250  fld,
251  scalarField(mesh_.nEdges(), 1.0), // uniform weighting
252  false // fld is not position.
253  );
254  const pointScalarField& avgFld = tavgFld();
255 
256  forAll(fld, pointI)
257  {
258  if (isAffectedPoint.get(pointI) == 1 && isInternalPoint(pointI))
259  {
260  newFld[pointI] = min
261  (
262  fld[pointI],
263  0.5*fld[pointI] + 0.5*avgFld[pointI]
264  );
265  }
266  }
267 
268  newFld.correctBoundaryConditions();
269  applyCornerConstraints(newFld);
270 }
271 
272 
273 // Scale on selected points
274 void Foam::motionSmoother::scaleField
275 (
276  const labelHashSet& pointLabels,
277  const scalar scale,
278  pointScalarField& fld
279 ) const
280 {
281  forAllConstIter(labelHashSet, pointLabels, iter)
282  {
283  if (isInternalPoint(iter.key()))
284  {
285  fld[iter.key()] *= scale;
286  }
287  }
288  fld.correctBoundaryConditions();
289  applyCornerConstraints(fld);
290 }
291 
292 
293 // Scale on selected points (usually patch points)
294 void Foam::motionSmoother::scaleField
295 (
296  const labelList& meshPoints,
297  const labelHashSet& pointLabels,
298  const scalar scale,
299  pointScalarField& fld
300 ) const
301 {
302  forAll(meshPoints, i)
303  {
304  label pointI = meshPoints[i];
305 
306  if (pointLabels.found(pointI))
307  {
308  fld[pointI] *= scale;
309  }
310  }
311 }
312 
313 
314 bool Foam::motionSmoother::isInternalPoint(const label pointI) const
315 {
316  return isInternalPoint_.get(pointI) == 1;
317 }
318 
319 
320 void Foam::motionSmoother::getAffectedFacesAndPoints
321 (
322  const label nPointIter,
323  const faceSet& wrongFaces,
324 
325  labelList& affectedFaces,
326  PackedBoolList& isAffectedPoint
327 ) const
328 {
329  isAffectedPoint.setSize(mesh_.nPoints());
330  isAffectedPoint = 0;
331 
332  faceSet nbrFaces(mesh_, "checkFaces", wrongFaces);
333 
334  // Find possible points influenced by nPointIter iterations of
335  // scaling and smoothing by doing pointCellpoint walk.
336  // Also update faces-to-be-checked to extend one layer beyond the points
337  // that will get updated.
338 
339  for (label i = 0; i < nPointIter; i++)
340  {
341  pointSet nbrPoints(mesh_, "grownPoints", getPoints(nbrFaces.toc()));
342 
343  forAllConstIter(pointSet, nbrPoints, iter)
344  {
345  const labelList& pCells = mesh_.pointCells(iter.key());
346 
347  forAll(pCells, pCellI)
348  {
349  const cell& cFaces = mesh_.cells()[pCells[pCellI]];
350 
351  forAll(cFaces, cFaceI)
352  {
353  nbrFaces.insert(cFaces[cFaceI]);
354  }
355  }
356  }
357  nbrFaces.sync(mesh_);
358 
359  if (i == nPointIter - 2)
360  {
361  forAllConstIter(faceSet, nbrFaces, iter)
362  {
363  const face& f = mesh_.faces()[iter.key()];
364  forAll(f, fp)
365  {
366  isAffectedPoint.set(f[fp], 1);
367  }
368  }
369  }
370  }
371 
372  affectedFaces = nbrFaces.toc();
373 }
374 
375 
376 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
377 
378 Foam::motionSmoother::motionSmoother
379 (
380  polyMesh& mesh,
381  pointMesh& pMesh,
383  const labelList& adaptPatchIDs,
384  const dictionary& paramDict
385 )
386 :
387  mesh_(mesh),
388  pMesh_(pMesh),
389  pp_(pp),
390  adaptPatchIDs_(adaptPatchIDs),
391  paramDict_(paramDict),
392  displacement_
393  (
394  IOobject
395  (
396  "displacement",
397  mesh_.time().timeName(),
398  mesh_,
401  ),
402  pMesh_
403  ),
404  scale_
405  (
406  IOobject
407  (
408  "scale",
409  mesh_.time().timeName(),
410  mesh_,
413  ),
414  pMesh_,
415  dimensionedScalar("scale", dimless, 1.0)
416  ),
417  oldPoints_(mesh_.points()),
418  isInternalPoint_(mesh_.nPoints(), 1),
419  twoDCorrector_(mesh_)
420 {
421  updateMesh();
422 }
423 
424 
425 Foam::motionSmoother::motionSmoother
426 (
427  polyMesh& mesh,
429  const labelList& adaptPatchIDs,
430  const pointVectorField& displacement,
431  const dictionary& paramDict
432 )
433 :
434  mesh_(mesh),
435  pMesh_(const_cast<pointMesh&>(displacement.mesh())),
436  pp_(pp),
437  adaptPatchIDs_(adaptPatchIDs),
438  paramDict_(paramDict),
439  displacement_
440  (
441  IOobject
442  (
443  "displacement",
444  mesh_.time().timeName(),
445  mesh_,
448  ),
449  displacement
450  ),
451  scale_
452  (
453  IOobject
454  (
455  "scale",
456  mesh_.time().timeName(),
457  mesh_,
460  ),
461  pMesh_,
462  dimensionedScalar("scale", dimless, 1.0)
463  ),
464  oldPoints_(mesh_.points()),
465  isInternalPoint_(mesh_.nPoints(), 1),
466  twoDCorrector_(mesh_)
467 {
468  updateMesh();
469 }
470 
471 
472 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
473 
475 {}
476 
477 
478 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
479 
481 {
482  return mesh_;
483 }
484 
485 
487 {
488  return pMesh_;
489 }
490 
491 
493 {
494  return pp_;
495 }
496 
497 
499 {
500  return adaptPatchIDs_;
501 }
502 
503 
505 {
506  return paramDict_;
507 }
508 
509 
511 {
512  return displacement_;
513 }
514 
515 
517 {
518  return displacement_;
519 }
520 
521 
523 {
524  return scale_;
525 }
526 
527 
529 {
530  return oldPoints_;
531 }
532 
533 
535 {
536  oldPoints_ = mesh_.points();
537 
538  scale_ = 1.0;
539 
540  // No need to update twoDmotion corrector since only holds edge labels
541  // which will remain the same as before. So unless the mesh was distorted
542  // severely outside of motionSmoother there will be no need.
543 }
544 
545 
547 {
548  // See comment in .H file about shared points.
549  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
550 
551  forAll(patches, patchI)
552  {
553  const polyPatch& pp = patches[patchI];
554 
555  if (pp.coupled())
556  {
557  const labelList& meshPoints = pp.meshPoints();
558 
559  forAll(meshPoints, i)
560  {
561  displacement_[meshPoints[i]] = vector::zero;
562  }
563  }
564  }
565 
566  const labelList& ppMeshPoints = pp_.meshPoints();
567 
568  // Set internal point data from displacement on combined patch points.
569  forAll(ppMeshPoints, patchPointI)
570  {
571  displacement_[ppMeshPoints[patchPointI]] = patchDisp[patchPointI];
572  }
573 
574  // Adapt the fixedValue bc's (i.e. copy internal point data to
575  // boundaryField for all affected patches)
576  forAll(adaptPatchIDs_, i)
577  {
578  label patchI = adaptPatchIDs_[i];
579 
580  displacement_.boundaryField()[patchI] ==
581  displacement_.boundaryField()[patchI].patchInternalField();
582  }
583 
584  // Make consistent with non-adapted bc's by evaluating those now and
585  // resetting the displacement from the values.
586  // Note that we're just doing a correctBoundaryConditions with
587  // fixedValue bc's first.
588  labelHashSet adaptPatchSet(adaptPatchIDs_);
589 
590  const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
591 
592  forAll(patchSchedule, patchEvalI)
593  {
594  label patchI = patchSchedule[patchEvalI].patch;
595 
596  if (!adaptPatchSet.found(patchI))
597  {
598  if (patchSchedule[patchEvalI].init)
599  {
600  displacement_.boundaryField()[patchI]
601  .initEvaluate(Pstream::scheduled);
602  }
603  else
604  {
605  displacement_.boundaryField()[patchI]
606  .evaluate(Pstream::scheduled);
607  }
608  }
609  }
610 
611  // Multi-patch constraints
612  applyCornerConstraints(displacement_);
613 
614  // Correct for problems introduced by corner constraints
616  (
617  mesh_,
618  displacement_,
619  maxMagEqOp(), // combine op
620  vector::zero, // null value
621  false // no separation
622  );
623 
624  // Adapt the fixedValue bc's (i.e. copy internal point data to
625  // boundaryField for all affected patches) to take the changes caused
626  // by multi-corner constraints into account.
627  forAll(adaptPatchIDs_, i)
628  {
629  label patchI = adaptPatchIDs_[i];
630 
631  displacement_.boundaryField()[patchI] ==
632  displacement_.boundaryField()[patchI].patchInternalField();
633  }
634 
635  if (debug)
636  {
637  OFstream str(mesh_.db().path()/"changedPoints.obj");
638  label nVerts = 0;
639  forAll(ppMeshPoints, patchPointI)
640  {
641  const vector& newDisp = displacement_[ppMeshPoints[patchPointI]];
642 
643  if (mag(newDisp-patchDisp[patchPointI]) > SMALL)
644  {
645  const point& pt = mesh_.points()[ppMeshPoints[patchPointI]];
646 
647  meshTools::writeOBJ(str, pt);
648  nVerts++;
649  //Pout<< "Point:" << pt
650  // << " oldDisp:" << patchDisp[patchPointI]
651  // << " newDisp:" << newDisp << endl;
652  }
653  }
654  Pout<< "Written " << nVerts << " points that are changed to file "
655  << str.name() << endl;
656  }
657 
658  // Now reset input displacement
659  forAll(ppMeshPoints, patchPointI)
660  {
661  patchDisp[patchPointI] = displacement_[ppMeshPoints[patchPointI]];
662  }
663 }
664 
665 
666 // correctBoundaryConditions with fixedValue bc's first.
668 (
669  pointVectorField& displacement
670 ) const
671 {
672  labelHashSet adaptPatchSet(adaptPatchIDs_);
673 
674  const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
675 
676  // 1. evaluate on adaptPatches
677  forAll(patchSchedule, patchEvalI)
678  {
679  label patchI = patchSchedule[patchEvalI].patch;
680 
681  if (adaptPatchSet.found(patchI))
682  {
683  if (patchSchedule[patchEvalI].init)
684  {
685  displacement.boundaryField()[patchI]
686  .initEvaluate(Pstream::blocking);
687  }
688  else
689  {
690  displacement.boundaryField()[patchI]
691  .evaluate(Pstream::blocking);
692  }
693  }
694  }
695 
696 
697  // 2. evaluate on non-AdaptPatches
698  forAll(patchSchedule, patchEvalI)
699  {
700  label patchI = patchSchedule[patchEvalI].patch;
701 
702  if (!adaptPatchSet.found(patchI))
703  {
704  if (patchSchedule[patchEvalI].init)
705  {
706  displacement.boundaryField()[patchI]
707  .initEvaluate(Pstream::blocking);
708  }
709  else
710  {
711  displacement.boundaryField()[patchI]
712  .evaluate(Pstream::blocking);
713  }
714  }
715  }
716 
717  // Multi-patch constraints
718  applyCornerConstraints(displacement);
719 
720  // Correct for problems introduced by corner constraints
722  (
723  mesh_,
724  displacement,
725  maxMagEqOp(), // combine op
726  vector::zero, // null value
727  false // no separation
728  );
729 }
730 
731 
733 (
734  pointField& newPoints
735 )
736 {
737  // Correct for 2-D motion
738  if (twoDCorrector_.required())
739  {
740  Info<< "Correct-ing 2-D mesh motion";
741 
742  if (mesh_.globalData().parallel())
743  {
744  WarningIn("motionSmoother::movePoints(pointField& newPoints)")
745  << "2D mesh-motion probably not correct in parallel" << endl;
746  }
747 
748  // We do not want to move 3D planes so project all points onto those
749  const pointField& oldPoints = mesh_.points();
750  const edgeList& edges = mesh_.edges();
751 
752  const labelList& neIndices = twoDCorrector().normalEdgeIndices();
753  const vector& pn = twoDCorrector().planeNormal();
754 
755  forAll(neIndices, i)
756  {
757  const edge& e = edges[neIndices[i]];
758 
759  point& pStart = newPoints[e.start()];
760 
761  pStart += pn*(pn & (oldPoints[e.start()] - pStart));
762 
763  point& pEnd = newPoints[e.end()];
764 
765  pEnd += pn*(pn & (oldPoints[e.end()] - pEnd));
766  }
767 
768  // Correct tangentially
769  twoDCorrector_.correctPoints(newPoints);
770  Info<< " ...done" << endl;
771  }
772 
773  if (debug)
774  {
775  Pout<< "motionSmoother::movePoints : testing sync of newPoints."
776  << endl;
777  testSyncField
778  (
779  newPoints,
780  minEqOp<point>(), // combine op
781  vector(GREAT,GREAT,GREAT), // null
782  true, // separation
783  1E-6*mesh_.bounds().mag()
784  );
785  }
786 
787  tmp<scalarField> tsweptVol = mesh_.movePoints(newPoints);
788 
789  pp_.movePoints(mesh_.points());
790 
791  return tsweptVol;
792 }
793 
794 
796 (
797  const scalar errorReduction
798 )
799 {
800  scalar oldErrorReduction = readScalar(paramDict_.lookup("errorReduction"));
801 
802  paramDict_.remove("errorReduction");
803  paramDict_.add("errorReduction", errorReduction);
804 
805  return oldErrorReduction;
806 }
807 
808 
810 (
811  labelList& checkFaces,
812  const bool smoothMesh,
813  const label nAllowableErrors
814 )
815 {
816  List<labelPair> emptyBaffles;
817  return scaleMesh
818  (
819  checkFaces,
820  emptyBaffles,
821  smoothMesh,
822  nAllowableErrors
823  );
824 }
825 
826 
828 (
829  labelList& checkFaces,
830  const List<labelPair>& baffles,
831  const bool smoothMesh,
832  const label nAllowableErrors
833 )
834 {
835  return scaleMesh
836  (
837  checkFaces,
838  baffles,
839  paramDict_,
840  paramDict_,
841  smoothMesh,
842  nAllowableErrors
843  );
844 }
845 
846 
848 (
849  labelList& checkFaces,
850  const List<labelPair>& baffles,
851  const dictionary& paramDict,
852  const dictionary& meshQualityDict,
853  const bool smoothMesh,
854  const label nAllowableErrors
855 )
856 {
857  if (!smoothMesh && adaptPatchIDs_.empty())
858  {
859  FatalErrorIn("motionSmoother::scaleMesh(const bool")
860  << "You specified both no movement on the internal mesh points"
861  << " (smoothMesh = false)" << nl
862  << "and no movement on the patch (adaptPatchIDs is empty)" << nl
863  << "Hence nothing to adapt."
864  << exit(FatalError);
865  }
866 
867  if (debug)
868  {
869  // Had a problem with patches moved non-synced. Check transformations.
870  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
871 
872  Pout<< "Entering scaleMesh : coupled patches:" << endl;
873  forAll(patches, patchI)
874  {
875  if (patches[patchI].coupled())
876  {
877  const coupledPolyPatch& pp =
878  refCast<const coupledPolyPatch>(patches[patchI]);
879 
880  Pout<< '\t' << patchI << '\t' << pp.name()
881  << " parallel:" << pp.parallel()
882  << " separated:" << pp.separated()
883  << " forwardT:" << pp.forwardT().size()
884  << endl;
885  }
886  }
887  }
888 
889  const scalar errorReduction =
890  readScalar(paramDict.lookup("errorReduction"));
891  const label nSmoothScale =
892  readLabel(paramDict.lookup("nSmoothScale"));
893 
894 
895  // Note: displacement_ should already be synced already from setDisplacement
896  // but just to make sure.
898  (
899  mesh_,
900  displacement_,
901  maxMagEqOp(),
902  vector::zero, // null value
903  false // no separation
904  );
905 
906  // Set newPoints as old + scale*displacement
907  pointField newPoints;
908  {
909  // Create overall displacement with same b.c.s as displacement_
910  pointVectorField totalDisplacement
911  (
912  IOobject
913  (
914  "totalDisplacement",
915  mesh_.time().timeName(),
916  mesh_,
919  false
920  ),
921  scale_*displacement_,
922  displacement_.boundaryField().types()
923  );
924  correctBoundaryConditions(totalDisplacement);
925 
926  if (debug)
927  {
928  Pout<< "scaleMesh : testing sync of totalDisplacement" << endl;
929  testSyncField
930  (
931  totalDisplacement,
932  maxMagEqOp(),
933  vector::zero, // null value
934  false, // separation
935  1E-6*mesh_.bounds().mag()
936  );
937  }
938 
939  newPoints = oldPoints_ + totalDisplacement.internalField();
940  }
941 
942  Info<< "Moving mesh using diplacement scaling :"
943  << " min:" << gMin(scale_.internalField())
944  << " max:" << gMax(scale_.internalField())
945  << endl;
946 
947 
948  // Move
949  movePoints(newPoints);
950 
951  // Check. Returns parallel number of incorrect faces.
952  faceSet wrongFaces(mesh_, "wrongFaces", mesh_.nFaces()/100+100);
953  checkMesh(false, mesh_, meshQualityDict, checkFaces, baffles, wrongFaces);
954 
955  if (returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
956  {
957  return true;
958  }
959  else
960  {
961  // Sync across coupled faces by extending the set.
962  wrongFaces.sync(mesh_);
963 
964  // Special case:
965  // if errorReduction is set to zero, extend wrongFaces
966  // to face-Cell-faces to ensure quick return to previously valid mesh
967 
968  if (mag(errorReduction) < SMALL)
969  {
970  labelHashSet newWrongFaces(wrongFaces);
971  forAllConstIter(labelHashSet, wrongFaces, iter)
972  {
973  label own = mesh_.faceOwner()[iter.key()];
974  const cell& ownFaces = mesh_.cells()[own];
975 
976  forAll(ownFaces, cfI)
977  {
978  newWrongFaces.insert(ownFaces[cfI]);
979  }
980 
981  if (iter.key() < mesh_.nInternalFaces())
982  {
983  label nei = mesh_.faceNeighbour()[iter.key()];
984  const cell& neiFaces = mesh_.cells()[nei];
985 
986  forAll(neiFaces, cfI)
987  {
988  newWrongFaces.insert(neiFaces[cfI]);
989  }
990  }
991  }
992  wrongFaces.transfer(newWrongFaces);
993  wrongFaces.sync(mesh_);
994  }
995 
996 
997  // Find out points used by wrong faces and scale displacement.
998  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
999 
1000  pointSet usedPoints(mesh_, "usedPoints", getPoints(wrongFaces));
1001  usedPoints.sync(mesh_);
1002 
1003 
1004 
1005  // Grow a few layers to determine
1006  // - points to be smoothed
1007  // - faces to be checked in next iteration
1008  PackedBoolList isAffectedPoint(mesh_.nPoints());
1009  getAffectedFacesAndPoints
1010  (
1011  nSmoothScale, // smoothing iterations
1012  wrongFaces, // error faces
1013  checkFaces,
1014  isAffectedPoint
1015  );
1016 
1017  if (debug)
1018  {
1019  Pout<< "Faces in error:" << wrongFaces.size()
1020  << " with points:" << usedPoints.size()
1021  << endl;
1022  }
1023 
1024  if (adaptPatchIDs_.size())
1025  {
1026  // Scale conflicting patch points
1027  scaleField(pp_.meshPoints(), usedPoints, errorReduction, scale_);
1028  }
1029  if (smoothMesh)
1030  {
1031  // Scale conflicting internal points
1032  scaleField(usedPoints, errorReduction, scale_);
1033  }
1034 
1035  for (label i = 0; i < nSmoothScale; i++)
1036  {
1037  if (adaptPatchIDs_.size())
1038  {
1039  // Smooth patch values
1040  pointScalarField oldScale(scale_);
1041  minSmooth
1042  (
1043  isAffectedPoint,
1044  pp_.meshPoints(),
1045  oldScale,
1046  scale_
1047  );
1048  checkFld(scale_);
1049  }
1050  if (smoothMesh)
1051  {
1052  // Smooth internal values
1053  pointScalarField oldScale(scale_);
1054  minSmooth(isAffectedPoint, oldScale, scale_);
1055  checkFld(scale_);
1056  }
1057  }
1058 
1060  (
1061  mesh_,
1062  scale_,
1063  maxEqOp<scalar>(),
1064  -GREAT, // null value
1065  false // no separation
1066  );
1067 
1068 
1069  if (debug)
1070  {
1071  Pout<< "scale_ after smoothing :"
1072  << " min:" << Foam::gMin(scale_)
1073  << " max:" << Foam::gMax(scale_)
1074  << endl;
1075  }
1076 
1077  return false;
1078  }
1079 }
1080 
1081 
1083 {
1084  const pointBoundaryMesh& patches = pMesh_.boundary();
1085 
1086  // Check whether displacement has fixed value b.c. on adaptPatchID
1087  forAll(adaptPatchIDs_, i)
1088  {
1089  label patchI = adaptPatchIDs_[i];
1090 
1091  if
1092  (
1093  !isA<fixedValuePointPatchVectorField>
1094  (
1095  displacement_.boundaryField()[patchI]
1096  )
1097  )
1098  {
1099  FatalErrorIn
1100  (
1101  "motionSmoother::motionSmoother"
1102  ) << "Patch " << patches[patchI].name()
1103  << " has wrong boundary condition "
1104  << displacement_.boundaryField()[patchI].type()
1105  << " on field " << displacement_.name() << nl
1106  << "Only type allowed is "
1107  << fixedValuePointPatchVectorField::typeName
1108  << exit(FatalError);
1109  }
1110  }
1111 
1112 
1113  // Determine internal points. Note that for twoD there are no internal
1114  // points so we use the points of adaptPatchIDs instead
1115 
1116  twoDCorrector_.updateMesh();
1117 
1118  const labelList& meshPoints = pp_.meshPoints();
1119 
1120  forAll(meshPoints, i)
1121  {
1122  isInternalPoint_.set(meshPoints[i], 0);
1123  }
1124 
1125  // Calculate master edge addressing
1126  isMasterEdge_ = syncTools::getMasterEdges(mesh_);
1127 
1128  makePatchPatchAddressing();
1129 }
1130 
1131 
1132 // Specialisation of applyCornerConstraints for scalars because
1133 // no constraint need be applied
1134 template<>
1135 void Foam::motionSmoother::applyCornerConstraints<Foam::scalar>
1137  GeometricField<scalar, pointPatchField, pointMesh>& pf
1138 ) const
1139 {}
1140 
1141 
1142 // ************************ vim: set sw=4 sts=4 et: ************************ //