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-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 "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 Foam::pairPotentialList::pairPotentialList
 150 (
 151     const dictionary& idListDict,
 152     const dictionary& pairPotentialDict,
 153     const polyMesh& mesh
 154 )
 155 :
 156     PtrList<pairPotential>(),
 157     idList_()
 158 {
 159     buildPotentials(idListDict, pairPotentialDict, mesh);
 160 }
 161
 162
 163 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 164
 165 Foam::pairPotentialList::~pairPotentialList()
 166 {}
 167
 168
 169 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 170
 171 void Foam::pairPotentialList::buildPotentials
 172 (
 173     const dictionary& idListDict,
 174     const dictionary& pairPotentialDict,
 175     const polyMesh& mesh
 176 )
 177 {
 178     idList_ = List<word>(idListDict.lookup("idList"));
 179
 180     setSize(((idList_.size()*(idList_.size() + 1))/2));
 181
 182     readPairPotentialDict(pairPotentialDict, mesh);
 183 }
 184
 185
 186 const Foam::pairPotential& Foam::pairPotentialList::pairPotentialFunction
 187 (
 188     const label a,
 189     const label b
 190 ) const
 191 {
 192     return (*this)[pairPotentialIndex (a, b)];
 193 }
 194
 195
 196 bool Foam::pairPotentialList::rCutSqr
 197 (
 198     const label a,
 199     const label b,
 200     const scalar rIJMagSqr
 201 ) const
 202 {
 203     if (rIJMagSqr <= rCutSqr (a, b))
 204     {
 205         return true;
 206     }
 207     else
 208     {
 209         return false;
 210     }
 211 }
 212
 213
 214 Foam::scalar Foam::pairPotentialList::rMin
 215 (
 216     const label a,
 217     const label b
 218 ) const
 219 {
 220     return (*this)[pairPotentialIndex (a, b)].rMin();
 221 }
 222
 223
 224 Foam::scalar Foam::pairPotentialList::dr
 225 (
 226     const label a,
 227     const label b
 228 ) const
 229 {
 230     return (*this)[pairPotentialIndex (a, b)].dr();
 231 }
 232
 233
 234 Foam::scalar Foam::pairPotentialList::rCutSqr
 235 (
 236     const label a,
 237     const label b
 238 ) const
 239 {
 240     return (*this)[pairPotentialIndex (a, b)].rCutSqr();
 241 }
 242
 243
 244 Foam::scalar Foam::pairPotentialList::rCut
 245 (
 246     const label a,
 247     const label b
 248 ) const
 249 {
 250     return (*this)[pairPotentialIndex (a, b)].rCut();
 251 }
 252
 253
 254 Foam::scalar Foam::pairPotentialList::force
 255 (
 256     const label a,
 257     const label b,
 258     const scalar rIJMag
 259 ) const
 260 {
 261     scalar f = (*this)[pairPotentialIndex (a, b)].forceLookup(rIJMag);
 262
 263     return f;
 264 }
 265
 266
 267 Foam::scalar Foam::pairPotentialList::energy
 268 (
 269     const label a,
 270     const label b,
 271     const scalar rIJMag
 272 ) const
 273 {
 274     scalar e = (*this)[pairPotentialIndex (a, b)].energyLookup(rIJMag);
 275
 276     return e;
 277 }
 278
 279 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 280
 281 // ************************************************************************* //