potential.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) 2011-2016 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 "potential.H"
27 
28 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
29 
30 void Foam::potential::setSiteIdList(const dictionary& moleculePropertiesDict)
31 {
32  DynamicList<word> siteIdList;
33  DynamicList<word> pairPotentialSiteIdList;
34 
35  forAll(idList_, i)
36  {
37  const word& id(idList_[i]);
38 
39  if (!moleculePropertiesDict.found(id))
40  {
42  << id << " molecule subDict not found"
43  << nl << abort(FatalError);
44  }
45 
46  const dictionary& molDict(moleculePropertiesDict.subDict(id));
47 
48  List<word> siteIdNames = molDict.lookup("siteIds");
49 
50  forAll(siteIdNames, sI)
51  {
52  const word& siteId = siteIdNames[sI];
53 
54  if (findIndex(siteIdList, siteId) == -1)
55  {
56  siteIdList.append(siteId);
57  }
58  }
59 
60  List<word> pairPotSiteIds = molDict.lookup("pairPotentialSiteIds");
61 
62  forAll(pairPotSiteIds, sI)
63  {
64  const word& siteId = pairPotSiteIds[sI];
65 
66  if (findIndex(siteIdNames, siteId) == -1)
67  {
69  << siteId << " in pairPotentialSiteIds is not in siteIds: "
70  << siteIdNames << nl << abort(FatalError);
71  }
72 
73  if (findIndex(pairPotentialSiteIdList, siteId) == -1)
74  {
75  pairPotentialSiteIdList.append(siteId);
76  }
77  }
78  }
79 
80  nPairPotIds_ = pairPotentialSiteIdList.size();
81 
82  forAll(siteIdList, aSIN)
83  {
84  const word& siteId = siteIdList[aSIN];
85 
86  if (findIndex(pairPotentialSiteIdList, siteId) == -1)
87  {
88  pairPotentialSiteIdList.append(siteId);
89  }
90  }
91 
92  siteIdList_.transfer(pairPotentialSiteIdList.shrink());
93 }
94 
95 
96 void Foam::potential::potential::readPotentialDict()
97 {
98  Info<< nl << "Reading potential dictionary:" << endl;
99 
100  IOdictionary idListDict
101  (
102  IOobject
103  (
104  "idList",
105  mesh_.time().constant(),
106  mesh_,
109  )
110  );
111 
112  idList_ = List<word>(idListDict.lookup("idList"));
113 
114  setSiteIdList
115  (
116  IOdictionary
117  (
118  IOobject
119  (
120  "moleculeProperties",
121  mesh_.time().constant(),
122  mesh_,
125  false
126  )
127  )
128  );
129 
130  List<word> pairPotentialSiteIdList
131  (
132  SubList<word>(siteIdList_, nPairPotIds_)
133  );
134 
135  Info<< nl << "Unique site ids found: " << siteIdList_
136  << nl << "Site Ids requiring a pair potential: "
137  << pairPotentialSiteIdList
138  << endl;
139 
140  List<word> tetherSiteIdList(0);
141 
142  if (idListDict.found("tetherSiteIdList"))
143  {
144  tetherSiteIdList = List<word>(idListDict.lookup("tetherSiteIdList"));
145  }
146 
147  IOdictionary potentialDict
148  (
149  IOobject
150  (
151  "potentialDict",
152  mesh_.time().system(),
153  mesh_,
156  )
157  );
158 
159  potentialEnergyLimit_ = readScalar
160  (
161  potentialDict.lookup("potentialEnergyLimit")
162  );
163 
164  if (potentialDict.found("removalOrder"))
165  {
166  List<word> remOrd = potentialDict.lookup("removalOrder");
167 
168  removalOrder_.setSize(remOrd.size());
169 
170  forAll(removalOrder_, rO)
171  {
172  removalOrder_[rO] = findIndex(idList_, remOrd[rO]);
173 
174  if (removalOrder_[rO] == -1)
175  {
177  << "removalOrder entry: " << remOrd[rO]
178  << " not found in idList."
179  << nl << abort(FatalError);
180  }
181  }
182  }
183 
184  // *************************************************************************
185  // Pair potentials
186 
187  if (!potentialDict.found("pair"))
188  {
190  << "pair potential specification subDict not found"
191  << abort(FatalError);
192  }
193 
194  const dictionary& pairDict = potentialDict.subDict("pair");
195 
196  pairPotentials_.buildPotentials
197  (
198  pairPotentialSiteIdList,
199  pairDict,
200  mesh_
201  );
202 
203  // *************************************************************************
204  // Tether potentials
205 
206  if (tetherSiteIdList.size())
207  {
208  if (!potentialDict.found("tether"))
209  {
211  << "tether potential specification subDict not found"
212  << abort(FatalError);
213  }
214 
215  const dictionary& tetherDict = potentialDict.subDict("tether");
216 
217  tetherPotentials_.buildPotentials
218  (
219  siteIdList_,
220  tetherDict,
221  tetherSiteIdList
222  );
223  }
224 
225  // *************************************************************************
226  // External Forces
227 
228  gravity_ = Zero;
229 
230  if (potentialDict.found("external"))
231  {
232  Info<< nl << "Reading external forces:" << endl;
233 
234  const dictionary& externalDict = potentialDict.subDict("external");
235 
236  // gravity
237  externalDict.readIfPresent("gravity", gravity_);
238  }
239 
240  Info<< nl << tab << "gravity = " << gravity_ << endl;
241 }
242 
243 
244 void Foam::potential::potential::readMdInitialiseDict
245 (
246  const IOdictionary& mdInitialiseDict,
247  IOdictionary& idListDict
248 )
249 {
250  IOdictionary moleculePropertiesDict
251  (
252  IOobject
253  (
254  "moleculeProperties",
255  mesh_.time().constant(),
256  mesh_,
259  false
260  )
261  );
262 
263  DynamicList<word> idList;
264 
265  DynamicList<word> tetherSiteIdList;
266 
267  forAll(mdInitialiseDict.toc(), zone)
268  {
269  const dictionary& zoneDict = mdInitialiseDict.subDict
270  (
271  mdInitialiseDict.toc()[zone]
272  );
273 
274  List<word> latticeIds
275  (
276  zoneDict.lookup("latticeIds")
277  );
278 
279  forAll(latticeIds, i)
280  {
281  const word& id = latticeIds[i];
282 
283  if (!moleculePropertiesDict.found(id))
284  {
286  << "Molecule type " << id
287  << " not found in moleculeProperties dictionary." << nl
288  << abort(FatalError);
289  }
290 
291  if (findIndex(idList,id) == -1)
292  {
293  idList.append(id);
294  }
295  }
296 
297  List<word> tetherSiteIds
298  (
299  zoneDict.lookup("tetherSiteIds")
300  );
301 
302  forAll(tetherSiteIds, t)
303  {
304  const word& tetherSiteId = tetherSiteIds[t];
305 
306  bool idFound = false;
307 
308  forAll(latticeIds, i)
309  {
310  if (idFound)
311  {
312  break;
313  }
314 
315  const word& id = latticeIds[i];
316 
317  List<word> siteIds
318  (
319  moleculePropertiesDict.subDict(id).lookup("siteIds")
320  );
321 
322  if (findIndex(siteIds, tetherSiteId) != -1)
323  {
324  idFound = true;
325  }
326  }
327 
328  if (idFound)
329  {
330  tetherSiteIdList.append(tetherSiteId);
331  }
332  else
333  {
335  << " not found as a site of any molecule in zone." << nl
336  << abort(FatalError);
337  }
338  }
339  }
340 
341  idList_.transfer(idList);
342 
343  tetherSiteIdList.shrink();
344 
345  idListDict.add("idList", idList_);
346 
347  idListDict.add("tetherSiteIdList", tetherSiteIdList);
348 
349  setSiteIdList(moleculePropertiesDict);
350 }
351 
352 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
353 
354 Foam::potential::potential(const polyMesh& mesh)
355 :
356  mesh_(mesh)
357 {
358  readPotentialDict();
359 }
360 
361 
362 Foam::potential::potential
363 (
364  const polyMesh& mesh,
365  const IOdictionary& mdInitialiseDict,
366  IOdictionary& idListDict
367 )
368 :
369  mesh_(mesh)
370 {
371  readMdInitialiseDict(mdInitialiseDict, idListDict);
372 }
373 
374 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
375 
377 {}
378 
379 
380 // ************************************************************************* //
const Time & time() const
Return time.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
static const char tab
Definition: Ostream.H:261
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const List< word > & siteIdList() const
Definition: potentialI.H:40
~potential()
Destructor.
Definition: potential.C:376
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
const word & constant() const
Return constant name.
Definition: TimePaths.H:124
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:306
static const zero Zero
Definition: zero.H:91
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static const char nl
Definition: Ostream.H:262
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
void buildPotentials(const List< word > &siteIdList, const dictionary &tetherPotentialDict, const List< word > &tetherSiteIdList)
void setSize(const label)
Reset size of List.
Definition: List.C:295
void buildPotentials(const List< word > &idList, const dictionary &pairPotentialDict, const polyMesh &mesh)
messageStream Info
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
const List< word > & idList() const
Definition: potentialI.H:34
const word & system() const
Return system name.
Definition: TimePaths.H:114
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:365
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451