src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Class
  26     mesquiteMotionSolver
  27
  28 Description
  29     Thin interface to the Mesquite Mesh Improvement library.
  30     Also provides surface-mesh smoothing capabilities.
  31
  32 SourceFiles
  33     mesquiteMotionSolver.C
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #ifndef mesquiteMotionSolver_H
  38 #define mesquiteMotionSolver_H
  39
  40 #include "motionSolver.H"
  41 #include "Time.H"
  42 #include "Map.H"
  43 #include "Switch.H"
  44 #include "edgeList.H"
  45 #include "labelPair.H"
  46 #include "pointIOField.H"
  47 #include "MeshObject.H"
  48 #include "HashSet.H"
  49
  50 // Have gcc ignore certain warnings while including mesquite headers
  51 #if defined(__GNUC__) && !defined(__INTEL_COMPILER)
  52 #   pragma GCC diagnostic ignored "-Wold-style-cast"
  53 #endif
  54
  55 #include "Mesquite_all_headers.hpp"
  56
  57 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  58
  59 namespace Foam
  60 {
  61
  62 // Class forward declarations
  63 class polyMesh;
  64
  65 /*---------------------------------------------------------------------------*\
  66                       Class mesquiteMotionSolver Declaration
  67 \*---------------------------------------------------------------------------*/
  68
  69 class mesquiteMotionSolver
  70 :
  71     public motionSolver,
  72     public MeshObject<polyMesh, mesquiteMotionSolver>
  73 {
  74
  75     // Private data
  76
  77         //- The mesh to be moved
  78         const polyMesh& Mesh_;
  79
  80         //- Mesh characteristics [2D/3D]
  81         Switch twoDMesh_;
  82
  83         //- Flag to determine if arrays have been initialized
  84         bool arraysInitialized_;
  85
  86         //- Number of points
  87         unsigned long nPoints_;
  88
  89         //- Number of elements
  90         unsigned long nCells_;
  91
  92         //- Number of auxiliary points
  93         labelList nAuxPoints_;
  94
  95         //- Number of auxiliary elements
  96         labelList nAuxCells_;
  97
  98         //- Switch to toggle surface smoothing
  99         Switch surfaceSmoothing_;
 100
 101         //- Switch to toggle volume correction
 102         Switch volumeCorrection_;
 103
 104         //- Specify tolerance for volume correction
 105         scalar volCorrTolerance_;
 106
 107         //- Specify max volume correction iterations
 108         label volCorrMaxIter_;
 109
 110         // Specify tolerance for the CG solver
 111         scalar tolerance_;
 112
 113         //- Specify multiple mesh-motion sweeps
 114         label nSweeps_;
 115
 116         //- Specify the surface-smoothing interval
 117         label surfInterval_;
 118
 119         //- Specify a relaxation factor for surface-smoothing
 120         scalar relax_;
 121
 122         //- Vertex coordinate array passed into Mesquite
 123         mutable List<double> vtxCoords_;
 124
 125         //- Connectivity array passed into mesquite
 126         mutable List<unsigned long> cellToNode_;
 127
 128         //- Flag array for vertices (fixed/free)
 129         mutable List<int> fixFlags_;
 130
 131         //- Array with element types
 132         mutable List<Mesquite::EntityTopology> mixedTypes_;
 133
 134         //- Reference points
 135         mutable pointIOField refPoints_;
 136
 137         //- Original points prior to surface-smoothing
 138         mutable pointField origPoints_;
 139
 140         //- Quality metric list
 141         wordList qMetrics_;
 142
 143         //- Pointers to base element metric class
 144         HashTable<autoPtr<Mesquite::QualityMetric> > qMetricTable_;
 145
 146         //- Pointer to base objective function
 147         autoPtr<Mesquite::ObjectiveFunction> objFunction_;
 148
 149         //- Pointer to base optimization algorithm (vertex mover)
 150         autoPtr<Mesquite::VertexMover> optAlgorithm_;
 151
 152         //- Termination criteria
 153         Mesquite::TerminationCriterion tcInner_;
 154         Mesquite::TerminationCriterion tcOuter_;
 155
 156         //- Data specific to Laplacian surface smoothing
 157         labelList pIDs_;
 158         labelList offsets_;
 159         List<vectorField> pNormals_;
 160         List<vectorField> gradEdge_;
 161         List<vectorField> localPts_;
 162         List<scalarField> edgeMarker_;
 163
 164         //- Data for the auxiliary entities
 165         labelList procIndices_;
 166         scalarField pointFractions_;
 167         labelListList procPointLabels_;
 168         List<List<labelPair> > globalProcPoints_;
 169
 170         List<Map<label> > unmatchedRecvPoints_;
 171         List<vectorField> sendSurfFields_, recvSurfFields_;
 172         List<Map<label> > sendSurfPointMap_, recvSurfPointMap_;
 173
 174         List<Map<label> > sendPointMap_, recvPointMap_;
 175         List<labelList> sendCellToNode_, recvCellToNode_;
 176         List<vectorField> sendPointBuffer_, recvPointBuffer_;
 177
 178         //- Explicit patch coupling for surface smoothing
 179         Map<label> patchCoupling_;
 180
 181         // Fields for the CG solver
 182         vectorField bV_;
 183         vectorField xV_;
 184         vectorField pV_;
 185         vectorField rV_;
 186         vectorField wV_;
 187         vectorField bdy_;
 188         scalarField pointMarker_;
 189
 190         scalar oldVolume_;
 191
 192     // Private Member Functions
 193
 194         // Sparse Matrix multiply
 195         void A(const vectorField& p, vectorField& w);
 196
 197         // Dot-product
 198         scalar dot(const vectorField& f1, const vectorField& f2);
 199
 200         // CG solver
 201         label CG
 202         (
 203             const vectorField& b,
 204             vectorField& p,
 205             vectorField& r,
 206             vectorField& w,
 207             vectorField& x
 208         );
 209
 210         // Compute the normalization factor for the matrix
 211         scalar normFactor
 212         (
 213             const vectorField& x,
 214             const vectorField& b,
 215             const vectorField& w,
 216             vectorField& tmpField
 217         );
 218
 219         // Component-wise sumMag
 220         scalar cmptSumMag(const vectorField& field);
 221
 222         // Transfer buffers for surface point fields
 223         void transferBuffers(vectorField &field);
 224
 225         // Apply boundary conditions
 226         void applyBCs(vectorField &field);
 227
 228         // Private member function to compute connectivity data
 229         void initArrays();
 230
 231         // Private member function to construct parallel connectivity data
 232         void initParallelConnectivity();
 233
 234         // Parallel send
 235         void parWrite(const label toID, const label& data);
 236
 237         // Parallel send (for lists)
 238         template <class Type>
 239         void parWrite(const label toID, const UList<Type>& data);
 240
 241         // Parallel receive
 242         void parRead(const label fromID, label& data);
 243
 244         // Parallel receive (for lists)
 245         template <class Type>
 246         void parRead(const label fromID, UList<Type>& data);
 247
 248         // Identify coupled patches
 249         void identifyCoupledPatches();
 250
 251         // Prepare for parallel surface smoothing
 252         void initParallelSurfaceSmoothing();
 253
 254         // Prepare mesquite connectivity for parallel runs
 255         void initMesquiteParallelArrays();
 256
 257         // Copy auxiliary points to/from buffers
 258         void copyAuxiliaryPoints(bool firstCopy);
 259
 260         // Private member functions to read dictionary options
 261         void readOptions();
 262
 263         // Apply fixed-value boundary conditions, if any.
 264         void applyFixedValuePatches();
 265
 266         // Private member function to perform Laplacian surface smoothing
 267         void smoothSurfaces();
 268
 269         // Compute quality of a tetrahedral cell
 270         scalar tetQuality
 271         (
 272             const label cIndex,
 273             const pointField& pField,
 274             bool returnVolume = false
 275         );
 276
 277         // Private member function to check for invalid
 278         // cells and correct if necessary.
 279         void correctInvalidCells();
 280
 281         // Correct global volume
 282         void correctGlobalVolume();
 283
 284         // Enforce cylindrical constraints for slip-patches
 285         void enforceCylindricalConstraints();
 286
 287         // Prepare point-normals with updated point positions
 288         void preparePointNormals();
 289
 290         // Utility method to check validity of cells connected to a point.
 291         bool checkValidity
 292         (
 293             const vector& x,
 294             const labelList& jList,
 295             scalar& beta
 296         );
 297
 298         // Update on topology change
 299         void update(const mapPolyMesh&);
 300
 301 public:
 302
 303     //- Runtime type information
 304     TypeName("mesquiteMotionSolver");
 305
 306
 307     // Constructors
 308
 309         //- Construct from polyMesh
 310         mesquiteMotionSolver(const polyMesh& mesh);
 311
 312         mesquiteMotionSolver(const polyMesh& mesh, Istream& msData);
 313
 314     // Destructor
 315
 316         virtual ~mesquiteMotionSolver();
 317
 318
 319     // Member Functions
 320
 321         //- Return reference to the mesh to be moved
 322         const polyMesh& mesh() const
 323         {
 324             return Mesh_;
 325         }
 326
 327         //- Return point location obtained from the current motion field
 328         virtual tmp<pointField> curPoints() const;
 329
 330         //- Solve for motion
 331         virtual void solve();
 332
 333         //- Move points
 334         virtual bool movePoints() const;
 335
 336         //- Update topology
 337         virtual void updateMesh(const mapPolyMesh&);
 338
 339         //- Update topology (using meshObjectBase)
 340         virtual bool updateMesh(const mapPolyMesh&) const;
 341 };
 342
 343
 344 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 345
 346 } // End namespace Foam
 347
 348 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 349
 350 #endif
 351
 352 // ************************************************************************* //