FreeFOAM The Cross-Platform CFD Toolkit
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) 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 "potential.H"
27 
28 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
29 
30 void Foam::potential::setSiteIdList(const IOdictionary& moleculePropertiesDict)
31 {
32  DynamicList<word> siteIdList;
33 
34  DynamicList<word> pairPotentialSiteIdList;
35 
36  forAll(idList_, i)
37  {
38  const word& id(idList_[i]);
39 
40  if (!moleculePropertiesDict.found(id))
41  {
42  FatalErrorIn("potential.C") << nl
43  << id << " molecule subDict not found"
44  << nl << abort(FatalError);
45  }
46 
47  const dictionary& molDict(moleculePropertiesDict.subDict(id));
48 
49  List<word> siteIdNames = molDict.lookup("siteIds");
50 
51  forAll(siteIdNames, sI)
52  {
53  const word& siteId = siteIdNames[sI];
54 
55  if(findIndex(siteIdList, siteId) == -1)
56  {
57  siteIdList.append(siteId);
58  }
59  }
60 
61  List<word> pairPotSiteIds = molDict.lookup("pairPotentialSiteIds");
62 
63  forAll(pairPotSiteIds, sI)
64  {
65  const word& siteId = pairPotSiteIds[sI];
66 
67  if(findIndex(siteIdNames, siteId) == -1)
68  {
69  FatalErrorIn("potential.C") << nl
70  << siteId << " in pairPotentialSiteIds is not in siteIds: "
71  << siteIdNames << nl << abort(FatalError);
72  }
73 
74  if(findIndex(pairPotentialSiteIdList, siteId) == -1)
75  {
76  pairPotentialSiteIdList.append(siteId);
77  }
78  }
79  }
80 
81  nPairPotIds_ = pairPotentialSiteIdList.size();
82 
83  forAll(siteIdList, aSIN)
84  {
85  const word& siteId = siteIdList[aSIN];
86 
87  if(findIndex(pairPotentialSiteIdList, siteId) == -1)
88  {
89  pairPotentialSiteIdList.append(siteId);
90  }
91  }
92 
93  siteIdList_.transfer(pairPotentialSiteIdList.shrink());
94 }
95 
96 
97 void Foam::potential::potential::readPotentialDict()
98 {
99  Info<< nl << "Reading potential dictionary:" << endl;
100 
101  IOdictionary idListDict
102  (
103  IOobject
104  (
105  "idList",
106  mesh_.time().constant(),
107  mesh_,
110  )
111  );
112 
113  idList_ = List<word>(idListDict.lookup("idList"));
114 
115  IOdictionary moleculePropertiesDict
116  (
117  IOobject
118  (
119  "moleculeProperties",
120  mesh_.time().constant(),
121  mesh_,
124  false
125  )
126  );
127 
128  setSiteIdList(moleculePropertiesDict);
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  {
176  FatalErrorIn("potentials.C") << nl
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  {
189  FatalErrorIn("potentials.C") << nl
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  {
210  FatalErrorIn("potential.C") << nl
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_ = vector::zero;
229 
230  if (potentialDict.found("external"))
231  {
232 
233  Info << nl << "Reading external forces:" << endl;
234 
235  const dictionary& externalDict = potentialDict.subDict("external");
236 
237  // *********************************************************************
238  // gravity
239 
240  if (externalDict.found("gravity"))
241  {
242  gravity_ = externalDict.lookup("gravity");
243  }
244  }
245 
246  Info << nl << tab << "gravity = " << gravity_ << endl;
247 }
248 
249 
250 void Foam::potential::potential::readMdInitialiseDict
251 (
252  const IOdictionary& mdInitialiseDict,
253  IOdictionary& idListDict
254 )
255 {
256  IOdictionary moleculePropertiesDict
257  (
258  IOobject
259  (
260  "moleculeProperties",
261  mesh_.time().constant(),
262  mesh_,
265  false
266  )
267  );
268 
269  DynamicList<word> idList;
270 
271  DynamicList<word> tetherSiteIdList;
272 
273  forAll(mdInitialiseDict.toc(), zone)
274  {
275  const dictionary& zoneDict = mdInitialiseDict.subDict
276  (
277  mdInitialiseDict.toc()[zone]
278  );
279 
280  List<word> latticeIds
281  (
282  zoneDict.lookup("latticeIds")
283  );
284 
285  forAll(latticeIds, i)
286  {
287  const word& id = latticeIds[i];
288 
289  if(!moleculePropertiesDict.found(id))
290  {
291  FatalErrorIn("potential.C") << nl
292  << "Molecule type "
293  << id
294  << " not found in moleculeProperties dictionary."
295  << nl
296  << abort(FatalError);
297  }
298 
299  if (findIndex(idList,id) == -1)
300  {
301  idList.append(id);
302  }
303  }
304 
305  List<word> tetherSiteIds
306  (
307  zoneDict.lookup("tetherSiteIds")
308  );
309 
310  forAll(tetherSiteIds, t)
311  {
312  const word& tetherSiteId = tetherSiteIds[t];
313 
314  bool idFound = false;
315 
316  forAll(latticeIds, i)
317  {
318  if (idFound)
319  {
320  break;
321  }
322 
323  const word& id = latticeIds[i];
324 
325  List<word> siteIds
326  (
327  moleculePropertiesDict.subDict(id).lookup("siteIds")
328  );
329 
330  if (findIndex(siteIds, tetherSiteId) != -1)
331  {
332  idFound = true;
333  }
334  }
335 
336  if (idFound)
337  {
338  tetherSiteIdList.append(tetherSiteId);
339  }
340  else
341  {
342  FatalErrorIn("potential.C") << nl
343  << "Tether id "
344  << tetherSiteId
345  << " not found as a site of any molecule in zone."
346  << nl
347  << abort(FatalError);
348  }
349  }
350  }
351 
352  idList_.transfer(idList);
353 
354  tetherSiteIdList.shrink();
355 
356  idListDict.add("idList", idList_);
357 
358  idListDict.add("tetherSiteIdList", tetherSiteIdList);
359 
360  setSiteIdList(moleculePropertiesDict);
361 }
362 
363 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
364 
365 Foam::potential::potential(const polyMesh& mesh)
366 :
367  mesh_(mesh)
368 {
369  readPotentialDict();
370 }
371 
372 
373 Foam::potential::potential
374 (
375  const polyMesh& mesh,
376  const IOdictionary& mdInitialiseDict,
377  IOdictionary& idListDict
378 )
379 :
380  mesh_(mesh)
381 {
382  readMdInitialiseDict(mdInitialiseDict, idListDict);
383 }
384 
385 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
386 
388 {}
389 
390 
391 // ************************ vim: set sw=4 sts=4 et: ************************ //