diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C
deleted file mode 100644
index c6d54036f7..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C
+++ /dev/null
@@ -1,1265 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see
-void Foam::MoleculeCloud::buildConstProps()
-{
- const List& idList(pot_.idList());
-
- constPropList_.setSize(idList.size());
-
- const List& siteIdList(pot_.siteIdList());
-
- IOdictionary molPropertiesDict
- (
- IOobject
- (
- this->name() + "Properties",
- mesh_.time().constant(),
- mesh_,
- IOobject::MUST_READ_IF_MODIFIED,
- IOobject::NO_WRITE,
- false
- )
- );
-
- Info<< nl << "Reading " << molPropertiesDict.name() << endl;
-
- forAll(idList, i)
- {
- const word& id = idList[i];
-
- const dictionary& molDict(molPropertiesDict.subDict(id));
-
- List siteIdNames = molDict.lookup("siteIds");
-
- List siteIds(siteIdNames.size());
-
- forAll(siteIdNames, sI)
- {
- const word& siteId = siteIdNames[sI];
-
- siteIds[sI] = findIndex(siteIdList, siteId);
-
- if (siteIds[sI] == -1)
- {
- FatalErrorIn
- (
- "void Foam::MoleculeCloud::buildConstProps()"
- )
- << siteId << " site not found."
- << nl << abort(FatalError);
- }
- }
-
- constPropList_[i] = typename MoleculeType::constantProperties
- (
- molDict,
- siteIds
- );
- }
-}
-
-
-template
-void Foam::MoleculeCloud::setSiteSizesAndPositions()
-{
- forAllIter(typename MoleculeCloud, *this, mol)
- {
- const typename MoleculeType::constantProperties& cP
- (
- constProps(mol().id())
- );
-
- mol().setSiteSizes(cP.nSites());
-
- mol().setSitePositions(cP);
- }
-}
-
-
-template
-void Foam::MoleculeCloud::buildCellOccupancy()
-{
- forAll(cellOccupancy_, cO)
- {
- cellOccupancy_[cO].clear();
- }
-
- forAllIter(typename MoleculeCloud, *this, mol)
- {
- cellOccupancy_[mol().cell()].append(&mol());
- }
-
- forAll(cellOccupancy_, cO)
- {
- cellOccupancy_[cO].shrink();
- }
-}
-
-
-template
-void Foam::MoleculeCloud::calculatePairForce()
-{
- PstreamBuffers pBufs(Pstream::nonBlocking);
-
- // Start sending referred data
- il_.sendReferredData(cellOccupancy(), pBufs);
-
- MoleculeType* molI = NULL;
- MoleculeType* molJ = NULL;
-
- {
- // Real-Real interactions
-
- const labelListList& dil = il_.dil();
-
- forAll(dil, d)
- {
- forAll(cellOccupancy_[d],cellIMols)
- {
- molI = cellOccupancy_[d][cellIMols];
-
- forAll(dil[d], interactingCells)
- {
- List cellJ =
- cellOccupancy_[dil[d][interactingCells]];
-
- forAll(cellJ, cellJMols)
- {
- molJ = cellJ[cellJMols];
-
- evaluatePair(*molI, *molJ);
- }
- }
-
- forAll(cellOccupancy_[d], cellIOtherMols)
- {
- molJ = cellOccupancy_[d][cellIOtherMols];
-
- if (molJ > molI)
- {
- evaluatePair(*molI, *molJ);
- }
- }
- }
- }
- }
-
- // Receive referred data
- il_.receiveReferredData(pBufs);
-
- {
- // Real-Referred interactions
-
- const labelListList& ril = il_.ril();
-
- List >& referredMols = il_.referredParticles();
-
- forAll(ril, r)
- {
- const List& realCells = ril[r];
-
- IDLList& refMols = referredMols[r];
-
- forAllIter
- (
- typename IDLList,
- refMols,
- refMol
- )
- {
- forAll(realCells, rC)
- {
- List cellI = cellOccupancy_[realCells[rC]];
-
- forAll(cellI, cellIMols)
- {
- molI = cellI[cellIMols];
-
- evaluatePair(*molI, refMol());
- }
- }
- }
- }
- }
-}
-
-
-template
-void Foam::MoleculeCloud::calculateTetherForce()
-{
- const tetherPotentialList& tetherPot(pot_.tetherPotentials());
-
- forAllIter(typename MoleculeCloud, *this, mol)
- {
- if (mol().tethered())
- {
- vector rIT = mol().position() - mol().specialPosition();
-
- label idI = mol().id();
-
- scalar massI = constProps(idI).mass();
-
- vector fIT = tetherPot.force(idI, rIT);
-
- mol().a() += fIT/massI;
-
- mol().potentialEnergy() += tetherPot.energy(idI, rIT);
-
- // What to do here?
- // mol().rf() += rIT*fIT;
- }
- }
-}
-
-
-template
-void Foam::MoleculeCloud::calculateExternalForce()
-{
- forAllIter(typename MoleculeCloud, *this, mol)
- {
- mol().a() += pot_.gravity();
- }
-}
-
-
-template
-void Foam::MoleculeCloud::removeHighEnergyOverlaps()
-{
- Info<< nl << "Removing high energy overlaps, limit = "
- << pot_.potentialEnergyLimit()
- << nl << "Removal order:";
-
- forAll(pot_.removalOrder(), rO)
- {
- Info<< ' ' << pot_.idList()[pot_.removalOrder()[rO]];
- }
-
- Info<< nl ;
-
- label initialSize = this->size();
-
- buildCellOccupancy();
-
- // Real-Real interaction
-
- MoleculeType* molI = NULL;
- MoleculeType* molJ = NULL;
-
- {
- DynamicList molsToDelete;
-
- const labelListList& dil(il_.dil());
-
- forAll(dil, d)
- {
- forAll(cellOccupancy_[d], cellIMols)
- {
- molI = cellOccupancy_[d][cellIMols];
-
- forAll(dil[d], interactingCells)
- {
- List cellJ =
- cellOccupancy_[dil[d][interactingCells]];
-
- forAll(cellJ, cellJMols)
- {
- molJ = cellJ[cellJMols];
-
- if (evaluatePotentialLimit(*molI, *molJ))
- {
- label idI = molI->id();
-
- label idJ = molJ->id();
-
- if
- (
- idI == idJ
- || findIndex(pot_.removalOrder(), idJ)
- < findIndex(pot_.removalOrder(), idI)
- )
- {
- if (findIndex(molsToDelete, molJ) == -1)
- {
- molsToDelete.append(molJ);
- }
- }
- else if (findIndex(molsToDelete, molI) == -1)
- {
- molsToDelete.append(molI);
- }
- }
- }
- }
-
- forAll(cellOccupancy_[d], cellIOtherMols)
- {
- molJ = cellOccupancy_[d][cellIOtherMols];
-
- if (molJ > molI)
- {
- if (evaluatePotentialLimit(*molI, *molJ))
- {
- label idI = molI->id();
-
- label idJ = molJ->id();
-
- if
- (
- idI == idJ
- || findIndex(pot_.removalOrder(), idJ)
- < findIndex(pot_.removalOrder(), idI)
- )
- {
- if (findIndex(molsToDelete, molJ) == -1)
- {
- molsToDelete.append(molJ);
- }
- }
- else if (findIndex(molsToDelete, molI) == -1)
- {
- molsToDelete.append(molI);
- }
- }
- }
- }
- }
- }
-
- forAll(molsToDelete, mTD)
- {
- deleteParticle(*(molsToDelete[mTD]));
- }
- }
-
- buildCellOccupancy();
-
- PstreamBuffers pBufs(Pstream::nonBlocking);
-
- // Start sending referred data
- il_.sendReferredData(cellOccupancy(), pBufs);
-
- // Receive referred data
- il_.receiveReferredData(pBufs);
-
- // Real-Referred interaction
-
- {
- DynamicList molsToDelete;
-
- const labelListList& ril(il_.ril());
-
- List >& referredMols = il_.referredParticles();
-
- forAll(ril, r)
- {
- IDLList& refMols = referredMols[r];
-
- forAllIter
- (
- typename IDLList,
- refMols,
- refMol
- )
- {
- molJ = &refMol();
-
- const List& realCells = ril[r];
-
- forAll(realCells, rC)
- {
- label cellI = realCells[rC];
-
- List cellIMols = cellOccupancy_[cellI];
-
- forAll(cellIMols, cIM)
- {
- molI = cellIMols[cIM];
-
- if (evaluatePotentialLimit(*molI, *molJ))
- {
- label idI = molI->id();
-
- label idJ = molJ->id();
-
- if
- (
- findIndex(pot_.removalOrder(), idI)
- < findIndex(pot_.removalOrder(), idJ)
- )
- {
- if (findIndex(molsToDelete, molI) == -1)
- {
- molsToDelete.append(molI);
- }
- }
- else if
- (
- findIndex(pot_.removalOrder(), idI)
- == findIndex(pot_.removalOrder(), idJ)
- )
- {
- // Remove one of the molecules
- // arbitrarily, assuring that a
- // consistent decision is made for
- // both real-referred pairs.
-
- if (molI->origId() > molJ->origId())
- {
- if (findIndex(molsToDelete, molI) == -1)
- {
- molsToDelete.append(molI);
- }
- }
- }
- }
- }
- }
- }
- }
-
- forAll(molsToDelete, mTD)
- {
- deleteParticle(*(molsToDelete[mTD]));
- }
- }
-
- buildCellOccupancy();
-
- // Start sending referred data
- il_.sendReferredData(cellOccupancy(), pBufs);
-
- // Receive referred data
- il_.receiveReferredData(pBufs);
-
- label molsRemoved = initialSize - this->size();
-
- if (Pstream::parRun())
- {
- reduce(molsRemoved, sumOp());
- }
-
- Info<< tab << molsRemoved << " molecules removed" << endl;
-}
-
-
-template
-void Foam::MoleculeCloud::initialiseMolecules
-(
- const dictionary& mdInitialiseDict
-)
-{
- Info<< nl
- << "Initialising molecules in each zone specified in "
- << mdInitialiseDict.name()
- << endl;
-
- const cellZoneMesh& cellZones = mesh_.cellZones();
-
- if (!cellZones.size())
- {
- FatalErrorIn
- (
- "void Foam::MoleculeCloud::initialiseMolecules"
- )
- << "No cellZones found in the mesh."
- << abort(FatalError);
- }
-
- forAll(cellZones, z)
- {
- const cellZone& zone(cellZones[z]);
-
- if (zone.size())
- {
- if (!mdInitialiseDict.found(zone.name()))
- {
- Info<< "No specification subDictionary for zone "
- << zone.name() << " found, skipping." << endl;
- }
- else
- {
- const dictionary& zoneDict =
- mdInitialiseDict.subDict(zone.name());
-
- const scalar temperature
- (
- readScalar(zoneDict.lookup("temperature"))
- );
-
- const vector bulkVelocity(zoneDict.lookup("bulkVelocity"));
-
- List latticeIds
- (
- zoneDict.lookup("latticeIds")
- );
-
- List latticePositions
- (
- zoneDict.lookup("latticePositions")
- );
-
- if (latticeIds.size() != latticePositions.size())
- {
- FatalErrorIn
- (
- "Foam::MoleculeCloud::initialiseMolecules"
- )
- << "latticeIds and latticePositions must be the same "
- << " size." << nl
- << abort(FatalError);
- }
-
- diagTensor latticeCellShape
- (
- zoneDict.lookup("latticeCellShape")
- );
-
- scalar latticeCellScale = 0.0;
-
- if (zoneDict.found("numberDensity"))
- {
- scalar numberDensity = readScalar
- (
- zoneDict.lookup("numberDensity")
- );
-
- if (numberDensity < VSMALL)
- {
- WarningIn
- (
- "MoleculeCloud::initialiseMolecules"
- )
- << "numberDensity too small, not filling zone "
- << zone.name() << endl;
-
- continue;
- }
-
- latticeCellScale = pow
- (
- latticeIds.size()/(det(latticeCellShape)*numberDensity),
- (1.0/3.0)
- );
- }
- else if (zoneDict.found("massDensity"))
- {
- scalar unitCellMass = 0.0;
-
- forAll(latticeIds, i)
- {
- label id = findIndex(pot_.idList(), latticeIds[i]);
-
- const typename MoleculeType::constantProperties& cP
- (
- constProps(id)
- );
-
- unitCellMass += cP.mass();
- }
-
- scalar massDensity = readScalar
- (
- zoneDict.lookup("massDensity")
- );
-
- if (massDensity < VSMALL)
- {
- WarningIn
- (
- "MoleculeCloud::initialiseMolecules"
- )
- << "massDensity too small, not filling zone "
- << zone.name() << endl;
-
- continue;
- }
-
- latticeCellScale = pow
- (
- unitCellMass/(det(latticeCellShape)*massDensity),
- (1.0/3.0)
- );
- }
- else
- {
- FatalErrorIn
- (
- "Foam::MoleculeCloud::initialiseMolecules"
- )
- << "massDensity or numberDensity not specified " << nl
- << abort(FatalError);
- }
-
- latticeCellShape *= latticeCellScale;
-
- vector anchor(zoneDict.lookup("anchor"));
-
- bool tethered = false;
-
- if (zoneDict.found("tetherSiteIds"))
- {
- tethered = bool
- (
- List(zoneDict.lookup("tetherSiteIds")).size()
- );
- }
-
- const vector orientationAngles
- (
- zoneDict.lookup("orientationAngles")
- );
-
- scalar phi(orientationAngles.x()*pi/180.0);
-
- scalar theta(orientationAngles.y()*pi/180.0);
-
- scalar psi(orientationAngles.z()*pi/180.0);
-
- const tensor R
- (
- cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
- cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
- sin(psi)*sin(theta),
- - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
- - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
- cos(psi)*sin(theta),
- sin(theta)*sin(phi),
- - sin(theta)*cos(phi),
- cos(theta)
- );
-
- // Find the optimal anchor position. Finding the approximate
- // mid-point of the zone of cells and snapping to the nearest
- // lattice location.
-
- vector zoneMin = VGREAT*vector::one;
-
- vector zoneMax = -VGREAT*vector::one;
-
- forAll(zone, cell)
- {
- const point cellCentre = mesh_.cellCentres()[zone[cell]];
-
- if (cellCentre.x() > zoneMax.x())
- {
- zoneMax.x() = cellCentre.x();
- }
- if (cellCentre.x() < zoneMin.x())
- {
- zoneMin.x() = cellCentre.x();
- }
- if (cellCentre.y() > zoneMax.y())
- {
- zoneMax.y() = cellCentre.y();
- }
- if (cellCentre.y() < zoneMin.y())
- {
- zoneMin.y() = cellCentre.y();
- }
- if (cellCentre.z() > zoneMax.z())
- {
- zoneMax.z() = cellCentre.z();
- }
- if (cellCentre.z() < zoneMin.z())
- {
- zoneMin.z() = cellCentre.z();
- }
- }
-
- point zoneMid = 0.5*(zoneMin + zoneMax);
-
- point latticeMid = inv(latticeCellShape) & (R.T()
- & (zoneMid - anchor));
-
- point latticeAnchor
- (
- label(latticeMid.x() + 0.5*sign(latticeMid.x())),
- label(latticeMid.y() + 0.5*sign(latticeMid.y())),
- label(latticeMid.z() + 0.5*sign(latticeMid.z()))
- );
-
- anchor += (R & (latticeCellShape & latticeAnchor));
-
- // Continue trying to place molecule as long as at
- // least one molecule is placed in each iteration.
- // The "|| totalZoneMols == 0" condition means that the
- // algorithm will continue if the origin is outside the
- // zone.
-
- label n = 0;
-
- label totalZoneMols = 0;
-
- label molsPlacedThisIteration = 0;
-
- while
- (
- molsPlacedThisIteration != 0
- || totalZoneMols == 0
- )
- {
- label sizeBeforeIteration = this->size();
-
- bool partOfLayerInBounds = false;
-
- if (n == 0)
- {
- // Special treatment is required for the first position,
- // i.e. iteration zero.
-
- labelVector unitCellLatticePosition(0,0,0);
-
- forAll(latticePositions, p)
- {
- label id = findIndex(pot_.idList(), latticeIds[p]);
-
- const vector& latticePosition =
- vector
- (
- unitCellLatticePosition.x(),
- unitCellLatticePosition.y(),
- unitCellLatticePosition.z()
- )
- + latticePositions[p];
-
- point globalPosition =
- anchor
- + (R & (latticeCellShape & latticePosition));
-
- partOfLayerInBounds = mesh_.bounds().contains
- (
- globalPosition
- );
-
- label cell = -1;
- label tetFace = -1;
- label tetPt = -1;
-
- mesh_.findCellFacePt
- (
- globalPosition,
- cell,
- tetFace,
- tetPt
- );
-
- if (findIndex(zone, cell) != -1)
- {
- createMolecule
- (
- globalPosition,
- cell,
- tetFace,
- tetPt,
- id,
- tethered,
- temperature,
- bulkVelocity
- );
- }
- }
- }
- else
- {
- // Place top and bottom caps.
-
- labelVector unitCellLatticePosition(0,0,0);
-
- for
- (
- unitCellLatticePosition.z() = -n;
- unitCellLatticePosition.z() <= n;
- unitCellLatticePosition.z() += 2*n
- )
- {
- for
- (
- unitCellLatticePosition.y() = -n;
- unitCellLatticePosition.y() <= n;
- unitCellLatticePosition.y()++
- )
- {
- for
- (
- unitCellLatticePosition.x() = -n;
- unitCellLatticePosition.x() <= n;
- unitCellLatticePosition.x()++
- )
- {
- forAll(latticePositions, p)
- {
- label id = findIndex
- (
- pot_.idList(),
- latticeIds[p]
- );
-
- const vector& latticePosition =
- vector
- (
- unitCellLatticePosition.x(),
- unitCellLatticePosition.y(),
- unitCellLatticePosition.z()
- )
- + latticePositions[p];
-
- point globalPosition =
- anchor
- + (
- R
- & (
- latticeCellShape
- & latticePosition
- )
- );
-
- partOfLayerInBounds =
- mesh_.bounds().contains
- (
- globalPosition
- );
-
- label cell = -1;
- label tetFace = -1;
- label tetPt = -1;
-
- mesh_.findCellFacePt
- (
- globalPosition,
- cell,
- tetFace,
- tetPt
- );
-
- if (findIndex(zone, cell) != -1)
- {
- createMolecule
- (
- globalPosition,
- cell,
- tetFace,
- tetPt,
- id,
- tethered,
- temperature,
- bulkVelocity
- );
- }
- }
- }
- }
- }
-
- for
- (
- unitCellLatticePosition.z() = -(n-1);
- unitCellLatticePosition.z() <= (n-1);
- unitCellLatticePosition.z()++
- )
- {
- for (label iR = 0; iR <= 2*n -1; iR++)
- {
- unitCellLatticePosition.x() = n;
-
- unitCellLatticePosition.y() = -n + (iR + 1);
-
- for (label iK = 0; iK < 4; iK++)
- {
- forAll(latticePositions, p)
- {
- label id = findIndex
- (
- pot_.idList(),
- latticeIds[p]
- );
-
- const vector& latticePosition =
- vector
- (
- unitCellLatticePosition.x(),
- unitCellLatticePosition.y(),
- unitCellLatticePosition.z()
- )
- + latticePositions[p];
-
- point globalPosition =
- anchor
- + (
- R
- & (
- latticeCellShape
- & latticePosition
- )
- );
-
- partOfLayerInBounds =
- mesh_.bounds().contains
- (
- globalPosition
- );
-
- label cell = -1;
- label tetFace = -1;
- label tetPt = -1;
-
- mesh_.findCellFacePt
- (
- globalPosition,
- cell,
- tetFace,
- tetPt
- );
-
- if (findIndex(zone, cell) != -1)
- {
- createMolecule
- (
- globalPosition,
- cell,
- tetFace,
- tetPt,
- id,
- tethered,
- temperature,
- bulkVelocity
- );
- }
- }
-
- unitCellLatticePosition =
- labelVector
- (
- - unitCellLatticePosition.y(),
- unitCellLatticePosition.x(),
- unitCellLatticePosition.z()
- );
- }
- }
- }
- }
-
- if
- (
- totalZoneMols == 0
- && !partOfLayerInBounds
- )
- {
- WarningIn
- (
- "Foam::MoleculeCloud::"
- "initialiseMolecules()"
- )
- << "A whole layer of unit cells was placed "
- << "outside the bounds of the mesh, but no "
- << "molecules have been placed in zone '"
- << zone.name()
- << "'. This is likely to be because the zone "
- << "has few cells ("
- << zone.size()
- << " in this case) and no lattice position "
- << "fell inside them. "
- << "Aborting filling this zone."
- << endl;
-
- break;
- }
-
- molsPlacedThisIteration =
- this->size() - sizeBeforeIteration;
-
- totalZoneMols += molsPlacedThisIteration;
-
- n++;
- }
- }
- }
- }
-}
-
-
-template
-void Foam::MoleculeCloud::createMolecule
-(
- const point& position,
- label cell,
- label tetFace,
- label tetPt,
- label id,
- bool tethered,
- scalar temperature,
- const vector& bulkVelocity
-)
-{
- if (cell == -1)
- {
- mesh_.findCellFacePt(position, cell, tetFace, tetPt);
- }
-
- if (cell == -1)
- {
- FatalErrorIn("Foam::MoleculeCloud::createMolecule")
- << "Position specified does not correspond to a mesh cell." << nl
- << abort(FatalError);
- }
-
- point specialPosition(vector::zero);
-
- label special = 0;
-
- if (tethered)
- {
- specialPosition = position;
-
- special = MoleculeType::SPECIAL_TETHERED;
- }
-
- const typename MoleculeType::constantProperties& cP(constProps(id));
-
- typename MoleculeType::trackingData td
- (
- *this,
- MoleculeType::trackingData::tpAccess
- );
-
- addParticle
- (
- new MoleculeType
- (
- mesh_,
- position,
- cell,
- tetFace,
- tetPt,
- temperature,
- bulkVelocity,
- specialPosition,
- cP,
- td,
- special,
- id
- )
- );
-}
-
-
-template
-Foam::label Foam::MoleculeCloud::nSites() const
-{
- label n = 0;
-
- forAllConstIter(typename MoleculeCloud, *this, mol)
- {
- n += constProps(mol().id()).nSites();
- }
-
- return n;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-
-template
-Foam::MoleculeCloud::MoleculeCloud
-(
- const word& cloudName,
- const polyMesh& mesh,
- const potential& pot,
- bool readFields
-)
-:
- Cloud(mesh, cloudName, false),
- moleculeCloud(),
- mesh_(mesh),
- pot_(pot),
- cellOccupancy_(mesh_.nCells()),
- il_(mesh_, pot_.pairPotentials().rCutMax(), false),
- constPropList_(),
- rndGen_(label(971501) + 1526*Pstream::myProcNo())
-{
- if (readFields)
- {
- MoleculeType::readFields(*this);
- }
-
- buildConstProps();
-
- setSiteSizesAndPositions();
-
- removeHighEnergyOverlaps();
-
- calculateForce();
-}
-
-
-template
-Foam::MoleculeCloud::MoleculeCloud
-(
- const word& cloudName,
- const polyMesh& mesh,
- const potential& pot,
- const dictionary& mdInitialiseDict,
- bool readFields
-)
-:
- Cloud(mesh, cloudName, false),
- moleculeCloud(),
- mesh_(mesh),
- pot_(pot),
- il_(mesh_),
- constPropList_(),
- rndGen_(label(971501) + 1526*Pstream::myProcNo())
-{
- if (readFields)
- {
- MoleculeType::readFields(*this);
- }
-
- this->clear();
-
- buildConstProps();
-
- initialiseMolecules(mdInitialiseDict);
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-template
-void Foam::MoleculeCloud::evolve()
-{
- typename MoleculeType::trackingData td0
- (
- *this,
- MoleculeType::trackingData::tpFirstVelocityHalfStep
- );
- Cloud::move(td0, mesh_.time().deltaTValue());
-
- typename MoleculeType::trackingData td1
- (
- *this,
- MoleculeType::trackingData::tpLinearTrack
- );
- Cloud::move(td1, mesh_.time().deltaTValue());
-
- typename MoleculeType::trackingData td2
- (
- *this,
- MoleculeType::trackingData::tpRotationalTrack
- );
- Cloud::move(td2, mesh_.time().deltaTValue());
-
- calculateForce();
-
- typename MoleculeType::trackingData td3
- (
- *this,
- MoleculeType::trackingData::tpSecondVelocityHalfStep
- );
- Cloud::move(td3, mesh_.time().deltaTValue());
-
- info();
-}
-
-
-template
-void Foam::MoleculeCloud::calculateForce()
-{
- buildCellOccupancy();
-
- // Set accumulated quantities to zero
- forAllIter(typename MoleculeCloud, *this, mol)
- {
- mol().siteForces() = vector::zero;
-
- mol().potentialEnergy() = 0.0;
-
- mol().rf() = tensor::zero;
- }
-
- calculatePairForce();
-
- calculateTetherForce();
-
- calculateExternalForce();
-}
-
-
-template
-void Foam::MoleculeCloud::info()
-{
- // Calculates and prints the mean momentum and energy in the system
- // and the number of molecules.
-
- typename MoleculeType::trackingData td
- (
- *this,
- MoleculeType::trackingData::tpAccess
- );
-
- MoleculeType::info(td);
-}
-
-
-template
-void Foam::MoleculeCloud::writeXYZ(const fileName& fName) const
-{
- OFstream os(fName);
-
- os << nSites() << nl
- << "MoleculeCloud site positions in angstroms" << nl;
-
- forAllConstIter(typename MoleculeCloud, *this, mol)
- {
- const typename MoleculeType::constantProperties& cP
- (
- constProps(mol().id())
- );
-
- forAll(mol().sitePositions(), i)
- {
- const point& sP = mol().sitePositions()[i];
-
- os << pot_.siteIdList()[cP.sites()[i].siteId()]
- << ' ' << sP.x()*1e10
- << ' ' << sP.y()*1e10
- << ' ' << sP.z()*1e10
- << nl;
- }
- }
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H
deleted file mode 100644
index 38c912a17d..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H
+++ /dev/null
@@ -1,252 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see
-class MoleculeCloud
-:
- public Cloud,
- public moleculeCloud
-{
-
-private:
-
- // Private data
-
- //-
- const polyMesh& mesh_;
-
- //-
- const potential& pot_;
-
- //-
- List > cellOccupancy_;
-
- //-
- InteractionLists il_;
-
- //-
- List constPropList_;
-
- //-
- Random rndGen_;
-
-
- // Private Member Functions
-
- //-
- void buildConstProps();
-
- //-
- void setSiteSizesAndPositions();
-
- //- Determine which molecules are in which cells
- void buildCellOccupancy();
-
- //-
- void calculatePairForce();
-
- //-
- inline void evaluatePair
- (
- MoleculeType& molI,
- MoleculeType& molJ
- );
-
- //-
- inline bool evaluatePotentialLimit
- (
- MoleculeType& molI,
- MoleculeType& molJ
- ) const;
-
- //-
- void calculateTetherForce();
-
- //-
- void calculateExternalForce();
-
- //-
- void removeHighEnergyOverlaps();
-
- //-
- void initialiseMolecules(const dictionary& mdInitialiseDict);
-
- //-
- void createMolecule
- (
- const point& position,
- label cell,
- label tetFace,
- label tetPt,
- label id,
- bool tethered,
- scalar temperature,
- const vector& bulkVelocity
- );
-
- //-
- label nSites() const;
-
- //- Disallow default bitwise copy construct
- MoleculeCloud(const MoleculeCloud&);
-
- //- Disallow default bitwise assignment
- void operator=(const MoleculeCloud&);
-
-
-public:
-
- // Constructors
-
- //- Construct given mesh and potential references
- MoleculeCloud
- (
- const word& cloudName,
- const polyMesh& mesh,
- const potential& pot,
- bool readFields = true
- );
-
- //- Construct given mesh, potential and mdInitialiseDict
- MoleculeCloud
- (
- const word& cloudName,
- const polyMesh& mesh,
- const potential& pot,
- const dictionary& mdInitialiseDict,
- bool readFields = true
- );
-
-
- // Member Functions
-
- //- Evolve the molecules (move, calculate forces, control state etc)
- void evolve();
-
- //-
- void calculateForce();
-
- //- Print cloud information
- void info();
-
-
- // Access
-
- //-
- inline const polyMesh& mesh() const;
-
- //-
- inline const potential& pot() const;
-
- //-
- inline const List >&
- cellOccupancy() const;
-
- //-
- inline const InteractionLists& il() const;
-
- //-
- inline const List
- constProps() const;
-
- //-
- inline const typename MoleculeType::constantProperties&
- constProps(label id) const;
-
- //-
- inline Random& rndGen();
-
- //-
- inline vector equipartitionLinearVelocity
- (
- scalar temperature,
- scalar mass
- );
-
- //-
- inline vector equipartitionAngularMomentum
- (
- scalar temperature,
- const typename MoleculeType::constantProperties& cP
- );
-
-
- // Member Operators
-
- //- Write molecule sites in XYZ format
- void writeXYZ(const fileName& fName) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "MoleculeCloudI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-# include "MoleculeCloud.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H
deleted file mode 100644
index 42c506024f..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H
+++ /dev/null
@@ -1,416 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see
-inline void Foam::MoleculeCloud::evaluatePair
-(
- MoleculeType& molI,
- MoleculeType& molJ
-)
-{
- const pairPotentialList& pairPot = pot_.pairPotentials();
-
- const pairPotential& electrostatic = pairPot.electrostatic();
-
- label idI = molI.id();
-
- label idJ = molJ.id();
-
- const typename MoleculeType::constantProperties& constPropI
- (
- constProps(idI)
- );
-
- const typename MoleculeType::constantProperties& constPropJ
- (
- constProps(idJ)
- );
-
- forAll(constPropI.pairPotSites(), pI)
- {
- label sI = constPropI.pairPotSites()[pI];
-
- label idsI = constPropI.sites()[sI].siteId();
-
- forAll(constPropJ.pairPotSites(), pJ)
- {
- label sJ = constPropJ.pairPotSites()[pJ];
-
- label idsJ = constPropJ.sites()[sJ].siteId();
-
- vector rsIsJ =
- molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
-
- scalar rsIsJMagSq = magSqr(rsIsJ);
-
- if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
- {
- scalar rsIsJMag = mag(rsIsJ);
-
- vector fsIsJ =
- (rsIsJ/rsIsJMag)
- *pairPot.force(idsI, idsJ, rsIsJMag);
-
- molI.siteForces()[sI] += fsIsJ;
-
- molJ.siteForces()[sJ] += -fsIsJ;
-
- scalar potentialEnergy
- (
- pairPot.energy(idsI, idsJ, rsIsJMag)
- );
-
- molI.potentialEnergy() += 0.5*potentialEnergy;
-
- molJ.potentialEnergy() += 0.5*potentialEnergy;
-
- vector rIJ = molI.position() - molJ.position();
-
- tensor virialContribution =
- (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
-
- molI.rf() += virialContribution;
-
- molJ.rf() += virialContribution;
- }
- }
- }
-
-
- forAll(constPropI.electrostaticSites(), pI)
- {
- label sI = constPropI.electrostaticSites()[pI];
-
- forAll(constPropJ.electrostaticSites(), pJ)
- {
- label sJ = constPropJ.electrostaticSites()[pJ];
-
- vector rsIsJ =
- molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
-
- scalar rsIsJMagSq = magSqr(rsIsJ);
-
- if (rsIsJMagSq <= electrostatic.rCutSqr())
- {
- scalar rsIsJMag = mag(rsIsJ);
-
- scalar chargeI = constPropI.sites()[sI].siteCharge();
-
- scalar chargeJ = constPropJ.sites()[sJ].siteCharge();
-
- vector fsIsJ =
- (rsIsJ/rsIsJMag)
- *chargeI*chargeJ*electrostatic.force(rsIsJMag);
-
- molI.siteForces()[sI] += fsIsJ;
-
- molJ.siteForces()[sJ] += -fsIsJ;
-
- scalar potentialEnergy =
- chargeI*chargeJ
- *electrostatic.energy(rsIsJMag);
-
- molI.potentialEnergy() += 0.5*potentialEnergy;
-
- molJ.potentialEnergy() += 0.5*potentialEnergy;
-
- vector rIJ = molI.position() - molJ.position();
-
- tensor virialContribution =
- (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
-
- molI.rf() += virialContribution;
-
- molJ.rf() += virialContribution;
- }
- }
- }
-}
-
-
-template
-inline bool Foam::MoleculeCloud::evaluatePotentialLimit
-(
- MoleculeType& molI,
- MoleculeType& molJ
-) const
-{
- const pairPotentialList& pairPot = pot_.pairPotentials();
-
- const pairPotential& electrostatic = pairPot.electrostatic();
-
- label idI = molI.id();
-
- label idJ = molJ.id();
-
- const typename MoleculeType::constantProperties& constPropI
- (
- constProps(idI)
- );
-
- const typename MoleculeType::constantProperties& constPropJ
- (
- constProps(idJ)
- );
-
- forAll(constPropI.pairPotSites(), pI)
- {
- label sI = constPropI.pairPotSites()[pI];
-
- label idsI = constPropI.sites()[sI].siteId();
-
- forAll(constPropJ.pairPotSites(), pJ)
- {
- label sJ = constPropJ.pairPotSites()[pJ];
-
- label idsJ = constPropJ.sites()[sJ].siteId();
-
- vector rsIsJ =
- molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
-
- scalar rsIsJMagSq = magSqr(rsIsJ);
-
- if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
- {
- scalar rsIsJMag = mag(rsIsJ);
-
- // Guard against pairPot.energy being evaluated
- // if rIJMag < SMALL. A floating point exception will
- // happen otherwise.
-
- if (rsIsJMag < SMALL)
- {
- WarningIn
- (
- "MoleculeCloud::"
- "removeHighEnergyOverlaps()"
- )
- << "Molecule site pair closer than "
- << SMALL
- << ": mag separation = " << rsIsJMag
- << ". These may have been placed on top of each"
- << " other by a rounding error in mdInitialise in"
- << " parallel or a block filled with moleculess"
- << " twice. Removing one of the molecules."
- << endl;
-
- return true;
- }
-
- // Guard against pairPot.energy being evaluated if rIJMag <
- // rMin. A tabulation lookup error will occur otherwise.
-
- if (rsIsJMag < pairPot.rMin(idsI, idsJ))
- {
- return true;
- }
-
- if
- (
- mag(pairPot.energy(idsI, idsJ, rsIsJMag))
- > pot_.potentialEnergyLimit()
- )
- {
- return true;
- };
- }
- }
- }
-
- forAll(constPropI.electrostaticSites(), pI)
- {
- label sI = constPropI.electrostaticSites()[pI];
-
- forAll(constPropJ.electrostaticSites(), pJ)
- {
- label sJ = constPropJ.electrostaticSites()[pJ];
-
- vector rsIsJ =
- molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
-
- scalar rsIsJMagSq = magSqr(rsIsJ);
-
- if (pairPot.rCutMaxSqr(rsIsJMagSq))
- {
- scalar rsIsJMag = mag(rsIsJ);
-
- // Guard against pairPot.energy being evaluated
- // if rIJMag < SMALL. A floating point exception will
- // happen otherwise.
-
- if (rsIsJMag < SMALL)
- {
- WarningIn
- (
- "MoleculeCloud::"
- "removeHighEnergyOverlaps()"
- )
- << "Molecule site pair closer than "
- << SMALL
- << ": mag separation = " << rsIsJMag
- << ". These may have been placed on top of each"
- << " other by a rounding error in mdInitialise in"
- << " parallel or a block filled with molecules"
- << " twice. Removing one of the molecules."
- << endl;
-
- return true;
- }
-
- if (rsIsJMag < electrostatic.rMin())
- {
- return true;
- }
-
- scalar chargeI = constPropI.sites()[sI].siteCharge();
-
- scalar chargeJ = constPropJ.sites()[sJ].siteCharge();
-
- if
- (
- mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
- > pot_.potentialEnergyLimit()
- )
- {
- return true;
- };
- }
- }
- }
-
- return false;
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-template
-inline const Foam::polyMesh& Foam::MoleculeCloud::mesh() const
-{
- return mesh_;
-}
-
-
-template
-inline const Foam::potential& Foam::MoleculeCloud::pot() const
-{
- return pot_;
-}
-
-
-template
-inline const Foam::List >&
-Foam::MoleculeCloud::cellOccupancy() const
-{
- return cellOccupancy_;
-}
-
-
-template
-inline const Foam::InteractionLists&
-Foam::MoleculeCloud::il() const
-{
- return il_;
-}
-
-
-template
-inline const Foam::List
-Foam::MoleculeCloud::constProps() const
-{
- return constPropList_;
-}
-
-
-template
-inline const typename MoleculeType::constantProperties&
-Foam::MoleculeCloud::constProps(label id) const
-{
- return constPropList_[id];
-}
-
-
-template
-inline Foam::Random& Foam::MoleculeCloud::rndGen()
-{
- return rndGen_;
-}
-
-
-template
-inline Foam::vector
-Foam::MoleculeCloud::equipartitionLinearVelocity
-(
- scalar temperature,
- scalar mass
-)
-{
- return sqrt(physicoChemical::k.value()*temperature/mass)*vector
- (
- rndGen_.GaussNormal(),
- rndGen_.GaussNormal(),
- rndGen_.GaussNormal()
- );
-}
-
-
-template
-inline Foam::vector
-Foam::MoleculeCloud::equipartitionAngularMomentum
-(
- scalar temperature,
- const typename MoleculeType::constantProperties& cP
-)
-{
- scalar sqrtKbT = sqrt(physicoChemical::k.value()*temperature);
-
- if (cP.linearMolecule())
- {
- return sqrtKbT*vector
- (
- 0.0,
- sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
- sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
- );
- }
- else
- {
- return sqrtKbT*vector
- (
- sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal(),
- sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
- sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
- );
- }
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.C
deleted file mode 100644
index 4fda04df7d..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.C
+++ /dev/null
@@ -1,48 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see monoatomicCloud;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.H b/src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.H
deleted file mode 100644
index bb902fe3a4..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.H
+++ /dev/null
@@ -1,52 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see polyatomicCloud;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C
deleted file mode 100644
index f0bad336ee..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C
+++ /dev/null
@@ -1,573 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see
- (
- stateController::New(time_, cloud, controllersIDict)
- );
-
- sCNames_[sC] = stateControllers_[sC]->type();
-
- sCIds_[sC] = sC;
- }
- }
-
- //- flux controllers
-
- if (!fluxControllers_.empty())
- {
- forAll(fluxControllers_, fC)
- {
- const entry& controllersI = fluxControllersList_[fC];
-
- const dictionary& controllersIDict = controllersI.dict();
-
- fluxControllers_[fC] = autoPtr
- (
- fluxController::New(time_, cloud, controllersIDict)
- );
-
- fCNames_[fC] = fluxControllers_[fC]->type();
- fCIds_[fC] = fC;
- }
- }
-
- // creating directories for state controllers
- if (!nStateControllers_.empty())
- {
- // case/controllers
- fileName controllersPath(time_.path()/"controllers");
-
- if (!isDir(controllersPath))
- {
- mkDir(controllersPath);
- }
-
- // case/controllers/
- fileName controllersPath(controllersPath/cloud.name());
-
- if (!isDir(controllersPath))
- {
- mkDir(controllersPath);
- }
-
- // case/controllers//stateControllers
- fileName stateControllersPath(controllersPath/"stateControllers");
-
- if (!isDir(stateControllersPath))
- {
- mkDir(stateControllersPath);
- }
-
- forAll(stateControllers_, sC)
- {
- if (stateControllers_[sC]->writeInCase())
- {
- // case/controllers//
- // stateControllers/
- fileName stateControllerPath(stateControllersPath/sCNames_[sC]);
-
- if (!isDir(stateControllerPath))
- {
- mkDir(stateControllerPath);
- }
-
- const word& regionName = stateControllers_[sC]->regionName();
-
- // case/controllers//
- // stateControllers//
- fileName zonePath(stateControllerPath/regionName);
-
- if (!isDir(zonePath))
- {
- mkDir(zonePath);
- }
-
- sCFixedPathNames_[sC] = zonePath;
- }
- }
- }
-
- // creating directories for flux controllers
- if (nFluxControllers_ > 0)
- {
- // case/controllers
- fileName controllersPath(time_.path()/"controllers");
-
- if ( !isDir(controllersPath) )
- {
- mkDir(controllersPath);
- }
-
- // case/controllers/
- fileName controllersPath(time_.path()/cloud.name());
-
- if ( !isDir(controllersPath) )
- {
- mkDir(controllersPath);
- }
-
- // case/controllers//fluxControllers
- fileName fluxControllersPath(controllersPath/"fluxControllers");
-
- if (!isDir(fluxControllersPath))
- {
- mkDir(fluxControllersPath);
- }
-
- forAll(fluxControllers_, fC)
- {
- if (fluxControllers_[fC]->writeInCase())
- {
- // case/controllers//
- // fluxControllers/
- fileName fluxControllerPath(fluxControllersPath/fCNames_[fC]);
-
- if (!isDir(fluxControllerPath))
- {
- mkDir(fluxControllerPath);
- }
-
- const word& regionName = fluxControllers_[fC]->regionName();
-
- // case/controllers//
- // fluxControllers//
- fileName zonePath(fluxControllerPath/regionName);
-
- if (!isDir(zonePath))
- {
- mkDir(zonePath);
- }
-
- fCFixedPathNames_[fC] = zonePath;
- }
- }
- }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-controllers::~controllers()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-void controllers::initialConfig()
-{
- forAll(stateControllers_, sC)
- {
- stateControllers_[sC]->initialConfiguration();
- }
-
- forAll(fluxControllers_, fC)
- {
- fluxControllers_[fC]->initialConfiguration();
- }
-}
-
-
-void controllers::updateTimeInfo()
-{
- forAll(stateControllers_, sC)
- {
- stateControllers_[sC]->updateTime();
- }
-
- forAll(fluxControllers_, fC)
- {
- fluxControllers_[fC]->updateTime();
- }
-}
-
-
-void controllers::controlState()
-{
- forAll(stateControllers_, sC)
- {
- stateControllers_[sC]->controlMols();
- }
-}
-
-
-void controllers::controlVelocitiesI()
-{
- forAll(stateControllers_, sC)
- {
- stateControllers_[sC]->controlMolsBeg();
- }
-}
-
-void controllers::controlVelocitiesII()
-{
- forAll(stateControllers_, sC)
- {
- stateControllers_[sC]->controlMolsEnd();
- }
-}
-
-
-void controllers::controlPriorToForces()
-{
- forAll(stateControllers_, sC)
- {
- stateControllers_[sC]->controlBeforeForces();
- }
-}
-
-
-void controllers::calculateStateProps()
-{
- forAll(stateControllers_, sC)
- {
- stateControllers_[sC]->calculateProperties();
- }
-
- forAll(fluxControllers_, fC)
- {
- fluxControllers_[fC]->calculateProperties();
- }
-}
-
-
-void controllers::outputStateResults()
-{
- const Time& runTime = time_;
-
- if (runTime.outputTime())
- {
- // creating a set of directories in the current time directory
- {
- List timePathNames(sCFixedPathNames_.size());
-
- if (nStateControllers_ > 0)
- {
- if (Pstream::master())
- {
- // case//uniform
- fileName uniformTimePath
- (
- runTime.path()/runTime.timeName()/"uniform"
- );
-
- if (!isDir(uniformTimePath))
- {
- mkDir(uniformTimePath);
- }
-
- if (!stateControllers_.empty())
- {
- // case//uniform/controllers
- fileName controllersTimePath
- (
- uniformTimePath/"controllers"
- );
-
- if (!isDir(controllersTimePath))
- {
- mkDir(controllersTimePath);
- }
-
- // case//uniform/controllers/
- fileName cloudTimePath
- (
- controllersTimePath/cloud.name()
- );
-
- if (!isDir(cloudTimePath))
- {
- mkDir(cloudTimePath);
- }
-
- // case//uniform/controllers//
- fileName stateControllersTimePath
- (
- cloudTimePath/"stateControllers"
- );
-
- if (!isDir(stateControllersTimePath))
- {
- mkDir(stateControllersTimePath);
- }
-
- forAll(stateControllers_, sC)
- {
- if (stateControllers_[sC]->writeInTimeDir())
- {
- // case//uniform/controllers/
- // /
- fileName sCTimePath
- (
- stateControllersTimePath/sCNames_[sC]
- );
-
- if (!isDir(sCTimePath))
- {
- mkDir(sCTimePath);
- }
-
- // Creating directory for different zones but
- // of the same model
- const word& regionName =
- stateControllers_[sC]->regionName();
-
- // case//uniform/controllers/
- // //
- //
- fileName zoneTimePath(sCTimePath/regionName);
-
- if (!isDir(zoneTimePath))
- {
- mkDir(zoneTimePath);
- }
-
- timePathNames[sC] = zoneTimePath;
- }
- }
- }
- }
- }
-
- // write out data
- forAll(stateControllers_, sC)
- {
- stateControllers_[sC]->output
- (
- sCFixedPathNames_[sC],
- timePathNames[sC]
- );
- }
- }
-
- {
- List timePathNames(fCFixedPathNames_.size());
-
- if (nFluxControllers_ > 0)
- {
- if (Pstream::master())
- {
- // case//uniform
- fileName uniformTimePath
- (
- runTime.path()/runTime.timeName()/"uniform"
- );
-
- if (!isDir(uniformTimePath))
- {
- mkDir(uniformTimePath);
- }
-
- if (!fluxControllers_.empty())
- {
- // case//uniform/controllers
- fileName controllersTimePath
- (
- uniformTimePath/"controllers"
- );
-
- if (!isDir(controllersTimePath))
- {
- mkDir(controllersTimePath);
- }
-
- // case//uniform/controllers/
- fileName cloudTimePath
- (
- controllersTimePath/cloud.name()
- );
-
- if (!isDir(cloudTimePath))
- {
- mkDir(cloudTimePath);
- }
-
- // case//uniform/fluxControllers
- fileName controllersTimePath
- (
- cloudTimePath/"fluxControllers"
- );
-
- if (!isDir(controllersTimePath))
- {
- mkDir(controllersTimePath);
- }
-
- forAll(fluxControllers_, fC)
- {
- if (stateControllers_[fC]->writeInTimeDir())
- {
- // case//uniform/controllers/
- // /
- fileName fCTimePath
- (
- controllersTimePath/fCNames_[fC]
- );
-
- if (!isDir(fCTimePath))
- {
- mkDir(fCTimePath);
- }
-
- const word& regionName =
- fluxControllers_[fC]->regionName();
-
- // case//uniform/controllers/
- // //
- //
- fileName zoneTimePath(fCTimePath/regionName);
-
- if (!isDir(zoneTimePath))
- {
- mkDir(zoneTimePath);
- }
-
- timePathNames[fC] = zoneTimePath;
- }
- }
- }
- }
- }
-
- // write out data
- forAll(fluxControllers_, fC)
- {
- fluxControllers_[fC]->output
- (
- fCFixedPathNames_[fC],
- timePathNames[fC]
- );
- }
- }
-
- // Re-read dictionaries for modified properties (run-time selection)
- {
- stateControllersList_.clear();
-
- stateControllersList_ = controllersDict_.lookup("stateControllers");
-
- forAll(stateControllers_, sC)
- {
- const entry& controllersI = stateControllersList_[sC];
- const dictionary& controllersIDict = controllersI.dict();
-
- stateControllers_[sC]->updateProperties(controllersIDict);
- }
- }
-
- {
- fluxControllersList_.clear();
-
- fluxControllersList_ = controllersDict_.lookup("fluxControllers");
-
- forAll(fluxControllers_, fC)
- {
- const entry& controllersI = fluxControllersList_[fC];
- const dictionary& controllersIDict = controllersI.dict();
-
- fluxControllers_[fC]->updateProperties(controllersIDict);
- }
- }
- }
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H
deleted file mode 100644
index 22e9ca7cca..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H
+++ /dev/null
@@ -1,163 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see stateControllersList_;
- List sCNames_;
- List sCIds_;
- List sCFixedPathNames_;
- List< autoPtr > stateControllers_;
-
- //- flux controllers
- PtrList fluxControllersList_;
- List fCNames_;
- List fCIds_;
- List fCFixedPathNames_;
- List< autoPtr > fluxControllers_;
-
-
-public:
-
- // Constructors
-
- //- Null Constructor
- controllers
- (
- const polyMesh& mesh
- );
-
- //- Constructor for with cloud
- controllers
- (
- const polyMesh& mesh,
- polyatomicCloud& cloud
- );
-
- //- Destructor
- ~controllers();
-
- // Member Functions
-
- //- Initial configuration call this function after the polyatomicCloud
- // is completely initialised
- void initialConfig();
-
- //- this function is to be called at the beginning of the MD time-step.
- // since we have placed a non-referenced time-data class in the
- // state-controller class.
- void updateTimeInfo();
-
- //- control molecular state -- call this after the intermolecular force
- // calulation
- void controlState();
-
- //-
- void controlVelocitiesI();
-
- //-
- void controlVelocitiesII();
-
- //-
- void controlPriorToForces();
-
- //- calculate properties -- call this at the end of the MD time-step.
- void calculateStateProps();
-
- //- output -- call this function at the end of the MD time-step
- void outputStateResults();
-
- // Access
-
- //-
- inline List< autoPtr >& stateControllers();
-
- //-
- inline const List< autoPtr >&
- stateControllers() const;
-
- //-
- inline List< autoPtr >& fluxControllers();
-
- //-
- inline const List< autoPtr >&
- fluxControllers() const;
-
- //-
- inline const List& stateControllersNames() const;
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "controllersI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H
deleted file mode 100644
index 66c10f0784..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H
+++ /dev/null
@@ -1,64 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see >&
-Foam::controllers::stateControllers()
-{
- return stateControllers_;
-}
-
-
-const Foam::List >&
-Foam::controllers::stateControllers() const
-{
- return stateControllers_;
-}
-
-
-Foam::List >&
-Foam::controllers::fluxControllers()
-{
- return fluxControllers_;
-}
-
-
-const Foam::List< autoPtr >&
-Foam::controllers::fluxControllers() const
-{
- return fluxControllers_;
-}
-
-
-const Foam::List& Foam::controllers::stateControllersNames() const
-{
- return sCNames_;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C b/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C
deleted file mode 100644
index 711b661758..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C
+++ /dev/null
@@ -1,524 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see (cloud.mesh())),
- cloud_(cloud),
- rndGen_(clock::getTime()),
- controllerDict_(dict.subDict("controllerProperties")),
- timeDict_(controllerDict_.subDict("timeProperties")),
- time_(t, timeDict_),
- regionName_(controllerDict_.lookup("zoneName")),
- regionId_(-1),
- zoneSurfaceArea_(0.0),
- internalFaces_(),
- processorFaces_(),
- control_(true),
- readStateFromFile_(true),
- singleValueController_(false),
- density_(0.0),
- velocity_(vector::zero),
- temperature_(0.0),
- pressure_(0.0),
- strainRate_(tensor::zero),
- tempGradient_(vector::zero),
- fieldController_(false),
- densities_(),
- velocities_(),
- temperatures_(),
- pressures_(),
- writeInTimeDir_(true),
- writeInCase_(true)
-{
- const faceZoneMesh& faceZones = mesh_.faceZones();
- regionId_ = faceZones.findZoneID(regionName_);
-
- if (regionId_ == -1)
- {
- FatalErrorIn("waterFluxController::waterFluxController()")
- << "Cannot find region (faceZone): " << regionName_ << nl << "in: "
- << t.time().system()/"controllersDict"
- << exit(FatalError);
- }
-
- control_ = Switch(controllerDict_.lookup("controlSwitch"));
- readStateFromFile_ = Switch(controllerDict_.lookup("readStateFromFile"));
-
- setFacesInfo();
-}
-
-
-// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
-
-autoPtr waterFluxController::New
-(
- Time& t,
- polyatomicCloud& cloud,
- const dictionary& dict
-)
-{
- word waterFluxControllerName
- (
- dict.lookup("fluxControllerModel")
- );
-
- Info<< "Selecting fluxController "
- << waterFluxControllerName << endl;
-
- dictionaryConstructorTable::iterator cstrIter =
- dictionaryConstructorTablePtr_->find(waterFluxControllerName);
-
- if (cstrIter == dictionaryConstructorTablePtr_->end())
- {
- FatalError
- << "waterFluxController::New(const dictionary&) : " << endl
- << " unknown waterFluxController type "
- << waterFluxControllerName
- << ", constructor not in hash table" << endl << endl
- << " Valid injector types are :" << endl;
- Info<< dictionaryConstructorTablePtr_->toc() << abort(FatalError);
- }
-
- return autoPtr
- (
- cstrIter()(t, cloud, dict)
- );
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-waterFluxController::~waterFluxController()
-{}
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-// void waterFluxController::updateTime()
-// {
-// time_++;
-//
-// const scalar& t = time_.time().timeOutputValue();
-//
-// if ((t - initialTime_) < timePeriod_)
-// {
-// time_.controlTimeInterval().endTime() = false;
-// // control_ = false;
-// }
-// else
-// {
-// // control_ = true;
-// }
-// }
-
-
-void waterFluxController::setFacesInfo()
-{
- const labelList& faces = controlZone();
-
- if (Pstream::parRun())
- {
- DynamicList processorFaces(0);
-
- forAll(mesh_.boundaryMesh(), patchI)
- {
- const polyPatch& patch = mesh_.boundaryMesh()[patchI];
-
- if (isA(patch))
- {
- for (label p = 0; p < patch.size(); p++)
- {
- label patchFaceI = p + patch.start();
- label faceId = findIndex (faces, patchFaceI);
-
- if (faceId != -1)
- {
- processorFaces.append(patchFaceI);
- }
- }
- }
- }
-
- processorFaces.shrink();
-
- processorFaces_.setSize(processorFaces.size(), -1);
-
- forAll(processorFaces, f)
- {
- processorFaces_[f] = processorFaces[f];
- }
-
- label nInternalFaces = faces.size() - processorFaces.size();
- internalFaces_.setSize(nInternalFaces, -1);
-
- label counter = 0;
-
- forAll(faces, f)
- {
- const label& faceI = faces[f];
-
- if (findIndex(processorFaces, faceI) == -1)
- {
- internalFaces_[counter] = faceI;
- counter++;
- }
- }
-
-// Pout << "processorFaces: " << processorFaces_ << endl;
-// Pout << "internalFaces: " << internalFaces_ << endl;
-
- forAll(internalFaces_, f)
- {
- const label& faceI = internalFaces_[f];
- zoneSurfaceArea_ += mag(mesh_.faceAreas()[faceI]);
- }
-
- // faces on a zone located on a processor cut belong to both processors
- // (hence the 0.5)
-
- forAll(processorFaces_, f)
- {
- const label& faceI = processorFaces_[f];
- zoneSurfaceArea_ += 0.5*mag(mesh_.faceAreas()[faceI]);
- }
-
-
- if (Pstream::parRun())
- {
- for (int p = 0; p < Pstream::nProcs(); p++)
- {
- if (p != Pstream::myProcNo())
- {
- const int proc = p;
- {
- OPstream toNeighbour(Pstream::blocking, proc);
- toNeighbour << zoneSurfaceArea_;
- }
- }
- }
-
- //- receiving
- for (int p = 0; p < Pstream::nProcs(); p++)
- {
- if (p != Pstream::myProcNo())
- {
- scalar zoneSurfaceAreaProc;
-
- const int proc = p;
- {
- IPstream fromNeighbour(Pstream::blocking, proc);
- fromNeighbour >> zoneSurfaceAreaProc;
- }
-
- zoneSurfaceArea_ += zoneSurfaceAreaProc;
- }
- }
- }
- }
- else
- {
- forAll(faces, f)
- {
- const label& faceI = faces[f];
-
- zoneSurfaceArea_ += mag(mesh_.faceAreas()[faceI]);
- }
- }
-}
-
-
-
-void waterFluxController::updateTime()
-{
- time_++;
-
-// const scalar& t = time_.time().timeOutputValue();
-//
-// if ((t - initialTime_) < timePeriod_)
-// {
-// time_.controlTimeInterval().endTime() = false;
-// // control_ = false;
-// }
-// else
-// {
-// // control_ = true;
-// }
-}
-
-void waterFluxController::updateFluxControllerProperties
-(
- const dictionary& newDict
-)
-{
- controllerDict_ = newDict.subDict("controllerProperties");
-
- //- you can reset the controlling zone from here. This essentially
- // means that the coupling zone can infact move arbitrarily. To make
- // this happen we probably need to devise a technique for automatically
- // changing the cellZone else where, and then calling this function to
- // reset the controlling zone in which the controller operates in.
-
- if (controllerDict_.found("controlSwitch"))
- {
- control_ = Switch(controllerDict_.lookup("controlSwitch"));
- }
-
- if (controllerDict_.found("readStateFromFile"))
- {
- readStateFromFile_ = Switch
- (
- controllerDict_.lookup("readStateFromFile")
- );
- }
-}
-
-const labelList& waterFluxController::controlZone() const
-{
- return mesh_.faceZones()[regionId_];
-}
-
-label waterFluxController::isFaceOnControlZone(const label& faceI)
-{
- const label f = findIndex(controlZone(), faceI);
-
- return f;
-}
-
-const word& waterFluxController::regionName() const
-{
- return regionName_;
-}
-
-const scalar& waterFluxController::density() const
-{
- return density_;
-}
-
-scalar& waterFluxController::density()
-{
- return density_;
-}
-
-const vector& waterFluxController::velocity() const
-{
- return velocity_;
-}
-
-vector& waterFluxController::velocity()
-{
- return velocity_;
-}
-
-const scalar& waterFluxController::temperature() const
-{
- return temperature_;
-}
-
-scalar& waterFluxController::temperature()
-{
- return temperature_;
-}
-
-const scalar& waterFluxController::pressure() const
-{
- return pressure_;
-}
-
-scalar& waterFluxController::pressure()
-{
- return pressure_;
-}
-
-const tensor& waterFluxController::strainRate() const
-{
- return strainRate_;
-}
-
-tensor& waterFluxController::strainRate()
-{
- return strainRate_;
-}
-
-const vector& waterFluxController::tempGradient() const
-{
- return tempGradient_;
-}
-
-vector& waterFluxController::tempGradient()
-{
- return tempGradient_;
-}
-
-
-const scalarField& waterFluxController::densityField() const
-{
- return densities_;
-}
-
-scalarField& waterFluxController::densityField()
-{
- return densities_;
-}
-
-const vectorField& waterFluxController::velocityField() const
-{
- return velocities_;
-}
-vectorField& waterFluxController::velocityField()
-{
- return velocities_;
-}
-
-const scalarField& waterFluxController::temperatureField() const
-{
- return temperatures_;
-}
-
-scalarField& waterFluxController::temperatureField()
-{
- return temperatures_;
-}
-
-const scalarField& waterFluxController::pressureField() const
-{
- return pressures_;
-}
-
-scalarField& waterFluxController::pressureField()
-{
- return pressures_;
-}
-
-
-const bool& waterFluxController::singleValueController() const
-{
- return singleValueController_;
-}
-
-bool& waterFluxController::singleValueController()
-{
- return singleValueController_;
-}
-
-const bool& waterFluxController::fieldController() const
-{
- return fieldController_;
-}
-
-bool& waterFluxController::fieldController()
-{
- return fieldController_;
-}
-
-
-const bool& waterFluxController::writeInTimeDir() const
-{
- return writeInTimeDir_;
-}
-
-const bool& waterFluxController::writeInCase() const
-{
- return writeInCase_;
-}
-
-
-// const scalar waterFluxController::avReqDensity() const
-// {
-// scalar totalDensity = 0.0;
-//
-// forAll(densities_, c)
-// {
-// totalDensity += densities_[c];
-// }
-//
-// if (cells_.size() > 0)
-// {
-// totalDensity /= scalar(cells_.size());
-// }
-//
-// return totalDensity;
-// }
-//
-// const vector waterFluxController::avReqVelocity() const
-// {
-// vector totalVel = vector::zero;
-//
-// forAll(velocities_, c)
-// {
-// totalVel += velocities_[c];
-// }
-//
-// if (cells_.size() > 0)
-// {
-// totalVel /= scalar(cells_.size());
-// }
-//
-// return totalVel;
-// }
-//
-// const scalar waterFluxController::avReqTemperature() const
-// {
-// scalar totalTemp = 0.0;
-//
-// forAll(densities_, c)
-// {
-// totalTemp += temperatures_[c];
-// }
-//
-// if (cells_.size() > 0)
-// {
-// totalTemp /= scalar(cells_.size());
-// }
-//
-// return totalTemp;
-// }
-
-
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H
deleted file mode 100644
index 9b8a375277..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H
+++ /dev/null
@@ -1,272 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see New
- (
- Time& t,
- polyatomicCloud& cloud,
- const dictionary& dict
- );
-
-
- // Destructor
-
- virtual ~waterFluxController();
-
-
- // Member Functions
-
- void updateTime();
-
- //- create an initial configuration
- virtual void initialConfiguration() = 0;
-
-
- //- calculate any required properties
- virtual void calculateProperties() = 0;
-
-
- //- control the polyatomic from the tracking function
- virtual void controlMol
- (
- polyatomic& mol,
- polyatomic::trackData& td
- ) = 0;
-
- //- output data
- virtual void output
- (
- const fileName& fixedPathName,
- const fileName& timePath
- ) = 0;
-
-
- //- E. update properties from a modified dictionary
- virtual void updateProperties(const dictionary&) = 0;
-
-
- void updateFluxControllerProperties(const dictionary&);
-
- // Access
-
- //- return the control zone cells
- const labelList& controlZone() const;
-
-
- label isFaceOnControlZone(const label& faceI);
-
- //- return the control zone name
- const word& regionName() const;
-
- //- return the targeted values
- const scalar& density() const;
- scalar& density();
-
- const vector& velocity() const;
- vector& velocity();
-
- const scalar& temperature() const;
- scalar& temperature();
-
- const scalar& pressure() const;
- scalar& pressure();
-
- const tensor& strainRate() const;
- tensor& strainRate();
-
- const vector& tempGradient() const;
- vector& tempGradient();
-
- //- return the targeted fields
- const scalarField& densityField() const;
- scalarField& densityField();
-
- const vectorField& velocityField() const;
- vectorField& velocityField();
-
- const scalarField& temperatureField() const;
- scalarField& temperatureField();
-
- const scalarField& pressureField() const;
- scalarField& pressureField();
-
-
- const bool& singleValueController() const;
- bool& singleValueController();
-
- const bool& fieldController() const;
- bool& fieldController();
-
- const bool& writeInTimeDir() const;
- const bool& writeInCase() const;
-
-// const scalar avReqDensity() const;
-// const vector avReqVelocity() const;
-// const scalar avReqTemperature() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C b/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C
deleted file mode 100644
index b98aa126a3..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C
+++ /dev/null
@@ -1,490 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see (cloud.mesh())),
- cloud_(cloud),
- rndGen_(clock::getTime()),
- controllerDict_(dict.subDict("controllerProperties")),
- timeDict_(controllerDict_.subDict("timeProperties")),
- time_(mesh_.time(), timeDict_),
- timePeriod_(readScalar(timeDict_.lookup("initialTimePeriod"))), //temp
- initialTime_(time_.time().startTime().value()),
- regionName_(controllerDict_.lookup("zoneName")),
- regionId_(-1),
- control_(true),
- readStateFromFile_(true),
- singleValueController_(false),
- density_(0.0),
- velocity_(vector::zero),
- temperature_(0.0),
- pressure_(0.0),
- strainRate_(tensor::zero),
- tempGradient_(vector::zero),
- fieldController_(false),
- densities_(),
- velocities_(),
- temperatures_(),
- pressures_(),
- writeInTimeDir_(true),
- writeInCase_(true)
-{
- const cellZoneMesh& cellZones = mesh_.cellZones();
-
- regionId_ = cellZones.findZoneID(regionName_);
-
- if (regionId_ == -1)
- {
- FatalErrorIn("stateController::stateController()")
- << "Cannot find region: " << regionName_ << nl << "in: "
- << time_.time().system()/"controllersDict"
- << exit(FatalError);
- }
-
- control_ = Switch(controllerDict_.lookup("controlSwitch"));
-
- readStateFromFile_ = Switch(controllerDict_.lookup("readStateFromFile"));
-
- const scalar& avTimeInterval = time_.averageTimeInterval().deltaT();
-
- if ((timePeriod_ < avTimeInterval) && (timePeriod_ > 0.0))
- {
- timePeriod_ = avTimeInterval;
- }
-}
-
-
-// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
-
-Foam::autoPtr Foam::stateController::New
-(
- polyatomicCloud& cloud,
- const dictionary& dict
-)
-{
- word stateControllerName
- (
- dict.lookup("stateControllerModel")
- );
-
- Info<< "Selecting stateController "
- << stateControllerName << endl;
-
- dictionaryConstructorTable::iterator cstrIter =
- dictionaryConstructorTablePtr_->find(stateControllerName);
-
- if (cstrIter == dictionaryConstructorTablePtr_->end())
- {
- FatalError
- << "stateController::New(const dictionary&) : " << endl
- << " unknown stateController type "
- << stateControllerName
- << ", constructor not in hash table" << endl << endl
- << " Valid types are :" << endl;
- Info<< dictionaryConstructorTablePtr_->toc() << abort(FatalError);
- }
-
- return autoPtr
- (
- cstrIter()(cloud, dict)
- );
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-Foam::stateController::~stateController()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-void Foam::stateController::updateTime()
-{
- time_++;
-
- const scalar& t = time_.time().timeOutputValue();
-
- if ((t - initialTime_) < timePeriod_)
- {
- time_.controlTimeInterval().endTime() = false;
- }
-}
-
-
-void Foam::stateController::updateStateControllerProperties
-(
- const dictionary& newDict
-)
-{
- controllerDict_ = newDict.subDict("controllerProperties");
-
- if (controllerDict_.found("controlSwitch"))
- {
- control_ = Switch(controllerDict_.lookup("controlSwitch"));
- }
-
- if (controllerDict_.found("readStateFromFile"))
- {
- readStateFromFile_ = Switch
- (
- controllerDict_.lookup("readStateFromFile")
- );
- }
-
- timeDict_ = controllerDict_.subDict("timeProperties");
-
- if (timeDict_.found("resetAtOutput"))
- {
- time_.resetFieldsAtOutput() = Switch(timeDict_.lookup("resetAtOutput"));
- }
-}
-
-
-const Foam::labelList& Foam::stateController::controlZone() const
-{
- return mesh_.cellZones()[regionId_];
-}
-
-const Foam::word& Foam::stateController::regionName() const
-{
- return regionName_;
-}
-
-
-Foam::scalar Foam::stateController::density() const
-{
- return density_;
-}
-
-
-Foam::scalar& Foam::stateController::density()
-{
- return density_;
-}
-
-
-const Foam::vector& Foam::stateController::velocity() const
-{
- return velocity_;
-}
-
-
-Foam::vector& Foam::stateController::velocity()
-{
- return velocity_;
-}
-
-
-Foam::scalar Foam::stateController::temperature() const
-{
- return temperature_;
-}
-
-
-Foam::scalar& Foam::stateController::temperature()
-{
- return temperature_;
-}
-
-const Foam::scalar& Foam::stateController::pressure() const
-{
- return pressure_;
-}
-
-
-Foam::scalar& Foam::stateController::pressure()
-{
- return pressure_;
-}
-
-
-const Foam::tensor& Foam::stateController::strainRate() const
-{
- return strainRate_;
-}
-
-
-Foam::tensor& Foam::stateController::strainRate()
-{
- return strainRate_;
-}
-
-
-const Foam::vector& Foam::stateController::tempGradient() const
-{
- return tempGradient_;
-}
-
-
-Foam::vector& Foam::stateController::tempGradient()
-{
- return tempGradient_;
-}
-
-
-const Foam::scalarField& Foam::stateController::densityField() const
-{
- return densities_;
-}
-
-Foam::scalarField& Foam::stateController::densityField()
-{
- return densities_;
-}
-
-
-const Foam::vectorField& Foam::stateController::velocityField() const
-{
- return velocities_;
-}
-
-
-Foam::vectorField& Foam::stateController::velocityField()
-{
- return velocities_;
-}
-
-
-const Foam::scalarField& Foam::stateController::temperatureField() const
-{
- return temperatures_;
-}
-
-
-Foam::scalarField& Foam::stateController::temperatureField()
-{
- return temperatures_;
-}
-
-
-const Foam::scalarField& Foam::stateController::pressureField() const
-{
- return pressures_;
-}
-
-
-Foam::scalarField& Foam::stateController::pressureField()
-{
- return pressures_;
-}
-
-
-bool Foam::stateController::singleValueController() const
-{
- return singleValueController_;
-}
-
-
-bool& Foam::stateController::singleValueController()
-{
- return singleValueController_;
-}
-
-
-bool Foam::stateController::fieldController() const
-{
- return fieldController_;
-}
-
-
-bool& Foam::stateController::fieldController()
-{
- return fieldController_;
-}
-
-
-bool Foam::stateController::writeInTimeDir() const
-{
- return writeInTimeDir_;
-}
-
-
-bool Foam::stateController::writeInCase() const
-{
- return writeInCase_;
-}
-
-
-Foam::scalar Foam::stateController::avReqDensity()
-{
- scalar totalDensity = 0.0;
-
- if (singleValueController_)
- {
- totalDensity = density_;
- }
- else if (fieldController_)
- {
- label controlCells = controlZone().size();
-
- forAll(densities_, c)
- {
- totalDensity += densities_[c];
- }
-
- if (Pstream::parRun())
- {
- reduce(totalDensity, sumOp());
-
- reduce(controlCells, sumOp());
- }
-
- if (controlCells > 0)
- {
- totalDensity /= scalar(controlCells);
- }
- }
-
- return totalDensity;
-}
-
-
-Foam::vector Foam::stateController::avReqVelocity()
-{
- vector totalVel = vector::zero;
-
- if (singleValueController_)
- {
- totalVel = velocity_;
- }
- else if (fieldController_)
- {
- label controlCells = controlZone().size();
-
- forAll(velocities_, c)
- {
- totalVel += velocities_[c];
- }
-
- if (Pstream::parRun())
- {
- reduce(totalVel, sumOp());
-
- reduce(controlCells, sumOp());
- }
-
- if (controlCells > 0)
- {
- totalVel /= scalar(controlCells);
- }
- }
-
- return totalVel;
-}
-
-
-Foam::scalar Foam::stateController::avReqTemperature()
-{
- scalar totalTemp = 0.0;
-
- if (singleValueController_)
- {
- totalTemp = temperature_;
- }
- else if (fieldController_)
- {
- label controlCells = controlZone().size();
-
- forAll(temperatures_, c)
- {
- totalTemp += temperatures_[c];
- }
-
- if (Pstream::parRun())
- {
- reduce(totalTemp, sumOp());
-
- reduce(controlCells, sumOp());
- }
-
- if (controlCells > 0)
- {
- totalTemp /= scalar(controlCells);
- }
- }
-
- return totalTemp;
-}
-
-
-Foam::scalar Foam::stateController::avReqPressure()
-{
- scalar totalPressure = 0.0;
-
- if (singleValueController_)
- {
- totalPressure = pressure_;
- }
- else if (fieldController_)
- {
- label controlCells = controlZone().size();
-
- forAll(pressures_, c)
- {
- totalPressure += pressures_[c];
- }
-
- if (Pstream::parRun())
- {
- reduce(totalPressure, sumOp());
-
- reduce(controlCells, sumOp());
- }
-
- if (controlCells > 0)
- {
- totalPressure /= scalar(controlCells);
- }
- }
-
- return totalPressure;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H b/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H
deleted file mode 100644
index 3e142c3962..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H
+++ /dev/null
@@ -1,270 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see New
- (
- polyatomicCloud& cloud,
- const dictionary& dict
- );
-
-
- // Destructor
-
- virtual ~stateController();
-
-
- // Member Functions
-
-
- void updateTime();
-
- //- create an initial configuration
- virtual void initialConfiguration() = 0;
-
- //- calculate any required properties
- virtual void calculateProperties() = 0;
-
- //- control molecules at different stages of the integration time-step
- virtual void controlMolsBeg() = 0;
-
- virtual void controlBeforeForces() = 0;
-
- virtual void controlMols() = 0;
-
- virtual void controlMolsEnd() = 0;
-
-
- //- output data
- virtual void output
- (
- const fileName& fixedPathName,
- const fileName& timePath
- ) = 0;
-
- //- update properties from a modified dictionary
- virtual void updateProperties(const dictionary&) = 0;
-
- void updateStateControllerProperties(const dictionary&);
-
- // Access
-
- //- return the control zone cells
- const labelList& controlZone() const;
-
- //- return the control zone name
- const word& regionName() const;
-
- //- return the targeted fields
- scalar density() const;
- scalar& density();
-
- const vector& velocity() const;
- vector& velocity();
-
- scalar temperature() const;
- scalar& temperature();
-
- scalar pressure() const;
- scalar& pressure();
-
- const tensor& strainRate() const;
- tensor& strainRate();
-
- const vector& tempGradient() const;
- vector& tempGradient();
-
- //- return the targeted fields
- const scalarField& densityField() const;
- scalarField& densityField();
-
- const vectorField& velocityField() const;
- vectorField& velocityField();
-
- const scalarField& temperatureField() const;
- scalarField& temperatureField();
-
- const scalarField& pressureField() const;
- scalarField& pressureField();
-
- bool singleValueController() const;
- bool& singleValueController();
-
- bool fieldController() const;
- bool& fieldController();
-
- bool writeInTimeDir() const;
- bool writeInCase() const;
-
- scalar avReqDensity();
- vector avReqVelocity();
- scalar avReqTemperature();
- scalar avReqPressure();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C b/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C
deleted file mode 100644
index 08ea263dd6..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C
+++ /dev/null
@@ -1,114 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see , 0);
-}
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-bool Foam::monoatomic::move
-(
- monoatomic::trackingData& td,
- const scalar trackTime
-)
-{
- td.switchProcessor = false;
- td.keepParticle = true;
-
- if (special_ != SPECIAL_FROZEN)
- {
- return td.keepParticle;
- }
-
- const constantProperties& constProps(td.cloud().constProps(id_));
-
- if (td.part() == trackingData::tpFirstVelocityHalfStep)
- {
- // First leapfrog velocity adjust part, required before tracking+force
- // part
-
- v_ += 0.5*trackTime*a_;
- }
- else if (td.part() == trackingData::tpLinearTrack)
- {
- // Leapfrog tracking part
-
- scalar tEnd = (1.0 - stepFraction())*trackTime;
- scalar dtMax = tEnd;
-
- while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
- {
- // set the lagrangian time-step
- scalar dt = min(dtMax, tEnd);
-
- dt *= trackToFace(position() + dt*v_, td);
-
- tEnd -= dt;
- stepFraction() = 1.0 - tEnd/trackTime;
- }
-
- setSitePositions(constProps);
- }
- else if (td.part() == trackingData::tpSecondVelocityHalfStep)
- {
- // Second leapfrog velocity adjust part, required after tracking+force
- // part
-
- a_ = siteForces_[0]/constProps.mass();
-
- v_ += 0.5*trackTime*a_;
- }
- else if (td.part() != trackingData::tpRotationalTrack)
- {
- FatalErrorIn("monoatomic::move(trackingData&, const scalar)") << nl
- << td.part() << " is an invalid part of the integration method."
- << abort(FatalError);
- }
-
- return td.keepParticle;
-}
-
-
-void Foam::monoatomic::transformProperties(const tensor& T)
-{
- particle::transformProperties(T);
-
- v_ = transform(T, v_);
-
- a_ = transform(T, a_);
-
- rf_ = transform(T, rf_);
-
- sitePositions_[0] = position_ + (T & (sitePositions_[0] - position_));
-
- siteForces_[0] = T & siteForces_[0];
-}
-
-
-void Foam::monoatomic::transformProperties(const vector& separation)
-{
- particle::transformProperties(separation);
-
- if (special_ == SPECIAL_TETHERED)
- {
- specialPosition_ += separation;
- }
-
- sitePositions_[0] += separation;
-}
-
-
-void Foam::monoatomic::setSitePositions(const constantProperties& constProps)
-{
- sitePositions_[0] = position_;
-}
-
-
-void Foam::monoatomic::setSiteSizes(label size)
-{
- // Nothing required, size controlled internally
-}
-
-
-bool Foam::monoatomic::hitPatch
-(
- const polyPatch&,
- trackingData&,
- const label,
- const scalar,
- const tetIndices&
-)
-{
- return false;
-}
-
-
-void Foam::monoatomic::hitProcessorPatch
-(
- const processorPolyPatch&,
- trackingData& td
-)
-{
- td.switchProcessor = true;
-}
-
-
-void Foam::monoatomic::hitWallPatch
-(
- const wallPolyPatch& wpp,
- trackingData& td,
- const tetIndices& tetIs
-)
-{
- // Use of the normal from tetIs is not required as
- // hasWallImpactDistance is false.
- vector nw = normal();
- nw /= mag(nw);
-
- scalar vn = v_ & nw;
-
- // Specular reflection
- if (vn > 0)
- {
- v_ -= 2*vn*nw;
- }
-}
-
-
-void Foam::monoatomic::hitPatch
-(
- const polyPatch&,
- trackingData& td
-)
-{
- td.keepParticle = false;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.H b/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.H
deleted file mode 100644
index d1e999435b..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.H
+++ /dev/null
@@ -1,418 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see sites_;
-
- //- Which sites require electrostatic interactions
- FixedList electrostaticSites_;
-
- //- Which sites require pair interactions
- FixedList pairPotSites_;
-
- //-
- scalar mass_;
-
-
- public:
-
- //-
- inline constantProperties();
-
- //- Construct from dictionary
- inline constantProperties
- (
- const dictionary& dict,
- const List& siteIds
- );
-
- // Member functions
-
- //-
- inline const FixedList& sites() const;
-
- //-
- inline const FixedList& pairPotSites() const;
-
- //-
- inline const FixedList& electrostaticSites() const;
-
- //-
- inline label degreesOfFreedom() const;
-
- //-
- inline scalar mass() const;
-
- //-
- inline label nSites() const;
- };
-
-
- //- Class used to pass tracking data to the trackToFace function
- class trackingData
- :
- public particle::TrackingData >
- {
- public:
-
- enum trackPart
- {
- tpFirstVelocityHalfStep,
- tpLinearTrack,
- tpRotationalTrack,
- tpSecondVelocityHalfStep,
- tpAccess
- };
-
- private:
-
- // label specifying which part of the integration algorithm is taking
- label part_;
-
-
- public:
-
- // Constructors
-
- trackingData(MoleculeCloud& cloud, trackPart part)
- :
- particle::TrackingData >(cloud),
- part_(part)
- {}
-
- // Member functions
-
- inline label part() const
- {
- return part_;
- }
- };
-
-
-private:
-
- // Private data
-
- //- Linear velocity of monoatomic
- vector v_;
-
- //- Total linear acceleration of monoatomic
- vector a_;
-
- //-
- vector specialPosition_;
-
- //-
- scalar potentialEnergy_;
-
- // - r_ij f_ij, stress dyad
- tensor rf_;
-
-// // - r_ij outer product f_ij: virial contribution
-// tensor rDotf_;
-
- //-
- label special_;
-
- //-
- label id_;
-
- //-
- FixedList siteForces_;
-
- //-
- FixedList sitePositions_;
-
-
-public:
-
- //- Runtime type information
- TypeName("monoatomic");
-
- friend class Cloud;
-
- // Constructors
-
- //- Construct with macroscopic description
- inline monoatomic
- (
- const polyMesh& mesh,
- const vector& position,
- const label cellI,
- const label tetFaceI,
- const label tetPtI,
- const scalar temperature,
- const vector& bulkVelocity,
- const vector& specialPosition,
- const constantProperties& constProps,
- trackingData& td,
- const label special,
- const label id
- );
-
- //- Construct from all components
- inline monoatomic
- (
- const polyMesh& mesh,
- const vector& position,
- const label cellI,
- const label tetFaceI,
- const label tetPtI,
- const vector& v,
- const vector& a,
- const vector& specialPosition,
- const constantProperties& constProps,
- const label special,
- const label id
- );
-
- //- Construct from Istream
- monoatomic
- (
- const polyMesh& mesh,
- Istream& is,
- bool readFields = true
- );
-
- //- Construct and return a clone
- autoPtr clone() const
- {
- return autoPtr(new monoatomic(*this));
- }
-
- //- Factory class to read-construct particles used for
- // parallel transfer
- class iNew
- {
- const polyMesh& mesh_;
-
- public:
-
- iNew(const polyMesh& mesh)
- :
- mesh_(mesh)
- {}
-
- autoPtr operator()(Istream& is) const
- {
- return autoPtr(new monoatomic(mesh_, is, true));
- }
- };
-
-
- // Member Functions
-
- // Tracking
-
- //-
- bool move(trackingData&, const scalar trackTime);
-
- //-
- virtual void transformProperties(const tensor& T);
-
- //-
- virtual void transformProperties(const vector& separation);
-
- //-
- void setSitePositions(const constantProperties& constProps);
-
- //-
- void setSiteSizes(label size);
-
-
- // Access
-
- //-
- inline const vector& v() const;
-
- //-
- inline vector& v();
-
- //-
- inline const vector& a() const;
-
- //-
- inline vector& a();
-
- //-
- inline const FixedList& siteForces() const;
-
- //-
- inline FixedList& siteForces();
-
- //-
- inline const FixedList& sitePositions() const;
-
- //-
- inline FixedList& sitePositions();
-
- //-
- inline const vector& specialPosition() const;
-
- //-
- inline vector& specialPosition();
-
- //-
- inline scalar potentialEnergy() const;
-
- //-
- inline scalar& potentialEnergy();
-
- //-
- inline const tensor& rf() const;
-
- //-
- inline tensor& rf();
-
- //-
- inline label special() const;
-
- //-
- inline bool tethered() const;
-
- //-
- inline label id() const;
-
-
- // Member Operators
-
- //- Overridable function to handle the particle hitting a patch
- // Executed before other patch-hitting functions
- bool hitPatch
- (
- const polyPatch&,
- trackingData& td,
- const label patchI,
- const scalar trackFraction,
- const tetIndices& tetIs
- );
-
- //- Overridable function to handle the particle hitting a processorPatch
- void hitProcessorPatch
- (
- const processorPolyPatch&,
- trackingData& td
- );
-
- //- Overridable function to handle the particle hitting a wallPatch
- void hitWallPatch
- (
- const wallPolyPatch&,
- trackingData& td,
- const tetIndices&
- );
-
- //- Overridable function to handle the particle hitting a polyPatch
- void hitPatch
- (
- const polyPatch&,
- trackingData& td
- );
-
-
- // I-O
-
- //- Read
- static void readFields(Cloud& mC);
-
- //- Write
- static void writeFields(const Cloud& mC);
-
- //- Show info
- static void info(trackingData& td);
-
-
- // IOstream Operators
-
- friend Ostream& operator<<(Ostream&, const monoatomic&);
-};
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "monoatomicI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicI.H b/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicI.H
deleted file mode 100644
index 4485d53b8d..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicI.H
+++ /dev/null
@@ -1,328 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see & siteIds
-)
-:
- sites_(),
- electrostaticSites_(-1),
- pairPotSites_(-1),
- mass_(0)
-{
- if (siteIds.size() != 1)
- {
- FatalErrorIn
- (
- "Foam::monoatomic::constantProperties::constantProperties"
- "("
- "const dictionary& dict, "
- "const List& siteIds"
- ")"
- )
- << "monoatomic, single site only, given: " << dict
- << nl << abort(FatalError);
- }
-
- FixedList siteIdNames = dict.lookup("siteIds");
-
- FixedList siteReferencePositions
- (
- dict.lookup("siteReferencePositions")
- );
-
- FixedList siteMasses(dict.lookup("siteMasses"));
-
- FixedList siteCharges(dict.lookup("siteCharges"));
-
- FixedList pairPotentialIds(dict.lookup("pairPotentialSiteIds"));
-
- constPropSite site = sites_[0];
-
- site = constPropSite
- (
- siteReferencePositions[0],
- siteMasses[0],
- siteCharges[0],
- siteIds[0],
- siteIdNames[0],
- (findIndex(pairPotentialIds, siteIdNames[0]) != -1), // pair
- (mag(siteCharges[0]) > VSMALL) // charge
- );
-
- mass_ = site.siteMass();
-
- if (site.pairPotentialSite())
- {
- pairPotSites_[0] = 0;
- }
-
- if (site.electrostaticSite())
- {
- electrostaticSites_[0] = 0;
- }
-
- if (!site.pairPotentialSite() && !site.electrostaticSite())
- {
- WarningIn
- (
- "Foam::monoatomic::constantProperties::constantProperties"
- "("
- "const dictionary& dict"
- ")"
- )
- << siteIdNames[0] << " is a non-interacting site." << endl;
- }
-
- // Single site monoatomic - no rotational motion.
- site.siteReferencePosition() = vector::zero;
-}
-
-
-inline Foam::monoatomic::monoatomic
-(
- const polyMesh& mesh,
- const vector& position,
- const label cellI,
- const label tetFaceI,
- const label tetPtI,
- const scalar temperature,
- const vector& bulkVelocity,
- const vector& specialPosition,
- const constantProperties& cP,
- monoatomic::trackingData& td,
- const label special,
- const label id
-)
-:
- particle(mesh, position, cellI, tetFaceI, tetPtI),
- v_(bulkVelocity),
- a_(vector::zero),
- specialPosition_(specialPosition),
- potentialEnergy_(0.0),
- rf_(tensor::zero),
- special_(special),
- id_(id),
- siteForces_(vector::zero),
- sitePositions_()
-{
- setSitePositions(cP);
-
- v_ += td.cloud().equipartitionLinearVelocity(temperature, cP.mass());
-}
-
-
-inline Foam::monoatomic::monoatomic
-(
- const polyMesh& mesh,
- const vector& position,
- const label cellI,
- const label tetFaceI,
- const label tetPtI,
- const vector& v,
- const vector& a,
- const vector& specialPosition,
- const constantProperties& cP,
- const label special,
- const label id
-
-)
-:
- particle(mesh, position, cellI, tetFaceI, tetPtI),
- v_(v),
- a_(a),
- specialPosition_(specialPosition),
- potentialEnergy_(0.0),
- rf_(tensor::zero),
- special_(special),
- id_(id),
- siteForces_(vector::zero),
- sitePositions_()
-{
- setSitePositions(cP);
-}
-
-// * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
-
-inline const Foam::FixedList&
-Foam::monoatomic::constantProperties::sites() const
-{
- return sites_;
-}
-
-
-inline const Foam::FixedList&
-Foam::monoatomic::constantProperties::pairPotSites() const
-{
- return pairPotSites_;
-}
-
-
-inline const Foam::FixedList&
-Foam::monoatomic::constantProperties::electrostaticSites() const
-{
- return electrostaticSites_;
-}
-
-
-inline Foam::label
-Foam::monoatomic::constantProperties::degreesOfFreedom() const
-{
- return 3;
-}
-
-
-inline Foam::scalar Foam::monoatomic::constantProperties::mass() const
-{
- return mass_;
-}
-
-
-inline Foam::label Foam::monoatomic::constantProperties::nSites() const
-{
- return 1;
-}
-
-
-// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
-
-inline const Foam::vector& Foam::monoatomic::v() const
-{
- return v_;
-}
-
-
-inline Foam::vector& Foam::monoatomic::v()
-{
- return v_;
-}
-
-
-inline const Foam::vector& Foam::monoatomic::a() const
-{
- return a_;
-}
-
-
-inline Foam::vector& Foam::monoatomic::a()
-{
- return a_;
-}
-
-
-inline const Foam::FixedList&
-Foam::monoatomic::siteForces() const
-{
- return siteForces_;
-}
-
-
-inline Foam::FixedList& Foam::monoatomic::siteForces()
-{
- return siteForces_;
-}
-
-
-inline const Foam::FixedList&
-Foam::monoatomic::sitePositions() const
-{
- return sitePositions_;
-}
-
-
-inline Foam::FixedList& Foam::monoatomic::sitePositions()
-{
- return sitePositions_;
-}
-
-
-inline const Foam::vector& Foam::monoatomic::specialPosition() const
-{
- return specialPosition_;
-}
-
-
-inline Foam::vector& Foam::monoatomic::specialPosition()
-{
- return specialPosition_;
-}
-
-
-inline Foam::scalar Foam::monoatomic::potentialEnergy() const
-{
- return potentialEnergy_;
-}
-
-
-inline Foam::scalar& Foam::monoatomic::potentialEnergy()
-{
- return potentialEnergy_;
-}
-
-
-inline const Foam::tensor& Foam::monoatomic::rf() const
-{
- return rf_;
-}
-
-
-inline Foam::tensor& Foam::monoatomic::rf()
-{
- return rf_;
-}
-
-
-inline Foam::label Foam::monoatomic::special() const
-{
- return special_;
-}
-
-
-inline bool Foam::monoatomic::tethered() const
-{
- return special_ == SPECIAL_TETHERED;
-}
-
-
-inline Foam::label Foam::monoatomic::id() const
-{
- return id_;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicIO.C b/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicIO.C
deleted file mode 100644
index a9a27bb992..0000000000
--- a/src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicIO.C
+++ /dev/null
@@ -1,305 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see (&v_),
- sizeof(v_)
- + sizeof(a_)
- + sizeof(specialPosition_)
- + sizeof(potentialEnergy_)
- + sizeof(rf_)
- + sizeof(special_)
- + sizeof(id_)
- );
-
- is >> siteForces_ >> sitePositions_;
- }
- }
-
- // Check state of Istream
- is.check
- (
- "Foam::monoatomic::monoatomic"
- "("
- "const polyMesh& mesh,"
- "Istream& is,"
- "bool readFields"
- ")"
- );
-}
-
-
-void Foam::monoatomic::readFields(Cloud& mC)
-{
- if (!mC.size())
- {
- return;
- }
-
- particle::readFields(mC);
-
- IOField v(mC.fieldIOobject("v", IOobject::MUST_READ));
- mC.checkFieldIOobject(mC, v);
-
- IOField a(mC.fieldIOobject("a", IOobject::MUST_READ));
- mC.checkFieldIOobject(mC, a);
-
- IOField specialPosition
- (
- mC.fieldIOobject("specialPosition", IOobject::MUST_READ)
- );
- mC.checkFieldIOobject(mC, specialPosition);
-
- IOField special(mC.fieldIOobject("special", IOobject::MUST_READ));
- mC.checkFieldIOobject(mC, special);
-
- IOField id(mC.fieldIOobject("id", IOobject::MUST_READ));
- mC.checkFieldIOobject(mC, id);
-
- label i = 0;
-
- forAllIter(typename Cloud, mC, iter)
- {
- monoatomic& mol = iter();
-
- mol.v_ = v[i];
- mol.a_ = a[i];
- mol.specialPosition_ = specialPosition[i];
- mol.special_ = special[i];
- mol.id_ = id[i];
- i++;
- }
-}
-
-
-void Foam::monoatomic::writeFields(const Cloud& mC)
-{
- particle::writeFields(mC);
-
- label np = mC.size();
-
- IOField v(mC.fieldIOobject("v", IOobject::NO_READ), np);
- IOField a(mC.fieldIOobject("a", IOobject::NO_READ), np);
- IOField specialPosition
- (
- mC.fieldIOobject("specialPosition", IOobject::NO_READ),
- np
- );
- IOField special(mC.fieldIOobject("special", IOobject::NO_READ), np);
- IOField id(mC.fieldIOobject("id", IOobject::NO_READ), np);
-
- label i = 0;
- forAllConstIter(typename Cloud, mC, iter)
- {
- const monoatomic& mol = iter();
-
- v[i] = mol.v_;
- a[i] = mol.a_;
- specialPosition[i] = mol.specialPosition_;
- special[i] = mol.special_;
- id[i] = mol.id_;
- i++;
- }
-
- v.write();
- a.write();
- specialPosition.write();
- special.write();
- id.write();
-}
-
-
-void Foam::monoatomic::info(monoatomic::trackingData& td)
-{
- vector totalLinearMomentum(vector::zero);
- scalar maxVelocityMag = 0.0;
- scalar totalMass = 0.0;
- scalar totalLinearKE = 0.0;
- scalar totalPE = 0.0;
- scalar totalrDotf = 0.0;
- label nMols = td.cloud().size();
-
- forAllConstIter(typename Cloud, td.cloud(), mol)
- {
- const label molId = mol().id();
- scalar molMass(td.cloud().constProps(molId).mass());
- totalMass += molMass;
- }
-
- forAllConstIter(typename Cloud, td.cloud(), mol)
- {
- const label molId = mol().id();
- const monoatomic::constantProperties cP
- (
- td.cloud().constProps(molId)
- );
- scalar molMass(cP.mass());
- const vector& molV(mol().v());
- totalLinearMomentum += molV * molMass;
- if (mag(molV) > maxVelocityMag)
- {
- maxVelocityMag = mag(molV);
- }
- totalLinearKE += 0.5*molMass*magSqr(molV);
- totalPE += mol().potentialEnergy();
- totalrDotf += tr(mol().rf());
- }
-
- scalar meshVolume = sum(td.cloud().mesh().cellVolumes());
-
- if (Pstream::parRun())
- {
- reduce(totalLinearMomentum, sumOp());
- reduce(maxVelocityMag, maxOp());
- reduce(totalMass, sumOp());
- reduce(totalLinearKE, sumOp());
- reduce(totalPE, sumOp