applications/solvers/solidMechanics/elasticNonLinULSolidFoam/moveMeshLeastSquares.H

   1 //--------------------------------------------------//
   2 //- move mesh
   3 //--------------------------------------------------//
   4 {
   5     //Info << "Moving mesh using least squares interpolation" << endl;
   6
   7     leastSquaresVolPointInterpolation pointInterpolation(mesh);
   8
   9     // Create point mesh
  10     pointMesh pMesh(mesh);
  11     //pointMesh pMesh(mesh);
  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             mesh
  26         ),
  27         pMesh,
  28         dimensionedVector("zero", dimLength, vector::zero),
  29         types
  30     );
  31
  32     pointInterpolation.interpolate(DU, pointDU);
  33
  34     //pointDU.write();
  35
  36     const vectorField& pointDUI = pointDU.internalField();
  37
  38     //- Move mesh
  39     //vectorField newPoints = mesh.allPoints();
  40     pointVectorField newPoints
  41     (
  42         IOobject
  43         (
  44             "newPoints",
  45             runTime.timeName(),
  46             mesh
  47         ),
  48         pMesh,
  49         dimensionedVector("zero", dimLength, vector::zero)
  50         //mesh.allPoints()
  51     );
  52     newPoints.internalField() = mesh.allPoints();
  53
  54     // note: allPoints will have more points than pointDU
  55     // if there are globalFaceZones
  56     forAll (pointDUI, pointI)
  57     {
  58         newPoints[pointI] += pointDUI[pointI];
  59     }
  60
  61     // Correct symmetryPlane points
  62
  63     forAll(mesh.boundaryMesh(), patchI)
  64     {
  65         if (isA<symmetryPolyPatch>(mesh.boundaryMesh()[patchI]))
  66         {
  67             const labelList& meshPoints =
  68                 mesh.boundaryMesh()[patchI].meshPoints();
  69
  70             vector avgN =
  71                 gAverage(mesh.boundaryMesh()[patchI].pointNormals());
  72
  73             vector i(1, 0, 0);
  74             vector j(0, 1, 0);
  75             vector k(0, 0, 1);
  76
  77             if (mag(avgN&i) > 0.95)
  78             {
  79                 forAll(meshPoints, pI)
  80                 {
  81                     newPoints[meshPoints[pI]].x() = 0;
  82                 }
  83             }
  84             else if (mag(avgN&j) > 0.95)
  85             {
  86                 forAll(meshPoints, pI)
  87                 {
  88                     newPoints[meshPoints[pI]].y() = 0;
  89                 }
  90             }
  91             else if (mag(avgN&k) > 0.95)
  92             {
  93                 forAll(meshPoints, pI)
  94                 {
  95                     newPoints[meshPoints[pI]].z() = 0;
  96                 }
  97             }
  98         }
  99     }
 100
 101
 102 #   include "calcUnusedNewPoints.H"
 103
 104 //     // now we make sure processor patches are exactly the same
 105 //     newPoints.correctBoundaryConditions();
 106
 107     twoDPointCorrector twoDCorrector(mesh);
 108     twoDCorrector.correctPoints(newPoints);
 109     mesh.movePoints(newPoints);
 110     mesh.V00();
 111     mesh.moving(false);
 112
 113     // Update n
 114     n = mesh.Sf()/mesh.magSf();
 115 }