FreeFOAM The Cross-Platform CFD Toolkit
pairPotentialList.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2008-2010 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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your 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, see <http://www.gnu.org/licenses/>.
23 
24 \*----------------------------------------------------------------------------*/
25 
26 #include "pairPotentialList.H"
27 #include <OpenFOAM/OFstream.H>
28 #include <OpenFOAM/Time.H>
29 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
32 void Foam::pairPotentialList::readPairPotentialDict
33 (
34  const List<word>& idList,
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),
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  if (!pairPotentialDict.found("electrostatic"))
139  {
140  FatalErrorIn("pairPotentialList::buildPotentials") << nl
141  << "Pair pairPotential specification subDict electrostatic"
142  << nl << abort(FatalError);
143  }
144 
145  electrostaticPotential_ = pairPotential::New
146  (
147  "electrostatic",
148  pairPotentialDict.subDict("electrostatic")
149  );
150 
151  if (electrostaticPotential_->rCut() > rCutMax_)
152  {
153  rCutMax_ = electrostaticPotential_->rCut();
154  }
155 
156  if (electrostaticPotential_->writeTables())
157  {
158  OFstream ppTabFile(mesh.time().path()/"electrostatic");
159 
160  if(!electrostaticPotential_->writeEnergyAndForceTables(ppTabFile))
161  {
162  FatalErrorIn("pairPotentialList::readPairPotentialDict")
163  << "Failed writing to "
164  << ppTabFile.name() << nl
165  << abort(FatalError);
166  }
167  }
168 
169  rCutMaxSqr_ = rCutMax_*rCutMax_;
170 }
171 
172 
173 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
174 
176 :
178 {}
179 
180 
182 (
183  const List<word>& idList,
184  const dictionary& pairPotentialDict,
185  const polyMesh& mesh
186 )
187 :
189 {
190  buildPotentials(idList, pairPotentialDict, mesh);
191 }
192 
193 
194 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
195 
197 {}
198 
199 
200 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
201 
203 (
204  const List<word>& idList,
205  const dictionary& pairPotentialDict,
206  const polyMesh& mesh
207 )
208 {
209  setSize(((idList.size()*(idList.size() + 1))/2));
210 
211  nIds_ = idList.size();
212 
213  readPairPotentialDict(idList, pairPotentialDict, mesh);
214 }
215 
216 
218 (
219  const label a,
220  const label b
221 ) const
222 {
223  return (*this)[pairPotentialIndex(a, b)];
224 }
225 
226 
227 bool Foam::pairPotentialList::rCutMaxSqr(const scalar rIJMagSqr) const
228 {
229  if (rIJMagSqr < rCutMaxSqr_)
230  {
231  return true;
232  }
233  else
234  {
235  return false;
236  }
237 }
238 
239 
241 (
242  const label a,
243  const label b,
244  const scalar rIJMagSqr
245 ) const
246 {
247  if (rIJMagSqr < rCutSqr(a, b))
248  {
249  return true;
250  }
251  else
252  {
253  return false;
254  }
255 }
256 
257 
259 (
260  const label a,
261  const label b
262 ) const
263 {
264  return (*this)[pairPotentialIndex(a, b)].rMin();
265 }
266 
267 
268 Foam::scalar Foam::pairPotentialList::dr
269 (
270  const label a,
271  const label b
272 ) const
273 {
274  return (*this)[pairPotentialIndex(a, b)].dr();
275 }
276 
277 
279 (
280  const label a,
281  const label b
282 ) const
283 {
284  return (*this)[pairPotentialIndex(a, b)].rCutSqr();
285 }
286 
287 
289 (
290  const label a,
291  const label b
292 ) const
293 {
294  return (*this)[pairPotentialIndex(a, b)].rCut();
295 }
296 
297 
299 (
300  const label a,
301  const label b,
302  const scalar rIJMag
303 ) const
304 {
305  scalar f = (*this)[pairPotentialIndex(a, b)].force(rIJMag);
306 
307  return f;
308 }
309 
310 
312 (
313  const label a,
314  const label b,
315  const scalar rIJMag
316 ) const
317 {
318  scalar e = (*this)[pairPotentialIndex(a, b)].energy(rIJMag);
319 
320  return e;
321 }
322 
323 
324 // ************************ vim: set sw=4 sts=4 et: ************************ //