pairPotentialList.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
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 "OFstream.H"
28 #include "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  {
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  {
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  {
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  fileHandler().mkDir(mesh.time().path());
120  autoPtr<Ostream> ppTabFile
121  (
122  fileHandler().NewOFstream
123  (
124  mesh.time().path()/pairPotentialName
125  )
126  );
127 
128  if
129  (
130  !(*this)[pairPotentialIndex(a, b)].writeEnergyAndForceTables
131  (
132  ppTabFile()
133  )
134  )
135  {
137  << "Failed writing to "
138  << ppTabFile().name() << nl
139  << abort(FatalError);
140  }
141  }
142  }
143  }
144 
145  if (!pairPotentialDict.found("electrostatic"))
146  {
148  << "Pair pairPotential specification subDict electrostatic"
149  << nl << abort(FatalError);
150  }
151 
152  electrostaticPotential_ = pairPotential::New
153  (
154  "electrostatic",
155  pairPotentialDict.subDict("electrostatic")
156  );
157 
158  if (electrostaticPotential_->rCut() > rCutMax_)
159  {
160  rCutMax_ = electrostaticPotential_->rCut();
161  }
162 
163  if (electrostaticPotential_->writeTables())
164  {
165  fileHandler().mkDir(mesh.time().path());
166  autoPtr<Ostream> ppTabFile
167  (
168  fileHandler().NewOFstream
169  (
170  mesh.time().path()/"electrostatic"
171  )
172  );
173 
174  if (!electrostaticPotential_->writeEnergyAndForceTables(ppTabFile()))
175  {
177  << "Failed writing to "
178  << ppTabFile().name() << nl
179  << abort(FatalError);
180  }
181  }
182 
183  rCutMaxSqr_ = rCutMax_*rCutMax_;
184 }
185 
186 
187 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
188 
190 :
192 {}
193 
194 
196 (
197  const List<word>& idList,
198  const dictionary& pairPotentialDict,
199  const polyMesh& mesh
200 )
201 :
203 {
204  buildPotentials(idList, pairPotentialDict, mesh);
205 }
206 
207 
208 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
209 
211 {}
212 
213 
214 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
215 
217 (
218  const List<word>& idList,
219  const dictionary& pairPotentialDict,
220  const polyMesh& mesh
221 )
222 {
223  setSize(((idList.size()*(idList.size() + 1))/2));
224 
225  nIds_ = idList.size();
226 
227  readPairPotentialDict(idList, pairPotentialDict, mesh);
228 }
229 
230 
232 (
233  const label a,
234  const label b
235 ) const
236 {
237  return (*this)[pairPotentialIndex(a, b)];
238 }
239 
240 
241 bool Foam::pairPotentialList::rCutMaxSqr(const scalar rIJMagSqr) const
242 {
243  if (rIJMagSqr < rCutMaxSqr_)
244  {
245  return true;
246  }
247  else
248  {
249  return false;
250  }
251 }
252 
253 
255 (
256  const label a,
257  const label b,
258  const scalar rIJMagSqr
259 ) const
260 {
261  if (rIJMagSqr < rCutSqr(a, b))
262  {
263  return true;
264  }
265  else
266  {
267  return false;
268  }
269 }
270 
271 
273 (
274  const label a,
275  const label b
276 ) const
277 {
278  return (*this)[pairPotentialIndex(a, b)].rMin();
279 }
280 
281 
282 Foam::scalar Foam::pairPotentialList::dr
283 (
284  const label a,
285  const label b
286 ) const
287 {
288  return (*this)[pairPotentialIndex(a, b)].dr();
289 }
290 
291 
293 (
294  const label a,
295  const label b
296 ) const
297 {
298  return (*this)[pairPotentialIndex(a, b)].rCutSqr();
299 }
300 
301 
303 (
304  const label a,
305  const label b
306 ) const
307 {
308  return (*this)[pairPotentialIndex(a, b)].rCut();
309 }
310 
311 
313 (
314  const label a,
315  const label b,
316  const scalar rIJMag
317 ) const
318 {
319  scalar f = (*this)[pairPotentialIndex(a, b)].force(rIJMag);
320 
321  return f;
322 }
323 
324 
326 (
327  const label a,
328  const label b,
329  const scalar rIJMag
330 ) const
331 {
332  scalar e = (*this)[pairPotentialIndex(a, b)].energy(rIJMag);
333 
334  return e;
335 }
336 
337 
338 // ************************************************************************* //
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
static autoPtr< pairPotential > New(const word &name, const dictionary &pairPotentialProperties)
Return a reference to the selected viscosity model.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
virtual bool mkDir(const fileName &, mode_t=0777) const =0
Make directory.
bool rCutSqr(const label a, const label b, const scalar rIJMagSqr) const
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const pairPotential & pairPotentialFunction(const label a, const label b) const
const fileOperation & fileHandler()
Get current file handler.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:131
static const char nl
Definition: Ostream.H:260
labelList f(nPoints)
void buildPotentials(const List< word > &idList, const dictionary &pairPotentialDict, const polyMesh &mesh)
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
scalar rMin(const label a, const label b) const
messageStream Info
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
scalar dr(const label a, const label b) const
scalar rCut(const label a, const label b) const
~pairPotentialList()
Destructor.
scalar force(const label a, const label b, const scalar rIJMag) const
scalar energy(const label a, const label b, const scalar rIJMag) const