Merge branch 'cvm' of /home/noisy3/OpenFOAM/OpenFOAM-dev into cvm

This commit is contained in:
mattijs
2011-06-29 11:22:20 +01:00
100 changed files with 1793 additions and 1496 deletions

View File

@ -1,7 +1,6 @@
EXE_INC = \
-I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
@ -11,6 +10,4 @@ EXE_LIBS = \
-lfiniteVolume \
-llagrangian \
-lmolecule \
-lpotential \
-lmolecularMeasurements
-lpotential

View File

@ -0,0 +1,17 @@
Info<< "\nReading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
potential pot(mesh);
polyatomicCloud molecules(mesh, pot);

View File

@ -30,59 +30,25 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "md.H"
#include "polyatomicCloud.H"
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nReading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
potential pot(mesh);
moleculeCloud molecules(mesh, pot);
#include "temperatureAndPressureVariables.H"
label nAveragingSteps = 0;
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
nAveragingSteps++;
Info<< "Time = " << runTime.timeName() << endl;
molecules.evolve();
#include "meanMomentumEnergyAndNMols.H"
#include "temperatureAndPressure.H"
runTime.write();
if (runTime.outputTime())
{
nAveragingSteps = 0;
}
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;

View File

@ -4,7 +4,6 @@ makeType=${1:-libso}
set -x
wmake $makeType potential
wmake $makeType molecularMeasurements
wmake $makeType molecule
# ----------------------------------------------------------------- end-of-file

View File

@ -1,9 +1,8 @@
reducedUnits/reducedUnits.C
reducedUnits/reducedUnitsIO.C
constPropSite/constPropSite.C
molecule/molecule.C
molecule/moleculeIO.C
polyatomic/polyatomic/polyatomic.C
polyatomic/polyatomic/polyatomicIO.C
moleculeCloud/moleculeCloud.C
polyatomic/polyatomicCloud/polyatomicCloud.C
LIB = $(FOAM_LIBBIN)/libmolecule

View File

@ -0,0 +1,76 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
Class
constPropSite
Description
\*----------------------------------------------------------------------------*/
#include "constPropSite.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::constPropSite::constPropSite()
:
siteReferencePosition_(vector::zero),
siteMass_(0.0),
siteCharge_(0.0),
siteId_(0),
name_(),
pairPotentialSite_(false),
electrostaticSite_(false)
{}
Foam::constPropSite::constPropSite
(
const vector& siteReferencePosition,
const scalar& siteMass,
const scalar& siteCharge,
const label& siteId,
const word& name,
const bool& pairPotentialSite,
const bool& electrostaticSite
)
:
siteReferencePosition_(siteReferencePosition),
siteMass_(siteMass),
siteCharge_(siteCharge),
siteId_(siteId),
name_(name),
pairPotentialSite_(pairPotentialSite),
electrostaticSite_(electrostaticSite)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::constPropSite::~constPropSite()
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -0,0 +1,162 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
Class
Description
SourceFiles
constPropSiteI.H
constPropSite.C
\*---------------------------------------------------------------------------*/
#ifndef constPropSite_H
#define constPropSite_H
#include "vector.H"
#include "IFstream.H"
#include "OFstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class constPropSite Declaration
\*---------------------------------------------------------------------------*/
class constPropSite
{
// Private data
//-
vector siteReferencePosition_;
//-
scalar siteMass_;
//-
scalar siteCharge_;
//-
label siteId_;
//-
word name_;
//-
bool pairPotentialSite_;
//-
bool electrostaticSite_;
public:
// Constructors
//- Construct null
constPropSite();
//- Construct from components
constPropSite
(
const vector& siteReferencePosition,
const scalar& siteMass,
const scalar& siteCharge,
const label& siteId,
const word& name,
const bool& pairPotentialSite,
const bool& electrostaticSite
);
// Selectors
// Destructor
~constPropSite();
// Member Functions
// Access
//-
inline const vector& siteReferencePosition() const;
//-
inline vector& siteReferencePosition();
//-
inline const scalar& siteMass() const;
//-
inline scalar& siteMass();
//-
inline const scalar& siteCharge() const;
//-
inline scalar& siteCharge();
//-
inline const label& siteId() const;
//-
inline label& siteId();
//-
inline const word& name() const;
//-
inline word& name();
//-
inline const bool& pairPotentialSite() const;
//-
inline bool& pairPotentialSite();
//-
inline const bool& electrostaticSite() const;
//-
inline bool& electrostaticSite();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "constPropSiteI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,117 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
Class
Description
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::vector& Foam::constPropSite::siteReferencePosition() const
{
return siteReferencePosition_;
}
inline Foam::vector& Foam::constPropSite::siteReferencePosition()
{
return siteReferencePosition_;
}
inline const Foam::scalar& Foam::constPropSite::siteMass() const
{
return siteMass_;
}
inline Foam::scalar& Foam::constPropSite::siteMass()
{
return siteMass_;
}
inline const Foam::scalar& Foam::constPropSite::siteCharge() const
{
return siteCharge_;
}
inline Foam::scalar& Foam::constPropSite::siteCharge()
{
return siteCharge_;
}
inline const Foam::label& Foam::constPropSite::siteId() const
{
return siteId_;
}
inline Foam::label& Foam::constPropSite::siteId()
{
return siteId_;
}
inline const Foam::word& Foam::constPropSite::name() const
{
return name_;
}
inline Foam::word& Foam::constPropSite::name()
{
return name_;
}
inline const bool& Foam::constPropSite::pairPotentialSite() const
{
return pairPotentialSite_;
}
inline bool& Foam::constPropSite::pairPotentialSite()
{
return pairPotentialSite_;
}
inline const bool& Foam::constPropSite::electrostaticSite() const
{
return electrostaticSite_;
}
inline bool& Foam::constPropSite::electrostaticSite()
{
return electrostaticSite_;
}
// ************************************************************************* //

View File

@ -1,606 +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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
inline Foam::molecule::constantProperties::constantProperties()
:
siteReferencePositions_(Field<vector>(0)),
siteMasses_(List<scalar>(0)),
siteCharges_(List<scalar>(0)),
siteIds_(List<label>(0)),
pairPotentialSites_(List<bool>(false)),
electrostaticSites_(List<bool>(false)),
momentOfInertia_(diagTensor(0, 0, 0)),
mass_(0)
{}
inline Foam::molecule::constantProperties::constantProperties
(
const dictionary& dict
)
:
siteReferencePositions_(dict.lookup("siteReferencePositions")),
siteMasses_(dict.lookup("siteMasses")),
siteCharges_(dict.lookup("siteCharges")),
siteIds_(List<word>(dict.lookup("siteIds")).size(), -1),
pairPotentialSites_(),
electrostaticSites_(),
momentOfInertia_(),
mass_()
{
checkSiteListSizes();
setInteracionSiteBools
(
List<word>(dict.lookup("siteIds")),
List<word>(dict.lookup("pairPotentialSiteIds"))
);
mass_ = sum(siteMasses_);
vector centreOfMass(vector::zero);
// Calculate the centre of mass of the body and subtract it from each
// position
forAll(siteReferencePositions_, i)
{
centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
}
centreOfMass /= mass_;
siteReferencePositions_ -= centreOfMass;
if (siteIds_.size() == 1)
{
// Single site molecule - no rotational motion.
siteReferencePositions_[0] = vector::zero;
momentOfInertia_ = diagTensor(-1, -1, -1);
}
else if (linearMoleculeTest())
{
// Linear molecule.
Info<< nl << "Linear molecule." << endl;
vector dir = siteReferencePositions_[1] - siteReferencePositions_[0];
dir /= mag(dir);
tensor Q = rotationTensor(dir, vector(1,0,0));
siteReferencePositions_ = (Q & siteReferencePositions_);
// The rotation was around the centre of mass but remove any
// components that have crept in due to floating point errors
centreOfMass = vector::zero;
forAll(siteReferencePositions_, i)
{
centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
}
centreOfMass /= mass_;
siteReferencePositions_ -= centreOfMass;
diagTensor momOfInertia = diagTensor::zero;
forAll(siteReferencePositions_, i)
{
const vector& p(siteReferencePositions_[i]);
momOfInertia +=
siteMasses_[i]*diagTensor(0, p.x()*p.x(), p.x()*p.x());
}
momentOfInertia_ = diagTensor
(
-1,
momOfInertia.yy(),
momOfInertia.zz()
);
}
else
{
// Fully 6DOF molecule
// Calculate the inertia tensor in the current orientation
tensor momOfInertia(tensor::zero);
forAll(siteReferencePositions_, i)
{
const vector& p(siteReferencePositions_[i]);
momOfInertia += siteMasses_[i]*tensor
(
p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
-p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
-p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
);
}
if (eigenValues(momOfInertia).x() < VSMALL)
{
FatalErrorIn("molecule::constantProperties::constantProperties")
<< "An eigenvalue of the inertia tensor is zero, the molecule "
<< dict.name() << " is not a valid 6DOF shape."
<< nl << abort(FatalError);
}
// Normalise the inertia tensor magnitude to avoid SMALL numbers in the
// components causing problems
momOfInertia /= eigenValues(momOfInertia).x();
tensor e = eigenVectors(momOfInertia);
// Calculate the transformation between the principle axes to the
// global axes
tensor Q =
vector(1,0,0)*e.x() + vector(0,1,0)*e.y() + vector(0,0,1)*e.z();
// Transform the site positions
siteReferencePositions_ = (Q & siteReferencePositions_);
// Recalculating the moment of inertia with the new site positions
// The rotation was around the centre of mass but remove any
// components that have crept in due to floating point errors
centreOfMass = vector::zero;
forAll(siteReferencePositions_, i)
{
centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
}
centreOfMass /= mass_;
siteReferencePositions_ -= centreOfMass;
// Calculate the moment of inertia in the principle component
// reference frame
momOfInertia = tensor::zero;
forAll(siteReferencePositions_, i)
{
const vector& p(siteReferencePositions_[i]);
momOfInertia += siteMasses_[i]*tensor
(
p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
-p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
-p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
);
}
momentOfInertia_ = diagTensor
(
momOfInertia.xx(),
momOfInertia.yy(),
momOfInertia.zz()
);
}
}
inline Foam::molecule::molecule
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const tensor& Q,
const vector& v,
const vector& a,
const vector& pi,
const vector& tau,
const vector& specialPosition,
const constantProperties& constProps,
const label special,
const label id
)
:
particle(mesh, position, cellI, tetFaceI, tetPtI),
Q_(Q),
v_(v),
a_(a),
pi_(pi),
tau_(tau),
specialPosition_(specialPosition),
potentialEnergy_(0.0),
rf_(tensor::zero),
special_(special),
id_(id),
siteForces_(constProps.nSites(), vector::zero),
sitePositions_(constProps.nSites())
{
setSitePositions(constProps);
}
// * * * * * * * constantProperties Private Member Functions * * * * * * * * //
inline void Foam::molecule::constantProperties::checkSiteListSizes() const
{
if
(
siteIds_.size() != siteReferencePositions_.size()
|| siteIds_.size() != siteCharges_.size()
)
{
FatalErrorIn("molecule::constantProperties::checkSiteListSizes")
<< "Sizes of site id, charge and "
<< "referencePositions are not the same. "
<< nl << abort(FatalError);
}
}
inline void Foam::molecule::constantProperties::setInteracionSiteBools
(
const List<word>& siteIds,
const List<word>& pairPotSiteIds
)
{
pairPotentialSites_.setSize(siteIds_.size());
electrostaticSites_.setSize(siteIds_.size());
forAll(siteIds_, i)
{
const word& id(siteIds[i]);
pairPotentialSites_[i] = (findIndex(pairPotSiteIds, id) > -1);
electrostaticSites_[i] = (mag(siteCharges_[i]) > VSMALL);
}
}
inline bool Foam::molecule::constantProperties::linearMoleculeTest() const
{
if (siteIds_.size() == 2)
{
return true;
}
vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0];
refDir /= mag(refDir);
for
(
label i = 2;
i < siteReferencePositions_.size();
i++
)
{
vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1];
dir /= mag(dir);
if (mag(refDir & dir) < 1 - SMALL)
{
return false;
}
}
return true;
}
// * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
inline const Foam::Field<Foam::vector>&
Foam::molecule::constantProperties::siteReferencePositions() const
{
return siteReferencePositions_;
}
inline const Foam::List<Foam::scalar>&
Foam::molecule::constantProperties::siteMasses() const
{
return siteMasses_;
}
inline const Foam::List<Foam::scalar>&
Foam::molecule::constantProperties::siteCharges() const
{
return siteCharges_;
}
inline const Foam::List<Foam::label>&
Foam::molecule::constantProperties::siteIds() const
{
return siteIds_;
}
inline Foam::List<Foam::label>&
Foam::molecule::constantProperties::siteIds()
{
return siteIds_;
}
inline const Foam::List<bool>&
Foam::molecule::constantProperties::pairPotentialSites() const
{
return pairPotentialSites_;
}
inline bool Foam::molecule::constantProperties::pairPotentialSite
(
label sId
) const
{
label s = findIndex(siteIds_, sId);
if (s == -1)
{
FatalErrorIn("moleculeI.H") << nl
<< sId << " site not found."
<< nl << abort(FatalError);
}
return pairPotentialSites_[s];
}
inline const Foam::List<bool>&
Foam::molecule::constantProperties::electrostaticSites() const
{
return electrostaticSites_;
}
inline bool Foam::molecule::constantProperties::electrostaticSite
(
label sId
) const
{
label s = findIndex(siteIds_, sId);
if (s == -1)
{
FatalErrorIn
(
"molecule::constantProperties::electrostaticSite(label)"
) << sId << " site not found."
<< nl << abort(FatalError);
}
return electrostaticSites_[s];
}
inline const Foam::diagTensor&
Foam::molecule::constantProperties::momentOfInertia() const
{
return momentOfInertia_;
}
inline bool Foam::molecule::constantProperties::linearMolecule() const
{
return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
}
inline bool Foam::molecule::constantProperties::pointMolecule() const
{
return (momentOfInertia_.zz() < 0);
}
inline Foam::label Foam::molecule::constantProperties::degreesOfFreedom() const
{
if (linearMolecule())
{
return 5;
}
else if (pointMolecule())
{
return 3;
}
else
{
return 6;
}
}
inline Foam::scalar Foam::molecule::constantProperties::mass() const
{
return mass_;
}
inline Foam::label Foam::molecule::constantProperties::nSites() const
{
return siteIds_.size();
}
// * * * * * * * * * * * * molecule Member Functions * * * * * * * * * * * * //
inline const Foam::tensor& Foam::molecule::Q() const
{
return Q_;
}
inline Foam::tensor& Foam::molecule::Q()
{
return Q_;
}
inline const Foam::vector& Foam::molecule::v() const
{
return v_;
}
inline Foam::vector& Foam::molecule::v()
{
return v_;
}
inline const Foam::vector& Foam::molecule::a() const
{
return a_;
}
inline Foam::vector& Foam::molecule::a()
{
return a_;
}
inline const Foam::vector& Foam::molecule::pi() const
{
return pi_;
}
inline Foam::vector& Foam::molecule::pi()
{
return pi_;
}
inline const Foam::vector& Foam::molecule::tau() const
{
return tau_;
}
inline Foam::vector& Foam::molecule::tau()
{
return tau_;
}
inline const Foam::List<Foam::vector>& Foam::molecule::siteForces() const
{
return siteForces_;
}
inline Foam::List<Foam::vector>& Foam::molecule::siteForces()
{
return siteForces_;
}
inline const Foam::List<Foam::vector>& Foam::molecule::sitePositions() const
{
return sitePositions_;
}
inline Foam::List<Foam::vector>& Foam::molecule::sitePositions()
{
return sitePositions_;
}
inline const Foam::vector& Foam::molecule::specialPosition() const
{
return specialPosition_;
}
inline Foam::vector& Foam::molecule::specialPosition()
{
return specialPosition_;
}
inline Foam::scalar Foam::molecule::potentialEnergy() const
{
return potentialEnergy_;
}
inline Foam::scalar& Foam::molecule::potentialEnergy()
{
return potentialEnergy_;
}
inline const Foam::tensor& Foam::molecule::rf() const
{
return rf_;
}
inline Foam::tensor& Foam::molecule::rf()
{
return rf_;
}
inline Foam::label Foam::molecule::special() const
{
return special_;
}
inline bool Foam::molecule::tethered() const
{
return special_ == SPECIAL_TETHERED;
}
inline Foam::label Foam::molecule::id() const
{
return id_;
}
// ************************************************************************* //

View File

@ -1,390 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "constants.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline void Foam::moleculeCloud::evaluatePair
(
molecule& molI,
molecule& molJ
)
{
const pairPotentialList& pairPot = pot_.pairPotentials();
const pairPotential& electrostatic = pairPot.electrostatic();
label idI = molI.id();
label idJ = molJ.id();
const molecule::constantProperties& constPropI(constProps(idI));
const molecule::constantProperties& constPropJ(constProps(idJ));
List<label> siteIdsI = constPropI.siteIds();
List<label> siteIdsJ = constPropJ.siteIds();
List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
List<bool> electrostaticSitesI = constPropI.electrostaticSites();
List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
forAll(siteIdsI, sI)
{
label idsI(siteIdsI[sI]);
forAll(siteIdsJ, sJ)
{
label idsJ(siteIdsJ[sJ]);
if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
{
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;
}
}
if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
{
vector rsIsJ =
molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
scalar rsIsJMagSq = magSqr(rsIsJ);
if (rsIsJMagSq <= electrostatic.rCutSqr())
{
scalar rsIsJMag = mag(rsIsJ);
scalar chargeI = constPropI.siteCharges()[sI];
scalar chargeJ = constPropJ.siteCharges()[sJ];
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;
}
}
}
}
}
inline bool Foam::moleculeCloud::evaluatePotentialLimit
(
molecule& molI,
molecule& molJ
) const
{
const pairPotentialList& pairPot = pot_.pairPotentials();
const pairPotential& electrostatic = pairPot.electrostatic();
label idI = molI.id();
label idJ = molJ.id();
const molecule::constantProperties& constPropI(constProps(idI));
const molecule::constantProperties& constPropJ(constProps(idJ));
List<label> siteIdsI = constPropI.siteIds();
List<label> siteIdsJ = constPropJ.siteIds();
List<bool> pairPotentialSitesI = constPropI.pairPotentialSites();
List<bool> electrostaticSitesI = constPropI.electrostaticSites();
List<bool> pairPotentialSitesJ = constPropJ.pairPotentialSites();
List<bool> electrostaticSitesJ = constPropJ.electrostaticSites();
forAll(siteIdsI, sI)
{
label idsI(siteIdsI[sI]);
forAll(siteIdsJ, sJ)
{
label idsJ(siteIdsJ[sJ]);
if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
{
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 molecules"
<< " 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;
};
}
}
if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
{
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.siteCharges()[sI];
scalar chargeJ = constPropJ.siteCharges()[sJ];
if
(
mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
> pot_.potentialEnergyLimit()
)
{
return true;
};
}
}
}
}
return false;
}
inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
(
scalar temperature,
scalar mass
)
{
return sqrt(physicoChemical::k.value()*temperature/mass)*vector
(
rndGen_.GaussNormal(),
rndGen_.GaussNormal(),
rndGen_.GaussNormal()
);
}
inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
(
scalar temperature,
const molecule::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()
);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::polyMesh& Foam::moleculeCloud::mesh() const
{
return mesh_;
}
inline const Foam::potential& Foam::moleculeCloud::pot() const
{
return pot_;
}
inline const Foam::List<Foam::DynamicList<Foam::molecule*> >&
Foam::moleculeCloud::cellOccupancy() const
{
return cellOccupancy_;
}
inline const Foam::InteractionLists<Foam::molecule>&
Foam::moleculeCloud::il() const
{
return il_;
}
inline const Foam::List<Foam::molecule::constantProperties>
Foam::moleculeCloud::constProps() const
{
return constPropList_;
}
inline const Foam::molecule::constantProperties&
Foam::moleculeCloud::constProps(label id) const
{
return constPropList_[id];
}
inline Foam::Random& Foam::moleculeCloud::rndGen()
{
return rndGen_;
}
// ************************************************************************* //

View File

@ -23,14 +23,14 @@ License
\*----------------------------------------------------------------------------*/
#include "moleculeCloud.H"
#include "molecule.H"
#include "polyatomicCloud.H"
#include "polyatomic.H"
#include "Random.H"
#include "Time.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::tensor Foam::molecule::rotationTensorX(scalar phi) const
Foam::tensor Foam::polyatomic::rotationTensorX(scalar phi) const
{
return tensor
(
@ -41,7 +41,7 @@ Foam::tensor Foam::molecule::rotationTensorX(scalar phi) const
}
Foam::tensor Foam::molecule::rotationTensorY(scalar phi) const
Foam::tensor Foam::polyatomic::rotationTensorY(scalar phi) const
{
return tensor
(
@ -52,7 +52,7 @@ Foam::tensor Foam::molecule::rotationTensorY(scalar phi) const
}
Foam::tensor Foam::molecule::rotationTensorZ(scalar phi) const
Foam::tensor Foam::polyatomic::rotationTensorZ(scalar phi) const
{
return tensor
(
@ -65,11 +65,20 @@ Foam::tensor Foam::molecule::rotationTensorZ(scalar phi) const
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime)
bool Foam::polyatomic::move
(
polyatomic::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() == 0)
@ -151,13 +160,15 @@ bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime)
tau_ = vector::zero;
forAll(siteForces_, s)
forAll(siteForces_, sI)
{
const vector& f = siteForces_[s];
const vector& f = siteForces_[sI];
a_ += f/m;
tau_ += (constProps.siteReferencePositions()[s] ^ (Q_.T() & f));
tau_ +=
constProps.sites()[sI].siteReferencePosition()
^ (Q_.T() & f);
}
v_ += 0.5*trackTime*a_;
@ -180,7 +191,7 @@ bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime)
}
else
{
FatalErrorIn("molecule::move(trackingData&, const scalar)") << nl
FatalErrorIn("polyatomic::move(trackingData&, const scalar)") << nl
<< td.part() << " is an invalid part of the integration method."
<< abort(FatalError);
}
@ -189,7 +200,7 @@ bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime)
}
void Foam::molecule::transformProperties(const tensor& T)
void Foam::polyatomic::transformProperties(const tensor& T)
{
particle::transformProperties(T);
@ -211,7 +222,7 @@ void Foam::molecule::transformProperties(const tensor& T)
}
void Foam::molecule::transformProperties(const vector& separation)
void Foam::polyatomic::transformProperties(const vector& separation)
{
particle::transformProperties(separation);
@ -224,13 +235,18 @@ void Foam::molecule::transformProperties(const vector& separation)
}
void Foam::molecule::setSitePositions(const constantProperties& constProps)
void Foam::polyatomic::setSitePositions(const constantProperties& constProps)
{
sitePositions_ = position_ + (Q_ & constProps.siteReferencePositions());
forAll(constProps.sites(), sI)
{
sitePositions_[sI] =
position_
+ (Q_ & constProps.sites()[sI].siteReferencePosition());
}
}
void Foam::molecule::setSiteSizes(label size)
void Foam::polyatomic::setSiteSizes(label size)
{
sitePositions_.setSize(size);
@ -238,7 +254,7 @@ void Foam::molecule::setSiteSizes(label size)
}
bool Foam::molecule::hitPatch
bool Foam::polyatomic::hitPatch
(
const polyPatch&,
trackingData&,
@ -251,7 +267,7 @@ bool Foam::molecule::hitPatch
}
void Foam::molecule::hitProcessorPatch
void Foam::polyatomic::hitProcessorPatch
(
const processorPolyPatch&,
trackingData& td
@ -261,7 +277,7 @@ void Foam::molecule::hitProcessorPatch
}
void Foam::molecule::hitWallPatch
void Foam::polyatomic::hitWallPatch
(
const wallPolyPatch& wpp,
trackingData& td,
@ -269,7 +285,7 @@ void Foam::molecule::hitWallPatch
)
{
// Use of the normal from tetIs is not required as
// hasWallImpactDistance for a moleculeCloud is false.
// hasWallImpactDistance for a polyatomicCloud is false.
vector nw = normal();
nw /= mag(nw);
@ -283,7 +299,7 @@ void Foam::molecule::hitWallPatch
}
void Foam::molecule::hitPatch
void Foam::polyatomic::hitPatch
(
const polyPatch&,
trackingData& td

View File

@ -22,25 +22,26 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::molecule
Foam::polyatomic
Description
Foam::molecule
Foam::polyatomic
SourceFiles
moleculeI.H
molecule.C
moleculeIO.C
polyatomicI.H
polyatomic.C
polyatomicIO.C
\*---------------------------------------------------------------------------*/
#ifndef molecule_H
#define molecule_H
#ifndef polyatomic_H
#define polyatomic_H
#include "particle.H"
#include "IOstream.H"
#include "autoPtr.H"
#include "diagTensor.H"
#include "constPropSite.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -48,13 +49,13 @@ namespace Foam
{
// Class forward declarations
class moleculeCloud;
class polyatomicCloud;
/*---------------------------------------------------------------------------*\
Class molecule Declaration
Class polyatomic Declaration
\*---------------------------------------------------------------------------*/
class molecule
class polyatomic
:
public particle
{
@ -73,23 +74,21 @@ public:
SPECIAL_USER = 1
};
//- Class to hold molecule constant properties
//- Class to hold polyatomic constant properties
class constantProperties
{
// Private data
Field<vector> siteReferencePositions_;
//- Sites of mass, charge or interaction
List<constPropSite> sites_;
List<scalar> siteMasses_;
//- Which sites require electrostatic interactions
List<label> electrostaticSites_;
List<scalar> siteCharges_;
List<label> siteIds_;
List<bool> pairPotentialSites_;
List<bool> electrostaticSites_;
//- Which sites require pair interactions
List<label> pairPotSites_;
//- Moment of intertia (in principal axis configiration)
diagTensor momentOfInertia_;
scalar mass_;
@ -97,8 +96,6 @@ public:
// Private Member Functions
void checkSiteListSizes() const;
void setInteracionSiteBools
(
const List<word>& siteIds,
@ -113,27 +110,19 @@ public:
inline constantProperties();
//- Construct from dictionary
inline constantProperties(const dictionary& dict);
inline constantProperties
(
const dictionary& dict,
const labelList& siteIds
);
// Member functions
inline const Field<vector>& siteReferencePositions() const;
inline const List<constPropSite>& sites() const;
inline const List<scalar>& siteMasses() const;
inline const List<label>& pairPotSites() const;
inline const List<scalar>& siteCharges() const;
inline const List<label>& siteIds() const;
inline List<label>& siteIds();
inline const List<bool>& pairPotentialSites() const;
inline bool pairPotentialSite(label sId) const;
inline const List<bool>& electrostaticSites() const;
inline bool electrostaticSite(label sId) const;
inline const List<label>& electrostaticSites() const;
inline const diagTensor& momentOfInertia() const;
@ -152,7 +141,7 @@ public:
//- Class used to pass tracking data to the trackToFace function
class trackingData
:
public particle::TrackingData<moleculeCloud>
public particle::TrackingData<polyatomicCloud>
{
// label specifying which part of the integration algorithm is taking
label part_;
@ -162,9 +151,9 @@ public:
// Constructors
trackingData(moleculeCloud& cloud, label part)
trackingData(polyatomicCloud& cloud, label part)
:
particle::TrackingData<moleculeCloud>(cloud),
particle::TrackingData<polyatomicCloud>(cloud),
part_(part)
{}
@ -181,14 +170,22 @@ private:
// Private data
//- Orientation, stored as the rotation tensor to transform
// from the polyatomic to the global reference frame, i.e.:
// globalVector = Q_ & bodyLocalVector
// bodyLocalVector = Q_.T() & globalVector
tensor Q_;
// Linear velocity of polyatomic
vector v_;
// Total linear acceleration of polyatomic
vector a_;
//- Angular momentum of polyatomic, in body local reference frame
vector pi_;
//- Total torque on polyatomic, in body local reference frame
vector tau_;
vector specialPosition_;
@ -198,6 +195,9 @@ private:
// - r_ij f_ij, stress dyad
tensor rf_;
// // - r_ij outer product f_ij: virial contribution
// tensor rDotf_;
label special_;
label id_;
@ -218,12 +218,12 @@ private:
public:
friend class Cloud<molecule>;
friend class Cloud<polyatomic>;
// Constructors
//- Construct from components
inline molecule
inline polyatomic
(
const polyMesh& mesh,
const vector& position,
@ -242,7 +242,7 @@ public:
);
//- Construct from Istream
molecule
polyatomic
(
const polyMesh& mesh,
Istream& is,
@ -252,7 +252,7 @@ public:
//- Construct and return a clone
autoPtr<particle> clone() const
{
return autoPtr<particle>(new molecule(*this));
return autoPtr<particle>(new polyatomic(*this));
}
//- Factory class to read-construct particles used for
@ -268,9 +268,9 @@ public:
mesh_(mesh)
{}
autoPtr<molecule> operator()(Istream& is) const
autoPtr<polyatomic> operator()(Istream& is) const
{
return autoPtr<molecule>(new molecule(mesh_, is, true));
return autoPtr<polyatomic>(new polyatomic(mesh_, is, true));
}
};
@ -367,14 +367,14 @@ public:
// I-O
static void readFields(Cloud<molecule>& mC);
static void readFields(Cloud<polyatomic>& mC);
static void writeFields(const Cloud<molecule>& mC);
static void writeFields(const Cloud<polyatomic>& mC);
// IOstream Operators
friend Ostream& operator<<(Ostream&, const molecule&);
friend Ostream& operator<<(Ostream&, const polyatomic&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -383,7 +383,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "moleculeI.H"
#include "polyatomicI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -0,0 +1,592 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
inline Foam::polyatomic::constantProperties::constantProperties()
:
sites_(),
electrostaticSites_(),
pairPotSites_(),
momentOfInertia_(diagTensor(0, 0, 0)),
mass_(0)
{}
inline Foam::polyatomic::constantProperties::constantProperties
(
const dictionary& dict,
const labelList& siteIds
)
:
sites_(),
electrostaticSites_(),
pairPotSites_(),
momentOfInertia_(),
mass_(0)
{
List<word> siteIdNames = dict.lookup("siteIds");
List<vector> siteReferencePositions(dict.lookup("siteReferencePositions"));
List<scalar> siteMasses(dict.lookup("siteMasses"));
List<scalar> siteCharges(dict.lookup("siteCharges"));
List<word> pairPotentialIds(dict.lookup("pairPotentialSiteIds"));
if
(
( siteIdNames.size() != sites_.size() )
|| ( siteReferencePositions.size() != sites_.size() )
|| ( siteCharges.size() != sites_.size() )
)
{
FatalErrorIn
(
"Foam::polyatomic::constantProperties::constantProperties"
"("
"const dictionary& dict"
")"
)
<< "Sizes of site id, charge and "
<< "referencePositions are not the same: " << sites_.size()
<< nl << abort(FatalError);
}
sites_.setSize(siteIdNames.size());
electrostaticSites_.setSize(sites_.size(), -1);
pairPotSites_.setSize(sites_.size(), -1);
label pairI = 0;
label electrostaticI = 0;
forAll(sites_, sI)
{
sites_[sI] = constPropSite
(
siteReferencePositions[sI],
siteMasses[sI],
siteCharges[sI],
siteIds[sI],
siteIdNames[sI],
(findIndex(pairPotentialIds, siteIdNames[sI]) != -1), // pair
(mag(siteCharges[sI]) > VSMALL) // charge
);
if (sites_[sI].pairPotentialSite())
{
pairPotSites_[pairI++] = sI;
}
if (sites_[sI].electrostaticSite())
{
electrostaticSites_[electrostaticI++] = sI;
}
if (sites_[sI].pairPotentialSite() && !sites_[sI].electrostaticSite())
{
WarningIn
(
"Foam::polyatomic::constantProperties::constantProperties"
"("
"const dictionary& dict"
")"
)
<< siteIdNames[sI] << " is a non-interacting site." << endl;
}
}
pairPotSites_.setSize(pairI);
electrostaticSites_.setSize(electrostaticI);
// Calculate the centre of mass of the body and subtract it from each
// position
vector centreOfMass(vector::zero);
forAll(sites_, sI)
{
mass_ += sites_[sI].siteMass();
centreOfMass +=
sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
}
centreOfMass /= mass_;
forAll(sites_, sI)
{
sites_[sI].siteReferencePosition() -= centreOfMass;
}
if (sites_.size() == 1)
{
// Single site polyatomic - no rotational motion.
sites_[0].siteReferencePosition() = vector::zero;
momentOfInertia_ = diagTensor(-1, -1, -1);
}
else if (linearMoleculeTest())
{
// Linear polyatomic.
Info<< nl << "Linear polyatomic." << endl;
vector dir =
sites_[1].siteReferencePosition()
- sites_[0].siteReferencePosition();
dir /= mag(dir);
tensor Q = rotationTensor(dir, vector(1,0,0));
forAll(sites_, sI)
{
sites_[sI].siteReferencePosition() =
(Q & sites_[sI].siteReferencePosition());
}
// The rotation was around the centre of mass but remove any
// components that have crept in due to floating point errors
centreOfMass = vector::zero;
forAll(sites_, sI)
{
centreOfMass +=
sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
}
centreOfMass /= mass_;
forAll(sites_, sI)
{
sites_[sI].siteReferencePosition() -= centreOfMass;
}
diagTensor momOfInertia = diagTensor::zero;
forAll(sites_, sI)
{
const vector& p(sites_[sI].siteReferencePosition());
momOfInertia +=
sites_[sI].siteMass()*diagTensor(0, p.x()*p.x(), p.x()*p.x());
}
momentOfInertia_ = diagTensor
(
-1,
momOfInertia.yy(),
momOfInertia.zz()
);
}
else
{
// Fully 6DOF polyatomic
// Calculate the inertia tensor in the current orientation
tensor momOfInertia(tensor::zero);
forAll(sites_, sI)
{
const vector& p(sites_[sI].siteReferencePosition());
momOfInertia += sites_[sI].siteMass()*tensor
(
p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
-p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
-p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
);
}
if (eigenValues(momOfInertia).x() < VSMALL)
{
FatalErrorIn("polyatomic::constantProperties::constantProperties")
<< "An eigenvalue of the inertia tensor is zero, the "
<< "polyatomic " << dict.name() << " is not a valid 6DOF shape."
<< nl << abort(FatalError);
}
// Normalise the inertia tensor magnitude to avoid SMALL numbers in the
// components causing problems
momOfInertia /= eigenValues(momOfInertia).x();
tensor e = eigenVectors(momOfInertia);
// Calculate the transformation between the principle axes to the
// global axes
tensor Q =
vector(1,0,0)*e.x() + vector(0,1,0)*e.y() + vector(0,0,1)*e.z();
// Transform the site positions
forAll(sites_, sI)
{
sites_[sI].siteReferencePosition() =
(Q & sites_[sI].siteReferencePosition());
}
// Recalculating the moment of inertia with the new site positions
// The rotation was around the centre of mass but remove any
// components that have crept in due to floating point errors
centreOfMass = vector::zero;
forAll(sites_, sI)
{
centreOfMass +=
sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
}
centreOfMass /= mass_;
forAll(sites_, sI)
{
sites_[sI].siteReferencePosition() -= centreOfMass;
}
// Calculate the moment of inertia in the principle component
// reference frame
momOfInertia = tensor::zero;
forAll(sites_, sI)
{
const vector& p(sites_[sI].siteReferencePosition());
momOfInertia += sites_[sI].siteMass()*tensor
(
p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
-p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
-p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
);
}
momentOfInertia_ = diagTensor
(
momOfInertia.xx(),
momOfInertia.yy(),
momOfInertia.zz()
);
}
}
inline Foam::polyatomic::polyatomic
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const tensor& Q,
const vector& v,
const vector& a,
const vector& pi,
const vector& tau,
const vector& specialPosition,
const constantProperties& constProps,
const label special,
const label id
)
:
particle(mesh, position, cellI, tetFaceI, tetPtI),
Q_(Q),
v_(v),
a_(a),
pi_(pi),
tau_(tau),
specialPosition_(specialPosition),
potentialEnergy_(0.0),
rf_(tensor::zero),
special_(special),
id_(id),
siteForces_(constProps.nSites(), vector::zero),
sitePositions_(constProps.nSites())
{
setSitePositions(constProps);
}
// * * * * * * * constantProperties Private Member Functions * * * * * * * * //
inline bool Foam::polyatomic::constantProperties::linearMoleculeTest() const
{
if (sites_.size() == 2)
{
return true;
}
vector refDir =
sites_[1].siteReferencePosition()
- sites_[0].siteReferencePosition();
refDir /= mag(refDir);
for
(
label i = 2;
i < sites_.size();
i++
)
{
vector dir =
sites_[i].siteReferencePosition()
- sites_[i - 1].siteReferencePosition();
dir /= mag(dir);
if (mag(refDir & dir) < 1 - SMALL)
{
return false;
}
}
return true;
}
// * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
inline const Foam::List<Foam::constPropSite>&
Foam::polyatomic::constantProperties::sites() const
{
return sites_;
}
inline const Foam::List<Foam::label>&
Foam::polyatomic::constantProperties::pairPotSites() const
{
return pairPotSites_;
}
inline const Foam::List<Foam::label>&
Foam::polyatomic::constantProperties::electrostaticSites() const
{
return electrostaticSites_;
}
inline const Foam::diagTensor&
Foam::polyatomic::constantProperties::momentOfInertia() const
{
return momentOfInertia_;
}
inline bool Foam::polyatomic::constantProperties::linearMolecule() const
{
return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
}
inline bool Foam::polyatomic::constantProperties::pointMolecule() const
{
return (momentOfInertia_.zz() < 0);
}
inline Foam::label
Foam::polyatomic::constantProperties::degreesOfFreedom() const
{
if (linearMolecule())
{
return 5;
}
else if (pointMolecule())
{
return 3;
}
else
{
return 6;
}
}
inline Foam::scalar Foam::polyatomic::constantProperties::mass() const
{
return mass_;
}
inline Foam::label Foam::polyatomic::constantProperties::nSites() const
{
return sites_.size();
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
inline const Foam::tensor& Foam::polyatomic::Q() const
{
return Q_;
}
inline Foam::tensor& Foam::polyatomic::Q()
{
return Q_;
}
inline const Foam::vector& Foam::polyatomic::v() const
{
return v_;
}
inline Foam::vector& Foam::polyatomic::v()
{
return v_;
}
inline const Foam::vector& Foam::polyatomic::a() const
{
return a_;
}
inline Foam::vector& Foam::polyatomic::a()
{
return a_;
}
inline const Foam::vector& Foam::polyatomic::pi() const
{
return pi_;
}
inline Foam::vector& Foam::polyatomic::pi()
{
return pi_;
}
inline const Foam::vector& Foam::polyatomic::tau() const
{
return tau_;
}
inline Foam::vector& Foam::polyatomic::tau()
{
return tau_;
}
inline const Foam::List<Foam::vector>& Foam::polyatomic::siteForces() const
{
return siteForces_;
}
inline Foam::List<Foam::vector>& Foam::polyatomic::siteForces()
{
return siteForces_;
}
inline const Foam::List<Foam::vector>& Foam::polyatomic::sitePositions() const
{
return sitePositions_;
}
inline Foam::List<Foam::vector>& Foam::polyatomic::sitePositions()
{
return sitePositions_;
}
inline const Foam::vector& Foam::polyatomic::specialPosition() const
{
return specialPosition_;
}
inline Foam::vector& Foam::polyatomic::specialPosition()
{
return specialPosition_;
}
inline Foam::scalar Foam::polyatomic::potentialEnergy() const
{
return potentialEnergy_;
}
inline Foam::scalar& Foam::polyatomic::potentialEnergy()
{
return potentialEnergy_;
}
inline const Foam::tensor& Foam::polyatomic::rf() const
{
return rf_;
}
inline Foam::tensor& Foam::polyatomic::rf()
{
return rf_;
}
inline Foam::label Foam::polyatomic::special() const
{
return special_;
}
inline bool Foam::polyatomic::tethered() const
{
return special_ == SPECIAL_TETHERED;
}
inline Foam::label Foam::polyatomic::id() const
{
return id_;
}
// ************************************************************************* //

View File

@ -23,13 +23,13 @@ License
\*---------------------------------------------------------------------------*/
#include "molecule.H"
#include "polyatomic.H"
#include "IOstreams.H"
#include "moleculeCloud.H"
#include "polyatomicCloud.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::molecule::molecule
Foam::polyatomic::polyatomic
(
const polyMesh& mesh,
Istream& is,
@ -91,13 +91,13 @@ Foam::molecule::molecule
// Check state of Istream
is.check
(
"Foam::molecule::molecule"
"(const Cloud<molecule>& cloud, Foam::Istream&), bool"
"Foam::polyatomic::polyatomic"
"(const Cloud<polyatomic>& cloud, Foam::Istream&), bool"
);
}
void Foam::molecule::readFields(Cloud<molecule>& mC)
void Foam::polyatomic::readFields(Cloud<polyatomic>& mC)
{
if (!mC.size())
{
@ -134,9 +134,9 @@ void Foam::molecule::readFields(Cloud<molecule>& mC)
mC.checkFieldIOobject(mC, id);
label i = 0;
forAllIter(moleculeCloud, mC, iter)
forAllIter(polyatomicCloud, mC, iter)
{
molecule& mol = iter();
polyatomic& mol = iter();
mol.Q_ = Q[i];
mol.v_ = v[i];
@ -151,7 +151,7 @@ void Foam::molecule::readFields(Cloud<molecule>& mC)
}
void Foam::molecule::writeFields(const Cloud<molecule>& mC)
void Foam::polyatomic::writeFields(const Cloud<polyatomic>& mC)
{
particle::writeFields(mC);
@ -203,9 +203,9 @@ void Foam::molecule::writeFields(const Cloud<molecule>& mC)
);
label i = 0;
forAllConstIter(moleculeCloud, mC, iter)
forAllConstIter(polyatomicCloud, mC, iter)
{
const molecule& mol = iter();
const polyatomic& mol = iter();
Q[i] = mol.Q_;
v[i] = mol.v_;
@ -244,13 +244,13 @@ void Foam::molecule::writeFields(const Cloud<molecule>& mC)
Info<< "writeFields " << mC.name() << endl;
if (isA<moleculeCloud>(mC))
if (isA<polyatomicCloud>(mC))
{
const moleculeCloud& m = dynamic_cast<const moleculeCloud&>(mC);
const polyatomicCloud& m = dynamic_cast<const polyatomicCloud&>(mC);
m.writeXYZ
(
m.mesh().time().timePath()/cloud::prefix/"moleculeCloud.xmol"
m.mesh().time().timePath()/cloud::prefix/"polyatomicCloud.xmol"
);
}
}
@ -258,7 +258,7 @@ void Foam::molecule::writeFields(const Cloud<molecule>& mC)
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<(Ostream& os, const molecule& mol)
Foam::Ostream& Foam::operator<<(Ostream& os, const polyatomic& mol)
{
if (os.format() == IOstream::ASCII)
{
@ -302,7 +302,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const molecule& mol)
os.check
(
"Foam::Ostream& Foam::operator<<"
"(Foam::Ostream&, const Foam::molecule&)"
"(Foam::Ostream&, const Foam::polyatomic&)"
);
return os;

View File

@ -23,7 +23,7 @@ License
\*----------------------------------------------------------------------------*/
#include "moleculeCloud.H"
#include "polyatomicCloud.H"
#include "fvMesh.H"
#include "mathematicalConstants.H"
@ -33,14 +33,14 @@ using namespace Foam::constant::mathematical;
namespace Foam
{
defineTemplateTypeNameAndDebug(Cloud<molecule>, 0);
defineTemplateTypeNameAndDebug(Cloud<polyatomic>, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::moleculeCloud::buildConstProps()
void Foam::polyatomicCloud::buildConstProps()
{
Info<< nl << "Reading moleculeProperties dictionary." << endl;
Info<< nl << "Reading polyatomicProperties dictionary." << endl;
const List<word>& idList(pot_.idList());
@ -48,11 +48,11 @@ void Foam::moleculeCloud::buildConstProps()
const List<word>& siteIdList(pot_.siteIdList());
IOdictionary moleculePropertiesDict
IOdictionary polyatomicPropertiesDict
(
IOobject
(
"moleculeProperties",
"polyatomicProperties",
mesh_.time().constant(),
mesh_,
IOobject::MUST_READ_IF_MODIFIED,
@ -64,7 +64,8 @@ void Foam::moleculeCloud::buildConstProps()
forAll(idList, i)
{
const word& id = idList[i];
const dictionary& molDict = moleculePropertiesDict.subDict(id);
const dictionary& molDict(polyatomicPropertiesDict.subDict(id));
List<word> siteIdNames = molDict.lookup("siteIds");
@ -78,26 +79,32 @@ void Foam::moleculeCloud::buildConstProps()
if (siteIds[sI] == -1)
{
FatalErrorIn("moleculeCloud::buildConstProps()")
FatalErrorIn
(
"Foam::polyatomic::constantProperties::constantProperties"
"("
"const dictionary& dict"
")"
)
<< siteId << " site not found."
<< nl << abort(FatalError);
}
}
molecule::constantProperties& constProp = constPropList_[i];
constProp = molecule::constantProperties(molDict);
constProp.siteIds() = siteIds;
constPropList_[i] = polyatomic::constantProperties
(
molDict,
siteIds
);
}
}
void Foam::moleculeCloud::setSiteSizesAndPositions()
void Foam::polyatomicCloud::setSiteSizesAndPositions()
{
forAllIter(moleculeCloud, *this, mol)
forAllIter(polyatomicCloud, *this, mol)
{
const molecule::constantProperties& cP = constProps(mol().id());
const polyatomic::constantProperties& cP = constProps(mol().id());
mol().setSiteSizes(cP.nSites());
@ -106,14 +113,14 @@ void Foam::moleculeCloud::setSiteSizesAndPositions()
}
void Foam::moleculeCloud::buildCellOccupancy()
void Foam::polyatomicCloud::buildCellOccupancy()
{
forAll(cellOccupancy_, cO)
{
cellOccupancy_[cO].clear();
}
forAllIter(moleculeCloud, *this, mol)
forAllIter(polyatomicCloud, *this, mol)
{
cellOccupancy_[mol().cell()].append(&mol());
}
@ -125,15 +132,15 @@ void Foam::moleculeCloud::buildCellOccupancy()
}
void Foam::moleculeCloud::calculatePairForce()
void Foam::polyatomicCloud::calculatePairForce()
{
PstreamBuffers pBufs(Pstream::nonBlocking);
// Start sending referred data
il_.sendReferredData(cellOccupancy(), pBufs);
molecule* molI = NULL;
molecule* molJ = NULL;
polyatomic* molI = NULL;
polyatomic* molJ = NULL;
{
// Real-Real interactions
@ -148,7 +155,7 @@ void Foam::moleculeCloud::calculatePairForce()
forAll(dil[d], interactingCells)
{
List<molecule*> cellJ =
List<polyatomic*> cellJ =
cellOccupancy_[dil[d][interactingCells]];
forAll(cellJ, cellJMols)
@ -180,24 +187,24 @@ void Foam::moleculeCloud::calculatePairForce()
const labelListList& ril = il_.ril();
List<IDLList<molecule> >& referredMols = il_.referredParticles();
List<IDLList<polyatomic> >& referredMols = il_.referredParticles();
forAll(ril, r)
{
const List<label>& realCells = ril[r];
IDLList<molecule>& refMols = referredMols[r];
IDLList<polyatomic>& refMols = referredMols[r];
forAllIter
(
IDLList<molecule>,
IDLList<polyatomic>,
refMols,
refMol
)
{
forAll(realCells, rC)
{
List<molecule*> cellI = cellOccupancy_[realCells[rC]];
List<polyatomic*> cellI = cellOccupancy_[realCells[rC]];
forAll(cellI, cellIMols)
{
@ -212,11 +219,11 @@ void Foam::moleculeCloud::calculatePairForce()
}
void Foam::moleculeCloud::calculateTetherForce()
void Foam::polyatomicCloud::calculateTetherForce()
{
const tetherPotentialList& tetherPot(pot_.tetherPotentials());
forAllIter(moleculeCloud, *this, mol)
forAllIter(polyatomicCloud, *this, mol)
{
if (mol().tethered())
{
@ -239,16 +246,16 @@ void Foam::moleculeCloud::calculateTetherForce()
}
void Foam::moleculeCloud::calculateExternalForce()
void Foam::polyatomicCloud::calculateExternalForce()
{
forAllIter(moleculeCloud, *this, mol)
forAllIter(polyatomicCloud, *this, mol)
{
mol().a() += pot_.gravity();
}
}
void Foam::moleculeCloud::removeHighEnergyOverlaps()
void Foam::polyatomicCloud::removeHighEnergyOverlaps()
{
Info<< nl << "Removing high energy overlaps, limit = "
<< pot_.potentialEnergyLimit()
@ -267,11 +274,11 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
// Real-Real interaction
molecule* molI = NULL;
molecule* molJ = NULL;
polyatomic* molI = NULL;
polyatomic* molJ = NULL;
{
DynamicList<molecule*> molsToDelete;
DynamicList<polyatomic*> molsToDelete;
const labelListList& dil(il_.dil());
@ -283,7 +290,7 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
forAll(dil[d], interactingCells)
{
List<molecule*> cellJ =
List<polyatomic*> cellJ =
cellOccupancy_[dil[d][interactingCells]];
forAll(cellJ, cellJMols)
@ -369,19 +376,19 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
// Real-Referred interaction
{
DynamicList<molecule*> molsToDelete;
DynamicList<polyatomic*> molsToDelete;
const labelListList& ril(il_.ril());
List<IDLList<molecule> >& referredMols = il_.referredParticles();
List<IDLList<polyatomic> >& referredMols = il_.referredParticles();
forAll(ril, r)
{
IDLList<molecule>& refMols = referredMols[r];
IDLList<polyatomic>& refMols = referredMols[r];
forAllIter
(
IDLList<molecule>,
IDLList<polyatomic>,
refMols,
refMol
)
@ -394,7 +401,7 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
{
label cellI = realCells[rC];
List<molecule*> cellIMols = cellOccupancy_[cellI];
List<polyatomic*> cellIMols = cellOccupancy_[cellI];
forAll(cellIMols, cIM)
{
@ -423,7 +430,7 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
== findIndex(pot_.removalOrder(), idJ)
)
{
// Remove one of the molecules
// Remove one of the polyatomics
// arbitrarily, assuring that a
// consistent decision is made for
// both real-referred pairs.
@ -463,24 +470,25 @@ void Foam::moleculeCloud::removeHighEnergyOverlaps()
reduce(molsRemoved, sumOp<label>());
}
Info<< tab << molsRemoved << " molecules removed" << endl;
Info<< tab << molsRemoved << " polyatomics removed" << endl;
}
void Foam::moleculeCloud::initialiseMolecules
void Foam::polyatomicCloud::initialiseMolecules
(
const IOdictionary& mdInitialiseDict
)
{
Info<< nl
<< "Initialising molecules in each zone specified in mdInitialiseDict."
<< "Initialising polyatomics in each zone specified in "
<< mdInitialiseDict.name()
<< endl;
const cellZoneMesh& cellZones = mesh_.cellZones();
if (!cellZones.size())
{
FatalErrorIn("void Foam::moleculeCloud::initialiseMolecules")
FatalErrorIn("void Foam::polyatomicCloud::initialiseMolecules")
<< "No cellZones found in the mesh."
<< abort(FatalError);
}
@ -520,7 +528,7 @@ void Foam::moleculeCloud::initialiseMolecules
if (latticeIds.size() != latticePositions.size())
{
FatalErrorIn("Foam::moleculeCloud::initialiseMolecules")
FatalErrorIn("Foam::polyatomicCloud::initialiseMolecules")
<< "latticeIds and latticePositions must be the same "
<< " size." << nl
<< abort(FatalError);
@ -542,7 +550,7 @@ void Foam::moleculeCloud::initialiseMolecules
if (numberDensity < VSMALL)
{
WarningIn("moleculeCloud::initialiseMolecules")
WarningIn("polyatomicCloud::initialiseMolecules")
<< "numberDensity too small, not filling zone "
<< zone.name() << endl;
@ -563,7 +571,10 @@ void Foam::moleculeCloud::initialiseMolecules
{
label id = findIndex(pot_.idList(), latticeIds[i]);
const molecule::constantProperties& cP(constProps(id));
const polyatomic::constantProperties& cP
(
constProps(id)
);
unitCellMass += cP.mass();
}
@ -575,7 +586,7 @@ void Foam::moleculeCloud::initialiseMolecules
if (massDensity < VSMALL)
{
WarningIn("moleculeCloud::initialiseMolecules")
WarningIn("polyatomicCloud::initialiseMolecules")
<< "massDensity too small, not filling zone "
<< zone.name() << endl;
@ -591,7 +602,7 @@ void Foam::moleculeCloud::initialiseMolecules
}
else
{
FatalErrorIn("Foam::moleculeCloud::initialiseMolecules")
FatalErrorIn("Foam::polyatomicCloud::initialiseMolecules")
<< "massDensity or numberDensity not specified " << nl
<< abort(FatalError);
}
@ -686,8 +697,8 @@ void Foam::moleculeCloud::initialiseMolecules
anchor += (R & (latticeCellShape & latticeAnchor));
// Continue trying to place molecules as long as at
// least one molecule is placed in each iteration.
// Continue trying to place polyatomics as long as at
// least one polyatomic is placed in each iteration.
// The "|| totalZoneMols == 0" condition means that the
// algorithm will continue if the origin is outside the
// zone.
@ -709,7 +720,7 @@ void Foam::moleculeCloud::initialiseMolecules
bool partOfLayerInBounds = false;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// start of placement of molecules
// start of placement of polyatomics
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (n == 0)
@ -949,7 +960,7 @@ void Foam::moleculeCloud::initialiseMolecules
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// end of placement of molecules
// end of placement of polyatomics
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if
@ -958,10 +969,13 @@ void Foam::moleculeCloud::initialiseMolecules
&& !partOfLayerInBounds
)
{
WarningIn("Foam::moleculeCloud::initialiseMolecules()")
WarningIn
(
"Foam::polyatomicCloud::initialiseMolecules()"
)
<< "A whole layer of unit cells was placed "
<< "outside the bounds of the mesh, but no "
<< "molecules have been placed in zone '"
<< "polyatomics have been placed in zone '"
<< zone.name()
<< "'. This is likely to be because the zone "
<< "has few cells ("
@ -987,7 +1001,7 @@ void Foam::moleculeCloud::initialiseMolecules
}
void Foam::moleculeCloud::createMolecule
void Foam::polyatomicCloud::createMolecule
(
const point& position,
label cell,
@ -1006,7 +1020,7 @@ void Foam::moleculeCloud::createMolecule
if (cell == -1)
{
FatalErrorIn("Foam::moleculeCloud::createMolecule")
FatalErrorIn("Foam::polyatomicCloud::createMolecule")
<< "Position specified does not correspond to a mesh cell." << nl
<< abort(FatalError);
}
@ -1019,10 +1033,10 @@ void Foam::moleculeCloud::createMolecule
{
specialPosition = position;
special = molecule::SPECIAL_TETHERED;
special = polyatomic::SPECIAL_TETHERED;
}
const molecule::constantProperties& cP(constProps(id));
const polyatomic::constantProperties& cP(constProps(id));
vector v = equipartitionLinearVelocity(temperature, cP.mass());
@ -1058,7 +1072,7 @@ void Foam::moleculeCloud::createMolecule
addParticle
(
new molecule
new polyatomic
(
mesh_,
position,
@ -1079,11 +1093,11 @@ void Foam::moleculeCloud::createMolecule
}
Foam::label Foam::moleculeCloud::nSites() const
Foam::label Foam::polyatomicCloud::nSites() const
{
label n = 0;
forAllConstIter(moleculeCloud, *this, mol)
forAllConstIter(polyatomicCloud, *this, mol)
{
n += constProps(mol().id()).nSites();
}
@ -1094,14 +1108,14 @@ Foam::label Foam::moleculeCloud::nSites() const
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::moleculeCloud::moleculeCloud
Foam::polyatomicCloud::polyatomicCloud
(
const polyMesh& mesh,
const potential& pot,
bool readFields
)
:
Cloud<molecule>(mesh, "moleculeCloud", false),
Cloud<polyatomic>(mesh, "polyatomicCloud", false),
mesh_(mesh),
pot_(pot),
cellOccupancy_(mesh_.nCells()),
@ -1111,7 +1125,7 @@ Foam::moleculeCloud::moleculeCloud
{
if (readFields)
{
molecule::readFields(*this);
polyatomic::readFields(*this);
}
buildConstProps();
@ -1124,7 +1138,7 @@ Foam::moleculeCloud::moleculeCloud
}
Foam::moleculeCloud::moleculeCloud
Foam::polyatomicCloud::polyatomicCloud
(
const polyMesh& mesh,
const potential& pot,
@ -1132,7 +1146,7 @@ Foam::moleculeCloud::moleculeCloud
bool readFields
)
:
Cloud<molecule>(mesh, "moleculeCloud", false),
Cloud<polyatomic>(mesh, "polyatomicCloud", false),
mesh_(mesh),
pot_(pot),
il_(mesh_, 0.0, false),
@ -1141,7 +1155,7 @@ Foam::moleculeCloud::moleculeCloud
{
if (readFields)
{
molecule::readFields(*this);
polyatomic::readFields(*this);
}
clear();
@ -1154,30 +1168,30 @@ Foam::moleculeCloud::moleculeCloud
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::moleculeCloud::evolve()
void Foam::polyatomicCloud::evolve()
{
molecule::trackingData td0(*this, 0);
Cloud<molecule>::move(td0, mesh_.time().deltaTValue());
polyatomic::trackingData td0(*this, 0);
Cloud<polyatomic>::move(td0, mesh_.time().deltaTValue());
molecule::trackingData td1(*this, 1);
Cloud<molecule>::move(td1, mesh_.time().deltaTValue());
polyatomic::trackingData td1(*this, 1);
Cloud<polyatomic>::move(td1, mesh_.time().deltaTValue());
molecule::trackingData td2(*this, 2);
Cloud<molecule>::move(td2, mesh_.time().deltaTValue());
polyatomic::trackingData td2(*this, 2);
Cloud<polyatomic>::move(td2, mesh_.time().deltaTValue());
calculateForce();
molecule::trackingData td3(*this, 3);
Cloud<molecule>::move(td3, mesh_.time().deltaTValue());
polyatomic::trackingData td3(*this, 3);
Cloud<polyatomic>::move(td3, mesh_.time().deltaTValue());
}
void Foam::moleculeCloud::calculateForce()
void Foam::polyatomicCloud::calculateForce()
{
buildCellOccupancy();
// Set accumulated quantities to zero
forAllIter(moleculeCloud, *this, mol)
forAllIter(polyatomicCloud, *this, mol)
{
mol().siteForces() = vector::zero;
@ -1194,7 +1208,7 @@ void Foam::moleculeCloud::calculateForce()
}
void Foam::moleculeCloud::applyConstraintsAndThermostats
void Foam::polyatomicCloud::applyConstraintsAndThermostats
(
const scalar targetTemperature,
const scalar measuredTemperature
@ -1214,7 +1228,7 @@ void Foam::moleculeCloud::applyConstraintsAndThermostats
<< "----------------------------------------"
<< endl;
forAllIter(moleculeCloud, *this, mol)
forAllIter(polyatomicCloud, *this, mol)
{
mol().v() *= temperatureCorrectionFactor;
@ -1223,21 +1237,22 @@ void Foam::moleculeCloud::applyConstraintsAndThermostats
}
void Foam::moleculeCloud::writeXYZ(const fileName& fName) const
void Foam::polyatomicCloud::writeXYZ(const fileName& fName) const
{
OFstream os(fName);
os << nSites() << nl << "moleculeCloud site positions in angstroms" << nl;
os << nSites() << nl
<< "polyatomicCloud site positions in angstroms" << nl;
forAllConstIter(moleculeCloud, *this, mol)
forAllConstIter(polyatomicCloud, *this, mol)
{
const molecule::constantProperties& cP = constProps(mol().id());
const polyatomic::constantProperties& cP = constProps(mol().id());
forAll(mol().sitePositions(), i)
{
const point& sP = mol().sitePositions()[i];
os << pot_.siteIdList()[cP.siteIds()[i]]
os << pot_.siteIdList()[cP.sites()[i].siteId()]
<< ' ' << sP.x()*1e10
<< ' ' << sP.y()*1e10
<< ' ' << sP.z()*1e10

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,22 +22,22 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::moleculeCloud
Foam::polyatomicCloud
Description
SourceFiles
moleculeCloudI.H
moleculeCloud.C
polyatomicCloudI.H
polyatomicCloud.C
\*---------------------------------------------------------------------------*/
#ifndef moleculeCloud_H
#define moleculeCloud_H
#ifndef polyatomicCloud_H
#define polyatomicCloud_H
#include "Cloud.H"
#include "molecule.H"
#include "polyatomic.H"
#include "IOdictionary.H"
#include "potential.H"
#include "InteractionLists.H"
@ -51,12 +51,12 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class moleculeCloud Declaration
Class polyatomicCloud Declaration
\*---------------------------------------------------------------------------*/
class moleculeCloud
class polyatomicCloud
:
public Cloud<molecule>
public Cloud<polyatomic>
{
private:
@ -67,11 +67,11 @@ private:
const potential& pot_;
List<DynamicList<molecule*> > cellOccupancy_;
List<DynamicList<polyatomic*> > cellOccupancy_;
InteractionLists<molecule> il_;
InteractionLists<polyatomic> il_;
List<molecule::constantProperties> constPropList_;
List<polyatomic::constantProperties> constPropList_;
Random rndGen_;
@ -82,21 +82,21 @@ private:
void setSiteSizesAndPositions();
//- Determine which molecules are in which cells
//- Determine which polyatomics are in which cells
void buildCellOccupancy();
void calculatePairForce();
inline void evaluatePair
(
molecule& molI,
molecule& molJ
polyatomic& molI,
polyatomic& molJ
);
inline bool evaluatePotentialLimit
(
molecule& molI,
molecule& molJ
polyatomic& molI,
polyatomic& molJ
) const;
void calculateTetherForce();
@ -133,14 +133,14 @@ private:
inline vector equipartitionAngularMomentum
(
scalar temperature,
const molecule::constantProperties& cP
const polyatomic::constantProperties& cP
);
//- Disallow default bitwise copy construct
moleculeCloud(const moleculeCloud&);
polyatomicCloud(const polyatomicCloud&);
//- Disallow default bitwise assignment
void operator=(const moleculeCloud&);
void operator=(const polyatomicCloud&);
public:
@ -148,7 +148,7 @@ public:
// Constructors
//- Construct given mesh and potential references
moleculeCloud
polyatomicCloud
(
const polyMesh& mesh,
const potential& pot,
@ -156,7 +156,7 @@ public:
);
//- Construct given mesh, potential and mdInitialiseDict
moleculeCloud
polyatomicCloud
(
const polyMesh& mesh,
const potential& pot,
@ -167,7 +167,7 @@ public:
// Member Functions
//- Evolve the molecules (move, calculate forces, control state etc)
//- Evolve the polyatomics (move, calculate forces, control state etc)
void evolve();
void calculateForce();
@ -185,13 +185,14 @@ public:
inline const potential& pot() const;
inline const List<DynamicList<molecule*> >& cellOccupancy() const;
inline const List<DynamicList<polyatomic*> >& cellOccupancy() const;
inline const InteractionLists<molecule>& il() const;
inline const InteractionLists<polyatomic>& il() const;
inline const List<molecule::constantProperties> constProps() const;
inline const List<polyatomic::constantProperties>
constProps() const;
inline const molecule::constantProperties&
inline const polyatomic::constantProperties&
constProps(label id) const;
inline Random& rndGen();
@ -199,7 +200,7 @@ public:
// Member Operators
//- Write molecule sites in XYZ format
//- Write polyatomic sites in XYZ format
void writeXYZ(const fileName& fName) const;
};
@ -210,7 +211,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "moleculeCloudI.H"
#include "polyatomicCloudI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -0,0 +1,383 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "constants.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline void Foam::polyatomicCloud::evaluatePair
(
polyatomic& molI,
polyatomic& molJ
)
{
const pairPotentialList& pairPot = pot_.pairPotentials();
const pairPotential& electrostatic = pairPot.electrostatic();
label idI = molI.id();
label idJ = molJ.id();
const polyatomic::constantProperties& constPropI(constProps(idI));
const polyatomic::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;
}
}
}
}
inline bool Foam::polyatomicCloud::evaluatePotentialLimit
(
polyatomic& molI,
polyatomic& molJ
) const
{
const pairPotentialList& pairPot = pot_.pairPotentials();
const pairPotential& electrostatic = pairPot.electrostatic();
label idI = molI.id();
label idJ = molJ.id();
const polyatomic::constantProperties& constPropI(constProps(idI));
const polyatomic::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("polyatomicCloud::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 polyatomics"
<< " twice. Removing one of the polyatomics."
<< 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("polyatomicCloud::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 polyatomics"
<< " twice. Removing one of the polyatomics."
<< 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;
}
inline Foam::vector Foam::polyatomicCloud::equipartitionLinearVelocity
(
scalar temperature,
scalar mass
)
{
return sqrt(physicoChemical::k.value()*temperature/mass)*vector
(
rndGen_.GaussNormal(),
rndGen_.GaussNormal(),
rndGen_.GaussNormal()
);
}
inline Foam::vector Foam::polyatomicCloud::equipartitionAngularMomentum
(
scalar temperature,
const polyatomic::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()
);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::polyMesh& Foam::polyatomicCloud::mesh() const
{
return mesh_;
}
inline const Foam::potential& Foam::polyatomicCloud::pot() const
{
return pot_;
}
inline const Foam::List<Foam::DynamicList<Foam::polyatomic*> >&
Foam::polyatomicCloud::cellOccupancy() const
{
return cellOccupancy_;
}
inline const Foam::InteractionLists<Foam::polyatomic>&
Foam::polyatomicCloud::il() const
{
return il_;
}
inline const Foam::List<Foam::polyatomic::constantProperties>
Foam::polyatomicCloud::constProps() const
{
return constPropList_;
}
inline const Foam::polyatomic::constantProperties&
Foam::polyatomicCloud::constProps(label id) const
{
return constPropList_[id];
}
inline Foam::Random& Foam::polyatomicCloud::rndGen()
{
return rndGen_;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

@ -67,7 +67,6 @@ SourceFiles
#include "wallPolyPatch.H"
#include "processorPolyPatch.H"
#include "zeroGradientFvPatchFields.H"
#include "globalIndex.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -780,8 +779,7 @@ private:
void indexDualVertices
(
pointField& pts,
PackedBoolList& boundaryPts,
const globalIndex& globalParallelPointIndices
PackedBoolList& boundaryPts
);
//- Re-index all of the the Delaunay cells

View File

@ -58,90 +58,73 @@ void Foam::conformalVoronoiMesh::calcDualMesh
++cit
)
{
if
(
!cit->vertex(0)->real()
|| !cit->vertex(1)->real()
|| !cit->vertex(2)->real()
|| !cit->vertex(3)->real()
)
{
cit->filterCount() = 100;
}
else
{
cit->filterCount() = 0;
}
}
// for
// (
// Delaunay::Finite_cells_iterator cit = finite_cells_begin();
// cit != finite_cells_end();
// ++cit
// )
// {
// if
// (
// !cit->vertex(0)->real()
// || !cit->vertex(1)->real()
// || !cit->vertex(2)->real()
// || !cit->vertex(3)->real()
// )
// {
// cit->filterCount() = 100;
// }
// else
// {
// cit->filterCount() = 0;
// }
// }
for
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
{
std::list<Cell_handle> cells;
incident_cells(vit, std::back_inserter(cells));
// for
// (
// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
// vit != finite_vertices_end();
// vit++
// )
// {
// std::list<Cell_handle> cells;
// incident_cells(vit, std::back_inserter(cells));
bool hasProcPt = false;
// bool hasProcPt = false;
for
(
std::list<Cell_handle>::iterator cit=cells.begin();
cit != cells.end();
++cit
)
{
if
(
!(*cit)->vertex(0)->real()
|| !(*cit)->vertex(1)->real()
|| !(*cit)->vertex(2)->real()
|| !(*cit)->vertex(3)->real()
)
{
hasProcPt = true;
// for
// (
// std::list<Cell_handle>::iterator cit=cells.begin();
// cit != cells.end();
// ++cit
// )
// {
// if
// (
// !(*cit)->vertex(0)->real()
// || !(*cit)->vertex(1)->real()
// || !(*cit)->vertex(2)->real()
// || !(*cit)->vertex(3)->real()
// )
// {
// hasProcPt = true;
break;
}
}
// break;
// }
// }
// if (hasProcPt)
// {
// for
// (
// std::list<Cell_handle>::iterator cit=cells.begin();
// cit != cells.end();
// ++cit
// )
// {
// (*cit)->filterCount() = 100;
// }
// }
// }
if (hasProcPt)
{
for
(
std::list<Cell_handle>::iterator cit=cells.begin();
cit != cells.end();
++cit
)
{
(*cit)->filterCount() = 100;
}
}
}
PackedBoolList boundaryPts(number_of_cells(), false);
globalIndex globalParallelPointIndices(number_of_cells());
indexDualVertices
(
points,
boundaryPts,
globalParallelPointIndices
);
indexDualVertices(points, boundaryPts);
{
// Ideally requires a no-risk face filtering to get rid of zero area
@ -199,12 +182,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
// Reindexing the Delaunay cells and regenerating the
// points resets the mesh to the starting condition.
indexDualVertices
(
points,
boundaryPts,
globalParallelPointIndices
);
indexDualVertices(points, boundaryPts);
{
Info<< nl << "Merging close points" << endl;
@ -1551,65 +1529,65 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
List<polyPatch*> patches(patchStarts.size());
// label nValidPatches = 0;
// forAll(patches, p)
// {
// if (patchTypes[p] == processorPolyPatch::typeName)
// {
// // Do not create empty processor patches
// if (patchSizes[p] > 0)
// {
// patches[nValidPatches] = new processorPolyPatch
// (
// patchNames[p],
// patchSizes[p],
// patchStarts[p],
// nValidPatches,
// pMesh.boundaryMesh(),
// Pstream::myProcNo(),
// procNeighbours[p]
// );
// nValidPatches++;
// }
// }
// else
// {
// patches[nValidPatches] = polyPatch::New
// (
// patchTypes[p],
// patchNames[p],
// patchSizes[p],
// patchStarts[p],
// nValidPatches,
// pMesh.boundaryMesh()
// ).ptr();
// nValidPatches++;
// }
// }
// patches.setSize(nValidPatches);
// pMesh.addPatches(patches);
Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
label nValidPatches = 0;
forAll(patches, p)
{
patches[p] = new polyPatch
(
patchNames[p],
patchSizes[p],
patchStarts[p],
p,
pMesh.boundaryMesh()
);
if (patchTypes[p] == processorPolyPatch::typeName)
{
// Do not create empty processor patches
if (patchSizes[p] > 0)
{
patches[nValidPatches] = new processorPolyPatch
(
patchNames[p],
patchSizes[p],
patchStarts[p],
nValidPatches,
pMesh.boundaryMesh(),
Pstream::myProcNo(),
procNeighbours[p]
);
nValidPatches++;
}
}
else
{
patches[nValidPatches] = polyPatch::New
(
patchTypes[p],
patchNames[p],
patchSizes[p],
patchStarts[p],
nValidPatches,
pMesh.boundaryMesh()
).ptr();
nValidPatches++;
}
}
pMesh.addPatches(patches, false);
patches.setSize(nValidPatches);
pMesh.addPatches(patches);
// Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
// forAll(patches, p)
// {
// patches[p] = new polyPatch
// (
// patchNames[p],
// patchSizes[p],
// patchStarts[p],
// p,
// pMesh.boundaryMesh()
// );
// }
// pMesh.addPatches(patches, false);
// pMesh.overrideCellCentres(cellCentres);
@ -1807,8 +1785,7 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
void Foam::conformalVoronoiMesh::indexDualVertices
(
pointField& pts,
PackedBoolList& boundaryPts,
const globalIndex& globalParallelPointIndices
PackedBoolList& boundaryPts
)
{
// Indexing Delaunay cells, which are the dual vertices
@ -1821,16 +1798,6 @@ void Foam::conformalVoronoiMesh::indexDualVertices
boundaryPts = false;
OFstream tetStr
(
runTime_.path()
/"processor_"
+ name(Pstream::myProcNo())
+ "_parallelTets.obj"
);
label pTetI = 0;
for
(
Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@ -1840,35 +1807,23 @@ void Foam::conformalVoronoiMesh::indexDualVertices
{
if (cit->internalOrBoundaryDualVertex())
{
if (cit->parallelDualVertex())
{
cit->cellIndex() = globalParallelPointIndices.toGlobal
(
dualVertI
);
drawDelaunayCell(tetStr, cit, pTetI++);
}
else
{
cit->cellIndex() = dualVertI;
if
(
!cit->vertex(0)->internalOrBoundaryPoint()
|| !cit->vertex(1)->internalOrBoundaryPoint()
|| !cit->vertex(2)->internalOrBoundaryPoint()
|| !cit->vertex(3)->internalOrBoundaryPoint()
)
{
// This is a boundary dual vertex
boundaryPts[dualVertI] = true;
}
}
cit->cellIndex() = dualVertI;
pts[dualVertI] = topoint(dual(cit));
if
(
!cit->vertex(0)->internalOrBoundaryPoint()
|| !cit->vertex(1)->internalOrBoundaryPoint()
|| !cit->vertex(2)->internalOrBoundaryPoint()
|| !cit->vertex(3)->internalOrBoundaryPoint()
)
{
// This is a boundary dual vertex
boundaryPts[dualVertI] = true;
}
dualVertI++;
}
else

View File

@ -423,10 +423,10 @@ void Foam::conformalVoronoiMesh::writeMesh
patches.setSize(nValidPatches);
// mesh.addFvPatches(patches);
mesh.addFvPatches(patches);
Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
mesh.addFvPatches(patches, false);
// Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
// mesh.addFvPatches(patches, false);
if (!mesh.write())
{
@ -456,7 +456,7 @@ void Foam::conformalVoronoiMesh::writeMesh
// cellCs.write();
// writeCellSizes(mesh);
writeCellSizes(mesh);
findRemainingProtrusionSet(mesh);
}