applications/solvers/solidMechanics/utilities/smoothMesh/smoothMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     smoothMesh
  26
  27 Description
  28     Smoothing mesh using Laplacian smoothing
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "argList.H"
  33 #include "timeSelector.H"
  34 #include "objectRegistry.H"
  35 #include "foamTime.H"
  36 #include "fvMesh.H"
  37 #include "boolList.H"
  38 #include "emptyPolyPatch.H"
  39 #include "twoDPointCorrector.H"
  40
  41 using namespace Foam;
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 int main(int argc, char *argv[])
  46 {
  47     timeSelector::addOptions();
  48
  49     Foam::argList::validOptions.insert("smoothPatches", "");
  50
  51 #   include "setRootCase.H"
  52 #   include "createTime.H"
  53     instantList timeDirs = timeSelector::select0(runTime, args);
  54 #   include "createMesh.H"
  55
  56     // bool smoothPatches = args.optionFound("smoothPatches");
  57     // if(smoothPatches)
  58     //   Info << "Smoothing patches" << nl << endl;
  59
  60     // Smooth internal points
  61
  62     const vectorField& oldPoints = mesh.points();
  63
  64     vectorField newPoints = oldPoints;
  65
  66     if (mesh.nGeometricD() == 3)
  67     {
  68         Info << "3-D mesh" << endl;
  69
  70         const labelListList& pointEdges = mesh.pointEdges();
  71         const edgeList& edges = mesh.edges();
  72
  73         boolList fixedPoints(newPoints.size(), false);
  74
  75         forAll(mesh.boundaryMesh(), patchI)
  76         {
  77             const labelList& meshPoints =
  78                 mesh.boundaryMesh()[patchI].meshPoints();
  79
  80             forAll(meshPoints, pointI)
  81             {
  82                 fixedPoints[meshPoints[pointI]] = true;
  83             }
  84         }
  85
  86         scalarField residual(newPoints.size(), 0);
  87         label counter = 0;
  88         do
  89         {
  90             counter++;
  91
  92             forAll(newPoints, pointI)
  93             {
  94                 if (!fixedPoints[pointI])
  95                 {
  96                     vector curNewPoint = vector::zero;
  97
  98                     scalar sumW = 0;
  99
 100                     forAll(pointEdges[pointI], eI)
 101                     {
 102                         label curEdgeIndex = pointEdges[pointI][eI];
 103
 104                         const edge& curEdge = edges[curEdgeIndex];
 105
 106                         vector d =
 107                             newPoints[curEdge.otherVertex(pointI)]
 108                           - newPoints[pointI];
 109
 110                         scalar w = 1.0;
 111
 112                         curNewPoint += w*d;
 113
 114                         sumW += w;
 115                     }
 116
 117                     curNewPoint /= sumW;
 118
 119                     curNewPoint += newPoints[pointI];
 120
 121                     residual[pointI] = mag(curNewPoint - newPoints[pointI]);
 122
 123                     newPoints[pointI] = curNewPoint;
 124                 }
 125             }
 126
 127         // smoothed patch points independently
 128         // but we do not smooth points on the boundary of the patch
 129         // we reset these points
 130         // problem: we need to smooth feature edges separately
 131         // if(smoothPatches)
 132         //   {
 133         //      forAll(mesh.boundaryMesh(), patchI)
 134         //        {
 135         //          const labelList& meshPoints =
 136         //            mesh.boundaryMesh()[patchI].meshPoints();
 137
 138         //          forAll(meshPoints, pointI)
 139         //            {
 140         //          label pointID = meshPoints[pointI];
 141         //          vector curNewPoint = vector::zero;
 142         //          scalar sumW = 0;
 143
 144         //          forAll(pointEdges[pointID], eI)
 145         //            {
 146         //              label curEdgeIndex = pointEdges[pointID][eI];
 147         //              const edge& curEdge = edges[curEdgeIndex];
 148
 149         //              // use only boundary points
 150         //              label otherPointID = curEdge.otherVertex(pointID);
 151         //              if(fixedPoints[otherPointID])
 152         //                {
 153         //              vector d =
 154         //                newPoints[otherPointID]
 155         //                - newPoints[pointID];
 156
 157         //              scalar w = 1.0;
 158         //              curNewPoint += w*d;
 159         //              sumW += w;
 160         //                }
 161         //            }
 162
 163         //          curNewPoint /= sumW;
 164         //          curNewPoint += newPoints[pointID];
 165         //          residual[pointID] = mag(curNewPoint - newPoints[pointID]);
 166         //          newPoints[pointID] = curNewPoint;
 167         //            }
 168
 169         //          // reset boundary points
 170         //          const labelList& boundaryPoints =
 171             // mesh.boundaryMesh()[patchI].boundaryPoints();
 172         //          forAll(boundaryPoints, bpi)
 173         //            {
 174         //          label boundaryPointID = meshPoints[boundaryPoints[bpi]];
 175         //          newPoints[boundaryPointID] = oldPoints[boundaryPointID];
 176         //          residual[boundaryPointID] = 0.0;
 177         //            }
 178         //        }
 179         //   }
 180
 181             residual /= max(mag(newPoints - oldPoints) + SMALL);
 182         }
 183         while(max(residual) > 1e-3);
 184
 185         Info << "Internal points, max residual: " << max(residual)
 186             << ", num of iterations: " << counter << endl;
 187
 188         twoDPointCorrector twoDCorrector(mesh);
 189         twoDCorrector.correctPoints(newPoints);
 190     }
 191     else // 2-D
 192     {
 193         Info << "2-D mesh" << endl;
 194
 195         forAll(mesh.boundaryMesh(), patchI)
 196         {
 197             if
 198             (
 199                 mesh.boundaryMesh()[patchI].type()
 200              == emptyPolyPatch::typeName
 201             )
 202             {
 203                 const labelList& bPoints =
 204                     mesh.boundaryMesh()[patchI].boundaryPoints();
 205
 206                 const vectorField& oldPatchPoints =
 207                     mesh.boundaryMesh()[patchI].localPoints();
 208
 209                 vectorField newPatchPoints = oldPatchPoints;
 210
 211                 const labelListList& patchPointEdges =
 212                     mesh.boundaryMesh()[patchI].pointEdges();
 213
 214                 const edgeList& patchEdges =
 215                     mesh.boundaryMesh()[patchI].edges();
 216
 217                 boolList fixedPoints(newPatchPoints.size(), false);
 218
 219                 forAll(bPoints, pointI)
 220                 {
 221                     fixedPoints[bPoints[pointI]] = true;
 222                 }
 223
 224                 scalarField residual(newPatchPoints.size(), 0);
 225                 label counter = 0;
 226                 do
 227                 {
 228                     counter++;
 229
 230                     forAll(newPatchPoints, pointI)
 231                     {
 232                         if (!fixedPoints[pointI])
 233                         {
 234                             vector curNewPatchPoint = vector::zero;
 235
 236                             scalar sumW = 0;
 237
 238                             forAll(patchPointEdges[pointI], eI)
 239                             {
 240                                 label curEdgeIndex =
 241                                     patchPointEdges[pointI][eI];
 242
 243                                 const edge& curEdge = patchEdges[curEdgeIndex];
 244
 245                                 vector d =
 246                                     newPatchPoints[curEdge.otherVertex(pointI)]
 247                                   - newPatchPoints[pointI];
 248
 249                                 scalar w = 1.0;
 250
 251                                 curNewPatchPoint += w*d;
 252
 253                                 sumW += w;
 254                             }
 255
 256                             curNewPatchPoint /= sumW;
 257
 258                             curNewPatchPoint += newPatchPoints[pointI];
 259
 260                             residual[pointI] =
 261                                 mag(curNewPatchPoint - newPatchPoints[pointI]);
 262
 263                             newPatchPoints[pointI] = curNewPatchPoint;
 264                         }
 265                     }
 266
 267                     residual /=
 268                         max(mag(newPatchPoints - oldPatchPoints) + SMALL);
 269                 }
 270                 while(max(residual) > 1e-6);
 271
 272                 Info << "Empty points, max residual: " << max(residual)
 273                     << ", num of iterations: " << counter << endl;
 274
 275                 const labelList& meshPoints =
 276                     mesh.boundaryMesh()[patchI].meshPoints();
 277
 278                 forAll(meshPoints, pointI)
 279                 {
 280                     newPoints[meshPoints[pointI]] = newPatchPoints[pointI];
 281                 }
 282             }
 283         }
 284     }
 285
 286
 287     twoDPointCorrector twoDCorrector(mesh);
 288     twoDCorrector.correctPoints(newPoints);
 289
 290     mesh.movePoints(newPoints);
 291
 292     runTime++;
 293
 294     runTime.write();
 295
 296     Info<< "End\n" << endl;
 297
 298     return 0;
 299 }
 300
 301
 302 // ************************************************************************* //