src/errorEstimation/evaluateError/evaluateError.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 "evaluateError.H"
  28 #include "volFields.H"
  29 #include "surfaceFields.H"
  30 #include "refineCell.H"
  31
  32
  33 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  34
  35 // Construct null
  36 Foam::evaluateError::evaluateError()
  37 :
  38     unsplitFaces_(),
  39     refCells_()
  40 {}
  41
  42
  43 // Construct from components
  44 Foam::evaluateError::evaluateError
  45 (
  46     const volScalarField& cellError,
  47     const volVectorField& gradTheta,
  48     const surfaceScalarField& faceError,
  49     const labelList& candidateFaces
  50 )
  51 :
  52     unsplitFaces_(candidateFaces.size()),
  53     refCells_()
  54 {
  55     const polyMesh& mesh = cellError.mesh();
  56
  57     // picks up the error field and the gradient of the variable
  58     // and appends lists of cells to refine/unrefine based on the width of
  59     // standard deviation of the error distribution
  60
  61     // calculate the average error
  62     scalar avgError = cellError.average().value();
  63
  64     scalar squareError = sqr(cellError)().average().value();
  65     scalar deviation = sqrt(squareError - sqr(avgError));
  66
  67     Info<< "avgError:" << avgError
  68         << "  squareError:" << squareError
  69         << "  deviation:" << deviation
  70         << endl;
  71
  72     scalar ref = avgError + deviation;
  73     scalar unref = avgError - deviation;
  74
  75     Info<< "evaluateError : refinement criterion : " << ref << endl
  76         << "                unrefinement criterion : " << unref << endl;
  77
  78     // Coarsen mesh first.
  79     // Find out set of candidateFaces where error is above crit.
  80
  81     // Construct to filter unrefinement pattern
  82 //    removeFaces faceRemover(mesh);
  83
  84     // Keep track of unrefinement pattern.
  85     boolList markedFace(mesh.nFaces(), false);
  86
  87     label unsplitFaceI = 0;
  88
  89     // Subset candidate faces and update refinement pattern interference pattern
  90     forAll(candidateFaces, candidateFaceI)
  91     {
  92         label faceI = candidateFaces[candidateFaceI];
  93
  94         if (markedFace[faceI])
  95         {
  96             Info<< "evaluateError : protected candidate face:" << faceI
  97                 << endl;
  98         }
  99         else
 100         {
 101 //            if (faceError[faceI] < unref)
 102             if (unsplitFaceI < (candidateFaces.size()/2 + 1))
 103             {
 104                 unsplitFaces_[unsplitFaceI++] = faceI;
 105
 106 //                faceRemover.markAffectedFaces(faceI, markedFace);
 107             }
 108         }
 109     }
 110
 111     unsplitFaces_.setSize(unsplitFaceI);
 112
 113     // Now we have:
 114     // -unsplitFaces_: all the faces that will be removed
 115     // -markedFace   : all the faces affected by this removal.
 116     // From markedFace protect the cells using them.
 117
 118     boolList markedCells(mesh.nCells(), false);
 119
 120 //    forAll(markedFace, faceI)
 121 //    {
 122 //        if (markedFace[faceI])
 123 //        {
 124 //            markedCells[mesh.faceOwner()[faceI]] = true;
 125 //
 126 //            if (mesh.isInternalFace(faceI))
 127 //            {
 128 //                markedCells[mesh.faceNeighbour()[faceI]] = true;
 129 //            }
 130 //        }
 131 //    }
 132
 133     // Select the cells that need to be split.
 134     // Two pass: count first, select later.
 135
 136     label refCellI = 0;
 137
 138     forAll(cellError, cellI)
 139     {
 140         if ((cellError[cellI] > ref) && !markedCells[cellI])
 141         {
 142             refCellI++;
 143         }
 144     }
 145
 146     refCells_.setSize(refCellI);
 147
 148     refCellI = 0;
 149
 150     forAll(cellError, cellI)
 151     {
 152         if ((cellError[cellI] > ref) && !markedCells[cellI])
 153         {
 154             refCells_[refCellI++] = refineCell(cellI, gradTheta[cellI]);
 155         }
 156     }
 157
 158     Info<< "evaluateError : selected " << unsplitFaces_.size()
 159         << " faces out of " << candidateFaces.size() << " for removal" << endl;
 160     Info<< "evaluateError : selected " << refCells_.size()
 161         << " cells out of " << cellError.size() << " for refinement" << endl;
 162 }
 163
 164
 165 // ************************************************************************* //