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 
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 // ************************************************************************* //
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: PtrList.H:75
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
virtual bool mkDir(const fileName &, mode_t=0777) const =0
Make directory.
scalar energy(const label a, const label b, const scalar rIJMag) const
scalar force(const label a, const label b, const scalar rIJMag) const
scalar rCut(const label a, const label b) const
bool rCutSqr(const label a, const label b, const scalar rIJMagSqr) const
scalar dr(const label a, const label b) const
void buildPotentials(const List< word > &idList, const dictionary &pairPotentialDict, const polyMesh &mesh)
const pairPotential & pairPotentialFunction(const label a, const label b) const
scalar rMin(const label a, const label b) const
static autoPtr< pairPotential > New(const word &name, const dictionary &pairPotentialProperties)
Return a reference to the selected viscosity model.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
volScalarField & b
Definition: createFields.H:25
const fileOperation & fileHandler()
Get current file handler.
const doubleScalar e
Definition: doubleScalar.H:106
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
errorManip< error > abort(error &err)
Definition: errorManip.H:131
messageStream Info
error FatalError
static const char nl
Definition: Ostream.H:266
points setSize(newPointi)
labelList f(nPoints)