src/lagrangian/molecularDynamics/potential/pairPotential/pairPotentialList/pairPotentialList.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 "pairPotentialList.H"
  28 #include "OFstream.H"
  29 #include "Time.H"
  30
  31 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  32
  33 void Foam::pairPotentialList::readPairPotentialDict
  34 (
  35     const dictionary& pairPotentialDict,
  36     const polyMesh& mesh
  37 )
  38 {
  39     Info<< nl << "Building pair potentials." << endl;
  40
  41     rCutMax_ = 0.0;
  42
  43     for (label a = 0; a < nIds(); ++a)
  44     {
  45         word idA = idList_[a];
  46
  47         for (label b = a; b < nIds(); ++b)
  48         {
  49             word idB = idList_[b];
  50
  51             word pairPotentialName;
  52
  53             if (a == b)
  54             {
  55                 if (pairPotentialDict.found(idA + "-" + idB))
  56                 {
  57                     pairPotentialName = idA + "-" + idB;
  58                 }
  59                 else
  60                 {
  61                     FatalErrorIn("pairPotentialList::buildPotentials") << nl
  62                             << "Pair pairPotential specification subDict "
  63                             << idA << "-" << idB << " not found"
  64                             << nl << abort(FatalError);
  65                 }
  66             }
  67             else
  68             {
  69                 if (pairPotentialDict.found(idA + "-" + idB))
  70                 {
  71                     pairPotentialName = idA + "-" + idB;
  72                 }
  73
  74                 else if (pairPotentialDict.found(idB + "-" + idA))
  75                 {
  76                     pairPotentialName = idB + "-" + idA;
  77                 }
  78
  79                 else
  80                 {
  81                     FatalErrorIn("pairPotentialList::buildPotentials") << nl
  82                         << "Pair pairPotential specification subDict "
  83                         << idA << "-" << idB << " or "
  84                         << idB << "-" << idA << " not found"
  85                         << nl << abort(FatalError);
  86                 }
  87
  88                 if
  89                 (
  90                     pairPotentialDict.found(idA + "-" + idB)
  91                  && pairPotentialDict.found(idB + "-" + idA)
  92                 )
  93                 {
  94                     FatalErrorIn("pairPotentialList::buildPotentials") << nl
  95                         << "Pair pairPotential specification subDict "
  96                         << idA << "-" << idB << " and "
  97                         << idB << "-" << idA << " found multiple definition"
  98                         << nl << abort(FatalError);
  99                 }
 100             }
 101
 102             (*this).set
 103             (
 104                 pairPotentialIndex(a, b),
 105                 pairPotential::New
 106                 (
 107                     pairPotentialName,
 108                     pairPotentialDict.subDict(pairPotentialName)
 109                 )
 110             );
 111
 112             if ((*this)[pairPotentialIndex(a, b)].rCut() > rCutMax_)
 113             {
 114                 rCutMax_ = (*this)[pairPotentialIndex(a, b)].rCut();
 115             }
 116
 117             if ((*this)[pairPotentialIndex(a, b)].writeTables())
 118             {
 119                 OFstream ppTabFile(mesh.time().path()/pairPotentialName);
 120
 121                 if
 122                 (
 123                     !(*this)[pairPotentialIndex(a, b)].writeEnergyAndForceTables
 124                     (
 125                         ppTabFile
 126                     )
 127                 )
 128                 {
 129                     FatalErrorIn("pairPotentialList::readPairPotentialDict")
 130                         << "Failed writing to "
 131                         << ppTabFile.name() << nl
 132                         << abort(FatalError);
 133                 }
 134             }
 135         }
 136     }
 137
 138     rCutMaxSqr_ = rCutMax_*rCutMax_;
 139 }
 140
 141 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 142
 143 Foam::pairPotentialList::pairPotentialList()
 144 :
 145     PtrList<pairPotential>(),
 146     idList_()
 147 {}
 148
 149
 150 Foam::pairPotentialList::pairPotentialList
 151 (
 152     const dictionary& idListDict,
 153     const dictionary& pairPotentialDict,
 154     const polyMesh& mesh
 155 )
 156 :
 157     PtrList<pairPotential>(),
 158     idList_()
 159 {
 160     buildPotentials(idListDict, pairPotentialDict, mesh);
 161 }
 162
 163
 164 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 165
 166 Foam::pairPotentialList::~pairPotentialList()
 167 {}
 168
 169
 170 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 171
 172 void Foam::pairPotentialList::buildPotentials
 173 (
 174     const dictionary& idListDict,
 175     const dictionary& pairPotentialDict,
 176     const polyMesh& mesh
 177 )
 178 {
 179     idList_ = List<word>(idListDict.lookup("idList"));
 180
 181     setSize(((idList_.size()*(idList_.size() + 1))/2));
 182
 183     readPairPotentialDict(pairPotentialDict, mesh);
 184 }
 185
 186
 187 const Foam::pairPotential& Foam::pairPotentialList::pairPotentialFunction
 188 (
 189     const label a,
 190     const label b
 191 ) const
 192 {
 193     return (*this)[pairPotentialIndex(a, b)];
 194 }
 195
 196
 197 bool Foam::pairPotentialList::rCutSqr
 198 (
 199     const label a,
 200     const label b,
 201     const scalar rIJMagSqr
 202 ) const
 203 {
 204     if (rIJMagSqr < rCutSqr(a, b))
 205     {
 206         return true;
 207     }
 208     else
 209     {
 210         return false;
 211     }
 212 }
 213
 214
 215 Foam::scalar Foam::pairPotentialList::rMin
 216 (
 217     const label a,
 218     const label b
 219 ) const
 220 {
 221     return (*this)[pairPotentialIndex(a, b)].rMin();
 222 }
 223
 224
 225 Foam::scalar Foam::pairPotentialList::dr
 226 (
 227     const label a,
 228     const label b
 229 ) const
 230 {
 231     return (*this)[pairPotentialIndex(a, b)].dr();
 232 }
 233
 234
 235 Foam::scalar Foam::pairPotentialList::rCutSqr
 236 (
 237     const label a,
 238     const label b
 239 ) const
 240 {
 241     return (*this)[pairPotentialIndex(a, b)].rCutSqr();
 242 }
 243
 244
 245 Foam::scalar Foam::pairPotentialList::rCut
 246 (
 247     const label a,
 248     const label b
 249 ) const
 250 {
 251     return (*this)[pairPotentialIndex(a, b)].rCut();
 252 }
 253
 254
 255 Foam::scalar Foam::pairPotentialList::force
 256 (
 257     const label a,
 258     const label b,
 259     const scalar rIJMag
 260 ) const
 261 {
 262     scalar f = (*this)[pairPotentialIndex(a, b)].forceLookup(rIJMag);
 263
 264     return f;
 265 }
 266
 267
 268 Foam::scalar Foam::pairPotentialList::energy
 269 (
 270     const label a,
 271     const label b,
 272     const scalar rIJMag
 273 ) const
 274 {
 275     scalar e = (*this)[pairPotentialIndex(a, b)].energyLookup(rIJMag);
 276
 277     return e;
 278 }
 279
 280
 281 // ************************************************************************* //