applications/solvers/solidMechanics/icoFsiElasticNonLinULSolidFoam/moveSolidMeshLeastSquares.H

   1 //--------------------------------------------------//
   2 //- move mesh
   3 //--------------------------------------------------//
   4 if(min(J.internalField()) > 0)
   5 {
   6     Info << "Moving mesh using least squares interpolation" << endl;
   7
   8     leastSquaresVolPointInterpolation pointInterpolation(stressMesh);
   9
  10     // Create point mesh
  11     pointMesh pMesh(stressMesh);
  12
  13     wordList types
  14     (
  15         pMesh.boundary().size(),
  16         calculatedFvPatchVectorField::typeName
  17     );
  18
  19     pointVectorField pointDU
  20     (
  21         IOobject
  22         (
  23             "pointDU",
  24             runTime.timeName(),
  25             stressMesh
  26         ),
  27         pMesh,
  28         dimensionedVector("zero", dimLength, vector::zero),
  29         types
  30     );
  31
  32     pointInterpolation.interpolate(DU, pointDU);
  33
  34     const vectorField& pointDUI =
  35         pointDU.internalField();
  36
  37     //- Move mesh
  38     vectorField newPoints = stressMesh.allPoints();
  39
  40     forAll (pointDUI, pointI)
  41     {
  42         newPoints[pointI] += pointDUI[pointI];
  43     }
  44
  45     // Correct symmetryPlane points
  46
  47     forAll(stressMesh.boundaryMesh(), patchI)
  48     {
  49         if (isA<symmetryPolyPatch>(stressMesh.boundaryMesh()[patchI]))
  50         {
  51             const labelList& meshPoints =
  52                 stressMesh.boundaryMesh()[patchI].meshPoints();
  53
  54             vector avgN =
  55                 gAverage(stressMesh.boundaryMesh()[patchI].pointNormals());
  56
  57             vector i(1, 0, 0);
  58             vector j(0, 1, 0);
  59             vector k(0, 0, 1);
  60
  61             if (mag(avgN&i) > 0.95)
  62             {
  63                 forAll(meshPoints, pI)
  64                 {
  65                     newPoints[meshPoints[pI]].x() = 0;
  66                 }
  67             }
  68             else if (mag(avgN&j) > 0.95)
  69             {
  70                 forAll(meshPoints, pI)
  71                 {
  72                     newPoints[meshPoints[pI]].y() = 0;
  73                 }
  74             }
  75             else if (mag(avgN&k) > 0.95)
  76             {
  77                 forAll(meshPoints, pI)
  78                 {
  79                     newPoints[meshPoints[pI]].z() = 0;
  80                 }
  81             }
  82         }
  83     }
  84
  85 #   include "calcUnusedNewPoints.H"
  86
  87     twoDPointCorrector twoDCorrector(stressMesh);
  88     twoDCorrector.correctPoints(newPoints);
  89     stressMesh.movePoints(newPoints);
  90     stressMesh.V00();
  91     stressMesh.moving(false);
  92 }
  93 else
  94 {
  95     FatalErrorIn(args.executable())
  96         << "Negative Jacobian"
  97         << exit(FatalError);
  98 }