applications/utilities/preProcessing/molConfig/createMolecules.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*----------------------------------------------------------------------------*/
  26
  27 #include "molConfig.H"
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 void Foam::molConfig::createMolecules()
  32 {
  33     Info<< nl << "Creating molecules from zone specifications\n" << endl;
  34
  35     DynamicList<vector> initialPositions(0);
  36
  37     DynamicList<label> initialIds(0);
  38
  39     DynamicList<scalar> initialMasses(0);
  40
  41     DynamicList<label> initialCelli(0);
  42
  43     DynamicList<vector> initialVelocities(0);
  44
  45     DynamicList<vector> initialAccelerations(0);
  46
  47     DynamicList<label> initialTethered(0);
  48
  49     DynamicList<vector> initialTetherPositions(0);
  50
  51     label totalMols = 0;
  52
  53     label idAssign;
  54
  55     Random rand(clock::getTime());
  56
  57 // * * * * * * * * Building the IdList * * * * * * * * * //
  58
  59 //Notes: - each processor will have an identical idList_.
  60 //       - The order of id's inside the idList_ depends on the order
  61 //         of subDicts inside the molConigDict.
  62
  63     Info<< "Building the idList: " ;
  64
  65     forAll(molConfigDescription_.toc(), cZs)
  66     {
  67         word subDictName (molConfigDescription_.toc()[cZs]);
  68
  69         word iD (molConfigDescription_.subDict(subDictName).lookup("id"));
  70
  71         if (findIndex(idList_,iD) == -1)
  72         {
  73             idList_.append(iD);
  74         }
  75     }
  76
  77     forAll(idList_, i)
  78     {
  79         Info << " " << idList_[i];
  80     }
  81
  82     Info << nl << endl;
  83
  84 // * * * * * * * * Filling the Mesh * * * * * * * * * //
  85
  86     const cellZoneMesh& cellZoneI = mesh_.cellZones();
  87
  88     if (cellZoneI.size())
  89     {
  90         Info<< "Filling the zones with molecules." << nl << endl;
  91     }
  92     else
  93     {
  94         FatalErrorIn("void createMolecules()\n")
  95             << "No cellZones found in mesh description."
  96             << abort(FatalError);
  97     }
  98
  99     forAll (cellZoneI, cZ)
 100     {
 101         if (cellZoneI[cZ].size())
 102         {
 103             if (!molConfigDescription_.found(cellZoneI[cZ].name()))
 104             {
 105                 Info << "Zone specification subDictionary: "
 106                     << cellZoneI[cZ].name() << " not found." << endl;
 107             }
 108             else
 109             {
 110                 label n = 0;
 111
 112                 label totalZoneMols = 0;
 113
 114                 label molsPlacedThisIteration;
 115
 116 #               include "readZoneSubDict.H"
 117
 118                 idAssign = findIndex(idList_,id);
 119
 120 #               include "startingPoint.H"
 121
 122                 // Continue trying to place molecules as long as at
 123                 // least one molecule is placed in each iteration.
 124                 // The "|| totalZoneMols == 0" condition means that the
 125                 // algorithm will continue if the origin is outside the
 126                 // zone - it will cause an infinite loop if no molecules
 127                 // are ever placed by the algorithm.
 128
 129                 if (latticeStructure != "empty")
 130                 {
 131
 132                     while
 133                     (
 134                         molsPlacedThisIteration != 0
 135                      || totalZoneMols == 0
 136                     )
 137                     {
 138                         molsPlacedThisIteration = 0;
 139
 140                         bool partOfLayerInBounds = false;
 141
 142 #                       include "createPositions.H"
 143
 144                         if
 145                         (
 146                             totalZoneMols == 0
 147                          && !partOfLayerInBounds
 148                         )
 149                         {
 150                             WarningIn("molConfig::createMolecules()")
 151                                 << "A whole layer of unit cells was placed "
 152                                 << "outside the bounds of the mesh, but no "
 153                                 << "molecules have been placed in zone '"
 154                                 << cellZoneI[cZ].name()
 155                                 << "'.  This is likely to be because the zone "
 156                                 << "has few cells ("
 157                                 << cellZoneI[cZ].size()
 158                                 << " in this case) and no lattice position "
 159                                 << "fell inside them.  "
 160                                 << "Aborting filling this zone."
 161                                 << endl;
 162
 163                             break;
 164                         }
 165
 166                         totalZoneMols += molsPlacedThisIteration;
 167
 168                         n++;
 169                     }
 170
 171                     label molN;
 172
 173                     for
 174                     (
 175                         molN = totalMols;
 176                         molN < totalMols + totalZoneMols;
 177                         molN++
 178                     )
 179                     {
 180                         initialIds.append(idAssign);
 181
 182                         initialMasses.append(mass);
 183
 184                         initialAccelerations.append(vector::zero);
 185
 186                         if (tethered)
 187                         {
 188                             initialTethered.append(1);
 189
 190                             initialTetherPositions.append
 191                             (
 192                                 initialPositions[molN]
 193                             );
 194                         }
 195
 196                         else
 197                         {
 198                             initialTethered.append(0);
 199
 200                             initialTetherPositions.append(vector::zero);
 201                         }
 202                     }
 203
 204 #                   include "createVelocities.H"
 205
 206 #                   include "correctVelocities.H"
 207
 208                 }
 209
 210                 totalMols += totalZoneMols;
 211             }
 212         }
 213     }
 214
 215     idList_.shrink();
 216
 217     positions_ = initialPositions;
 218
 219     positions_.setSize(initialPositions.size());
 220
 221     id_ = initialIds;
 222
 223     id_.setSize(initialIds.size());
 224
 225     mass_ = initialMasses;
 226
 227     mass_.setSize(initialMasses.size());
 228
 229     cells_ = initialCelli;
 230
 231     cells_.setSize(initialCelli.size());
 232
 233     U_ = initialVelocities;
 234
 235     U_.setSize(initialVelocities.size());
 236
 237     A_ = initialAccelerations;
 238
 239     A_.setSize(initialAccelerations.size());
 240
 241     tethered_ = initialTethered;
 242
 243     tethered_.setSize(initialTethered.size());
 244
 245     tetherPositions_ = initialTetherPositions;
 246
 247     tetherPositions_.setSize(initialTetherPositions.size());
 248
 249     nMol_ = totalMols;
 250 }
 251
 252
 253 // ************************************************************************* //