src/errorEstimation/errorDrivenRefinement/errorDrivenRefinement.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 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 \*---------------------------------------------------------------------------*/
  26
  27 #include "errorDrivenRefinement.H"
  28 #include "polyTopoChanger.H"
  29 #include "polyMesh.H"
  30 #include "primitiveMesh.H"
  31 #include "polyTopoChange.H"
  32 #include "addToRunTimeSelectionTable.H"
  33 #include "volFields.H"
  34 #include "surfaceFields.H"
  35 #include "evaluateError.H"
  36 #include "fvc.H"
  37 #include "mapPolyMesh.H"
  38 #include "topoCellLooper.H"
  39 #include "cellCuts.H"
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45     defineTypeNameAndDebug(errorDrivenRefinement, 0);
  46     addToRunTimeSelectionTable
  47     (
  48         polyMeshModifier,
  49         errorDrivenRefinement,
  50         dictionary
  51     );
  52 }
  53
  54
  55 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  56
  57
  58 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  59
  60 // Construct from dictionary
  61 Foam::errorDrivenRefinement::errorDrivenRefinement
  62 (
  63     const word& name,
  64     const dictionary& dict,
  65     const label index,
  66     const polyTopoChanger& mme
  67 )
  68 :
  69     polyMeshModifier(name, index, mme, false),
  70     refinementEngine_(topoChanger().mesh(), true),
  71     errorField_(dict.lookup("errorField"))
  72 {}
  73
  74
  75 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  76
  77 Foam::errorDrivenRefinement::~errorDrivenRefinement()
  78 {}
  79
  80
  81 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  82
  83
  84 bool Foam::errorDrivenRefinement::changeTopology() const
  85 {
  86     const Time& runTime = topoChanger().mesh().time();
  87
  88     if (runTime.foundObject<volVectorField>(errorField_))
  89     {
  90         if (debug)
  91         {
  92             Info<< "errorDrivenRefinement::changeTopology() : triggering topo"
  93                 << " change since found errorField "
  94                 << errorField_ << endl;
  95         }
  96
  97         return true;
  98     }
  99     else
 100     {
 101         if (debug)
 102         {
 103             Info<< "errorDrivenRefinement::changeTopology() : no topo"
 104                 << " change request from me since no errorField "
 105                 << errorField_ << endl;
 106         }
 107
 108         return false;
 109     }
 110 }
 111
 112
 113 void Foam::errorDrivenRefinement::setRefinement(polyTopoChange& ref) const
 114 {
 115     // Insert the coarsen/refinement instructions into the topological change
 116
 117     if (debug)
 118     {
 119         Info<< "errorDrivenRefinement::setRefinement(polyTopoChange& ref)"
 120             << endl;
 121     }
 122
 123     const polyMesh& mesh = topoChanger().mesh();
 124
 125     const Time& runTime = mesh.time();
 126
 127     if (debug)
 128     {
 129         Info<< "Looking up vector field with name " << errorField_ << endl;
 130     }
 131     const volVectorField& resError =
 132         runTime.lookupObject<volVectorField>(errorField_);
 133
 134     const volScalarField magResError = Foam::mag(resError);
 135
 136     scalar min = Foam::min(magResError).value();
 137     scalar max = Foam::max(magResError).value();
 138     scalar avg = Foam::average(magResError).value();
 139
 140     if (debug)
 141     {
 142         Info<< "Writing magResError" << endl;
 143         magResError.write();
 144
 145         Info<< "min:" << min << " max:" << max << " avg:" << avg << endl;
 146     }
 147
 148     // Get faces to remove and cells to refine based on error
 149     evaluateError refPattern
 150     (
 151         magResError,                        // Error on cells
 152         resError,                           // Error vector on cells
 153         fvc::interpolate(magResError),      // Error on faces
 154         refinementEngine_.getSplitFaces()   // Current live split faces
 155     );
 156
 157
 158     // Insert mesh refinement into polyTopoChange:
 159     // - remove split faces
 160     // - refine cells
 161
 162     // Give 'hint' of faces to remove to cell splitter.
 163     const labelList& candidates = refPattern.unsplitFaces();
 164     ////Hack:no unsplitting
 165     //labelList candidates;
 166
 167     labelList removedFaces(refinementEngine_.removeSplitFaces(candidates, ref));
 168
 169     // Now success will be for every candidates whether face has been removed.
 170     // Protect cells using face from refinement.
 171
 172     // List of protected cells
 173     boolList markedCell(mesh.nCells(), false);
 174
 175     forAll(removedFaces, i)
 176     {
 177         label faceI = removedFaces[i];
 178
 179         markedCell[mesh.faceOwner()[faceI]] = true;
 180
 181         if (mesh.isInternalFace(faceI))
 182         {
 183             markedCell[mesh.faceNeighbour()[faceI]] = true;
 184         }
 185     }
 186
 187     // Repack list of cells to refine.
 188     List<refineCell> refCells = refPattern.refCells();
 189
 190     label newRefCellI = 0;
 191
 192     forAll(refCells, refCellI)
 193     {
 194         label cellI = refCells[refCellI].cellNo();
 195
 196         if (!markedCell[cellI] && (newRefCellI != refCellI))
 197         {
 198             refCells[newRefCellI++] = refCells[refCellI];
 199         }
 200     }
 201
 202     if (debug)
 203     {
 204         Info<< "errorDrivenRefinement : shrinking refCells from "
 205             << refCells.size()
 206             << " to " << newRefCellI << endl;
 207     }
 208
 209     refCells.setSize(newRefCellI);
 210
 211     // Determine cut pattern using topological cell walker
 212     topoCellLooper cellWalker(mesh);
 213
 214     cellCuts cuts(mesh, cellWalker, refCells);
 215
 216     // Do actual splitting
 217     refinementEngine_.setRefinement(cuts, ref);
 218 }
 219
 220
 221 // Has the responsability of moving my newly introduced points onto the right
 222 // place. This is since the whole mesh might e.g. have been moved by another
 223 // meshmodifier. So using preMotionPoints is hack for if I am only meshModifier.
 224 // Good solution:
 225 // - remember new point label of introduced point and vertices
 226 // of edge it is created from (in setRefinement)
 227 // - in here reposition point at correct position between current vertex
 228 // position of edge endpoints.
 229 void Foam::errorDrivenRefinement::modifyMotionPoints
 230 (
 231     pointField& motionPoints
 232 ) const
 233 {
 234     if (debug)
 235     {
 236         Info<< "errorDrivenRefinement::modifyMotionPoints(*pointField&)" << endl;
 237     }
 238 }
 239
 240
 241 void Foam::errorDrivenRefinement::updateMesh(const mapPolyMesh& morphMap)
 242 {
 243     // Mesh has changed topologically. Update local topological data
 244     if (debug)
 245     {
 246         Info<< "errorDrivenRefinement::updateMesh"
 247             << "(const mapPolyMesh& morphMap)" << endl;
 248     }
 249     refinementEngine_.updateMesh(morphMap);
 250 }
 251
 252
 253 void Foam::errorDrivenRefinement::write(Ostream& os) const
 254 {
 255     os  << nl << type() << nl;
 256 }
 257
 258
 259 void Foam::errorDrivenRefinement::writeDict(Ostream& os) const
 260 {
 261     os  << nl << name() << nl << token::BEGIN_BLOCK << nl
 262         << "    type " << type()
 263         << token::END_STATEMENT << nl
 264         << token::END_BLOCK << endl;
 265 }
 266
 267
 268 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 269
 270
 271 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 272
 273
 274 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 275
 276
 277 // ************************************************************************* //