songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy 2013-10-31 15:00:26 +00:00
commit 154a0d9e73
37 changed files with 652 additions and 279 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
applications/solvers/multiphase/interFoam/alphaEqn.H
View File

@ -6,7 +6,8 @@
phic = min(interface.cAlpha()*phic, max(phic));
surfaceScalarField phir(phic*interface.nHatf());
if (pimple.firstIter() && MULESCorr)
//***HGW if (pimple.firstIter() && MULESCorr)
if (MULESCorr)
{
fvScalarMatrix alpha1Eqn
(

9
applications/test/ODE/Test-ODE.C
View File

@ -27,7 +27,7 @@ Description
#include "argList.H"
#include "IOmanip.H"
#include "ODE.H"
#include "ODESystem.H"
#include "ODESolver.H"
#include "RK.H"
@ -37,7 +37,7 @@ using namespace Foam;
class testODE
:
public ODE
public ODESystem
{
public:
@ -107,9 +107,13 @@ int main(int argc, char *argv[])
argList::validArgs.append("ODESolver");
argList args(argc, argv);
// Create the ODE system
testODE ode;
// Create the selected ODE system solver
autoPtr<ODESolver> odeSolver = ODESolver::New(args[1], ode);
// Initialise the ODE system fields
scalar xStart = 1.0;
scalarField yStart(ode.nEqns());
yStart[0] = ::Foam::j0(xStart);
@ -117,6 +121,7 @@ int main(int argc, char *argv[])
yStart[2] = ::Foam::jn(2, xStart);
yStart[3] = ::Foam::jn(3, xStart);
// Print the evolution of the solution and the time-step
scalarField dyStart(ode.nEqns());
ode.derivatives(xStart, yStart, dyStart);

6
applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
View File

@ -1253,6 +1253,9 @@ int main(int argc, char *argv[])
// Snap parameters
const snapParameters snapParams(snapDict);
// Layer addition parameters
const layerParameters layerParams(layerDict, mesh.boundaryMesh());
if (wantRefine)
{
@ -1346,9 +1349,6 @@ int main(int argc, char *argv[])
globalToSlavePatch
);
// Layer addition parameters
layerParameters layerParams(layerDict, mesh.boundaryMesh());
// Use the maxLocalCells from the refinement parameters
bool preBalance = returnReduce
(

4
applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
View File

@ -304,11 +304,13 @@ addLayersControls
// size of the refined cell outside layer (true) or absolute sizes (false).
relativeSizes true;
// Layer thickness specification. This can be specified in one of four ways
// Layer thickness specification. This can be specified in one of following
// ways:
// - expansionRatio and finalLayerThickness (cell nearest internal mesh)
// - expansionRatio and firstLayerThickness (cell on surface)
// - overall thickness and firstLayerThickness
// - overall thickness and finalLayerThickness
// - overall thickness and expansionRatio
// Expansion factor for layer mesh
expansionRatio 1.0;

14
applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
View File

@ -720,6 +720,20 @@ int main(int argc, char *argv[])
PtrList<faceSet> faceSets(fSetNames.size());
PtrList<pointSet> pointSets(pSetNames.size());
Info<< "Reconstructing sets:" << endl;
if (cSetNames.size())
{
Info<< " cellSets " << cSetNames.sortedToc() << endl;
}
if (fSetNames.size())
{
Info<< " faceSets " << fSetNames.sortedToc() << endl;
}
if (pSetNames.size())
{
Info<< " pointSets " << pSetNames.sortedToc() << endl;
}
// Load sets
forAll(procMeshes.meshes(), procI)
{

8
src/ODE/ODESolvers/KRR4/KRR4.C
View File

@ -55,7 +55,7 @@ const scalar
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::KRR4::KRR4(const ODE& ode)
Foam::KRR4::KRR4(const ODESystem& ode)
:
ODESolver(ode),
yTemp_(n_, 0.0),
@ -76,7 +76,7 @@ Foam::KRR4::KRR4(const ODE& ode)
void Foam::KRR4::solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,
@ -168,7 +168,7 @@ void Foam::KRR4::solve
(
"void Foam::KRR4::solve"
"("
"const ODE&, "
"const ODESystem&, "
"scalar&, "
"scalarField&, "
"scalarField&, "
@ -206,7 +206,7 @@ void Foam::KRR4::solve
(
"void Foam::KRR4::solve"
"("
"const ODE&, "
"const ODESystem&, "
"scalar&, "
"scalarField&, "
"scalarField&, "

4
src/ODE/ODESolvers/KRR4/KRR4.H
View File

@ -88,14 +88,14 @@ public:
// Constructors
//- Construct from ODE
KRR4(const ODE& ode);
KRR4(const ODESystem& ode);
// Member Functions
void solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,

6
src/ODE/ODESolvers/ODESolver/ODESolver.C
View File

@ -36,7 +36,7 @@ namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::ODESolver::ODESolver(const ODE& ode)
Foam::ODESolver::ODESolver(const ODESystem& ode)
:
n_(ode.nEqns()),
yScale_(n_),
@ -48,7 +48,7 @@ Foam::ODESolver::ODESolver(const ODE& ode)
void Foam::ODESolver::solve
(
const ODE& ode,
const ODESystem& ode,
const scalar xStart,
const scalar xEnd,
scalarField& y,
@ -102,7 +102,7 @@ void Foam::ODESolver::solve
FatalErrorIn
(
"ODESolver::solve"
"(const ODE& ode, const scalar xStart, const scalar xEnd,"
"(const ODESystem& ode, const scalar xStart, const scalar xEnd,"
"scalarField& yStart, const scalar eps, scalar& hEst) const"
) << "Too many integration steps"
<< exit(FatalError);

16
src/ODE/ODESolvers/ODESolver/ODESolver.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ Class
Foam::ODESolver
Description
Selection for ODE solver
Abstract base-class for ODE system solvers
SourceFiles
ODESolver.C
@ -35,7 +35,7 @@ SourceFiles
#ifndef ODESolver_H
#define ODESolver_H
#include "ODE.H"
#include "ODESystem.H"
#include "typeInfo.H"
#include "autoPtr.H"
@ -82,7 +82,7 @@ public:
autoPtr,
ODESolver,
ODE,
(const ODE& ode),
(const ODESystem& ode),
(ode)
);
@ -90,7 +90,7 @@ public:
// Constructors
//- Construct for given ODE
ODESolver(const ODE& ode);
ODESolver(const ODESystem& ode);
// Selectors
@ -99,7 +99,7 @@ public:
static autoPtr<ODESolver> New
(
const word& ODESolverTypeName,
const ODE& ode
const ODESystem& ode
);
@ -112,7 +112,7 @@ public:
virtual void solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,
@ -126,7 +126,7 @@ public:
virtual void solve
(
const ODE& ode,
const ODESystem& ode,
const scalar xStart,
const scalar xEnd,
scalarField& y,

7
src/ODE/ODESolvers/ODESolver/ODESolverNew.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -30,7 +30,7 @@ License
Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
(
const Foam::word& ODESolverTypeName,
const Foam::ODE& ode
const Foam::ODESystem& ode
)
{
Info<< "Selecting ODE solver " << ODESolverTypeName << endl;
@ -42,7 +42,8 @@ Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
{
FatalErrorIn
(
"ODESolver::New(const word& ODESolverTypeName, const ODE& ode)"
"ODESolver::New"
"(const word& ODESolverTypeName, const ODESystem& ode)"
) << "Unknown ODESolver type "
<< ODESolverTypeName << nl << nl
<< "Valid ODESolvers are : " << endl

6
src/ODE/ODESolvers/RK/RK.C
View File

@ -54,7 +54,7 @@ const scalar
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::RK::RK(const ODE& ode)
Foam::RK::RK(const ODESystem& ode)
:
ODESolver(ode),
yTemp_(n_, 0.0),
@ -72,7 +72,7 @@ Foam::RK::RK(const ODE& ode)
void Foam::RK::solve
(
const ODE& ode,
const ODESystem& ode,
const scalar x,
const scalarField& y,
const scalarField& dydx,
@ -142,7 +142,7 @@ void Foam::RK::solve
void Foam::RK::solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,

6
src/ODE/ODESolvers/RK/RK.H
View File

@ -81,14 +81,14 @@ public:
// Constructors
//- Construct from ODE
RK(const ODE& ode);
RK(const ODESystem& ode);
// Member Functions
void solve
(
const ODE& ode,
const ODESystem& ode,
const scalar x,
const scalarField& y,
const scalarField& dydx,
@ -100,7 +100,7 @@ public:
void solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,

4
src/ODE/ODESolvers/SIBS/SIBS.C
View File

@ -47,7 +47,7 @@ namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::SIBS::SIBS(const ODE& ode)
Foam::SIBS::SIBS(const ODESystem& ode)
:
ODESolver(ode),
a_(iMaxX_, 0.0),
@ -70,7 +70,7 @@ Foam::SIBS::SIBS(const ODE& ode)
void Foam::SIBS::solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,

6
src/ODE/ODESolvers/SIBS/SIBS.H
View File

@ -78,7 +78,7 @@ class SIBS
void SIMPR
(
const ODE& ode,
const ODESystem& ode,
const scalar xStart,
const scalarField& y,
const scalarField& dydx,
@ -110,14 +110,14 @@ public:
// Constructors
//- Construct from ODE
SIBS(const ODE& ode);
SIBS(const ODESystem& ode);
// Member Functions
void solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,

2
src/ODE/ODESolvers/SIBS/SIMPR.C
View File

@ -29,7 +29,7 @@ License
void Foam::SIBS::SIMPR
(
const ODE& ode,
const ODESystem& ode,
const scalar xStart,
const scalarField& y,
const scalarField& dydx,

25
src/ODE/ODE/ODE.H → src/ODE/ODESystem/ODESystem.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,15 +22,15 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::ODE
Foam::ODESystem
Description
Abstract base class for the ODE solvers.
Abstract base class for the systems of ordinary differential equations.
\*---------------------------------------------------------------------------*/
#ifndef ODE_H
#define ODE_H
#ifndef ODESystem_H
#define ODESystem_H
#include "scalarField.H"
#include "scalarMatrices.H"
@ -41,30 +41,32 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class ODE Declaration
Class ODESystem Declaration
\*---------------------------------------------------------------------------*/
class ODE
class ODESystem
{
public:
// Constructors
//- Construct null
ODE()
//- Construct null8
ODESystem()
{}
//- Destructor
virtual ~ODE()
virtual ~ODESystem()
{}
// Member Functions
//- Return the number of equations in the system
virtual label nEqns() const = 0;
//- Calculate the derivatives in dydx
virtual void derivatives
(
const scalar x,
@ -72,7 +74,8 @@ public:
scalarField& dydx
) const = 0;
//- Calculate the Jacobian of the system
// Need by the stiff-system solvers
virtual void jacobian
(
const scalar x,

7
src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
View File

@ -542,8 +542,11 @@ public:
//- Compact maps. Gets per field a bool whether it is used (locally)
// and works out itself what this side and sender side can remove
// from maps.
void compact(const boolList& elemIsUsed, const int tag);
void compact
(
const boolList& elemIsUsed,
const int tag = UPstream::msgType()
);
//- Distribute data. Note:schedule only used for Pstream::scheduled
// for now, all others just use send-to-all, receive-from-all.

2
src/OpenFOAM/primitives/Tuple2/Tuple2.H
View File

@ -25,7 +25,7 @@ Class
Foam::Tuple2
Description
A 2-tuple.
A 2-tuple for storing two objects of different types.
SeeAlso
Foam::Pair for storing two objects of identical types.

348
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
View File

@ -52,6 +52,7 @@ Description
#include "fixedValuePointPatchFields.H"
#include "calculatedPointPatchFields.H"
#include "cyclicSlipPointPatchFields.H"
#include "fixedValueFvPatchFields.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -110,6 +111,31 @@ void Foam::autoLayerDriver::dumpDisplacement
}
Foam::tmp<Foam::scalarField> Foam::autoLayerDriver::avgPointData
(
const indirectPrimitivePatch& pp,
const scalarField& pointFld
)
{
tmp<scalarField> tfaceFld(new scalarField(pp.size(), 0.0));
scalarField& faceFld = tfaceFld();
forAll(pp.localFaces(), faceI)
{
const face& f = pp.localFaces()[faceI];
if (f.size())
{
forAll(f, fp)
{
faceFld[faceI] += pointFld[f[fp]];
}
faceFld[faceI] /= f.size();
}
}
return tfaceFld;
}
// Check that primitivePatch is not multiply connected. Collect non-manifold
// points in pointSet.
void Foam::autoLayerDriver::checkManifold
@ -2391,19 +2417,21 @@ Foam::label Foam::autoLayerDriver::countExtrusion
}
// Collect layer faces and layer cells into bools for ease of handling
// Collect layer faces and layer cells into mesh fields for ease of handling
void Foam::autoLayerDriver::getLayerCellsFaces
(
const polyMesh& mesh,
const addPatchCellLayer& addLayer,
boolList& flaggedCells,
boolList& flaggedFaces
const scalarField& oldRealThickness,
labelList& cellNLayers,
scalarField& faceRealThickness
)
{
flaggedCells.setSize(mesh.nCells());
flaggedCells = false;
flaggedFaces.setSize(mesh.nFaces());
flaggedFaces = false;
cellNLayers.setSize(mesh.nCells());
cellNLayers = 0;
faceRealThickness.setSize(mesh.nFaces());
faceRealThickness = 0;
// Mark all faces in the layer
const labelListList& layerFaces = addLayer.layerFaces();
@ -2415,27 +2443,195 @@ void Foam::autoLayerDriver::getLayerCellsFaces
{
const labelList& added = addedCells[oldPatchFaceI];
forAll(added, i)
const labelList& layer = layerFaces[oldPatchFaceI];
if (layer.size())
{
flaggedCells[added[i]] = true;
forAll(added, i)
{
cellNLayers[added[i]] = layer.size()-1;
}
}
}
forAll(layerFaces, oldPatchFaceI)
{
const labelList& layer = layerFaces[oldPatchFaceI];
const scalar realThickness = oldRealThickness[oldPatchFaceI];
if (layer.size())
{
for (label i = 1; i < layer.size()-1; i++)
// Layer contains both original boundary face and new boundary
// face so is nLayers+1
forAll(layer, i)
{
flaggedFaces[layer[i]] = true;
faceRealThickness[layer[i]] = realThickness;
}
}
}
}
bool Foam::autoLayerDriver::writeLayerData
(
const fvMesh& mesh,
const labelList& patchIDs,
const labelList& cellNLayers,
const scalarField& faceWantedThickness,
const scalarField& faceRealThickness
) const
{
bool allOk = true;
{
label nAdded = 0;
forAll(cellNLayers, cellI)
{
if (cellNLayers[cellI] > 0)
{
nAdded++;
}
}
cellSet addedCellSet(mesh, "addedCells", nAdded);
forAll(cellNLayers, cellI)
{
if (cellNLayers[cellI] > 0)
{
addedCellSet.insert(cellI);
}
}
addedCellSet.instance() = meshRefiner_.timeName();
Info<< "Writing "
<< returnReduce(addedCellSet.size(), sumOp<label>())
<< " added cells to cellSet "
<< addedCellSet.name() << endl;
bool ok = addedCellSet.write();
allOk = allOk & ok;
}
{
label nAdded = 0;
forAll(faceRealThickness, faceI)
{
if (faceRealThickness[faceI] > 0)
{
nAdded++;
}
}
faceSet layerFacesSet(mesh, "layerFaces", nAdded);
forAll(faceRealThickness, faceI)
{
if (faceRealThickness[faceI] > 0)
{
layerFacesSet.insert(faceI);
}
}
layerFacesSet.instance() = meshRefiner_.timeName();
Info<< "Writing "
<< returnReduce(layerFacesSet.size(), sumOp<label>())
<< " faces inside added layer to faceSet "
<< layerFacesSet.name() << endl;
bool ok = layerFacesSet.write();
allOk = allOk & ok;
}
{
volScalarField fld
(
IOobject
(
"nSurfaceLayers",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
mesh,
dimensionedScalar("zero", dimless, 0),
fixedValueFvPatchScalarField::typeName
);
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
forAll(patchIDs, i)
{
label patchI = patchIDs[i];
const polyPatch& pp = pbm[patchI];
const labelList& faceCells = pp.faceCells();
scalarField pfld(faceCells.size());
forAll(faceCells, i)
{
pfld[i] = cellNLayers[faceCells[i]];
}
fld.boundaryField()[patchI] == pfld;
}
Info<< "Writing volScalarField " << fld.name()
<< " with actual number of layers" << endl;
bool ok = fld.write();
allOk = allOk & ok;
}
{
volScalarField fld
(
IOobject
(
"wantedThickness",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
mesh,
dimensionedScalar("zero", dimless, 0),
fixedValueFvPatchScalarField::typeName
);
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
forAll(patchIDs, i)
{
label patchI = patchIDs[i];
fld.boundaryField()[patchI] == pbm[patchI].patchSlice
(
faceWantedThickness
);
}
Info<< "Writing volScalarField " << fld.name()
<< " with wanted thickness" << endl;
bool ok = fld.write();
allOk = allOk & ok;
}
{
volScalarField fld
(
IOobject
(
"thickness",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
mesh,
dimensionedScalar("zero", dimless, 0),
fixedValueFvPatchScalarField::typeName
);
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
forAll(patchIDs, i)
{
label patchI = patchIDs[i];
fld.boundaryField()[patchI] == pbm[patchI].patchSlice
(
faceRealThickness
);
}
Info<< "Writing volScalarField " << fld.name()
<< " with layer thickness" << endl;
bool ok = fld.write();
allOk = allOk & ok;
}
return allOk;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::autoLayerDriver::autoLayerDriver
@ -2618,7 +2814,7 @@ void Foam::autoLayerDriver::addLayers
labelList patchNLayers(pp().nPoints(), 0);
// Ideal number of cells added
label nIdealAddedCells = 0;
label nIdealTotAddedCells = 0;
// Whether to add edge for all pp.localPoints.
List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
@ -2637,7 +2833,7 @@ void Foam::autoLayerDriver::addLayers
patchDisp,
patchNLayers,
extrudeStatus,
nIdealAddedCells
nIdealTotAddedCells
);
// Precalculate mesh edge labels for patch edges
@ -2835,11 +3031,20 @@ void Foam::autoLayerDriver::addLayers
// Saved old points
pointField oldPoints(mesh.points());
// Last set of topology changes. (changing mesh clears out polyTopoChange)
// Current set of topology changes. (changing mesh clears out
// polyTopoChange)
polyTopoChange savedMeshMod(mesh.boundaryMesh().size());
// Per cell 0 or number of layers in the cell column it is part of
labelList cellNLayers;
// Per face actual overall layer thickness
scalarField faceRealThickness;
// Per face wanted overall layer thickness
scalarField faceWantedThickness(mesh.nFaces(), 0.0);
{
UIndirectList<scalar>(faceWantedThickness, pp().addressing()) =
avgPointData(pp, thickness);
}
boolList flaggedCells;
boolList flaggedFaces;
for (label iteration = 0; iteration < layerParams.nLayerIter(); iteration++)
{
@ -3088,35 +3293,52 @@ void Foam::autoLayerDriver::addLayers
(
newMesh,
addLayer,
flaggedCells,
flaggedFaces
avgPointData(pp, mag(patchDisp))(), // current thickness
cellNLayers,
faceRealThickness
);
// Count number of added cells
label nAddedCells = 0;
forAll(cellNLayers, cellI)
{
if (cellNLayers[cellI] > 0)
{
nAddedCells++;
}
}
if (debug&meshRefinement::MESH)
{
Info<< "Writing layer mesh to time " << meshRefiner_.timeName()
<< endl;
newMesh.write();
cellSet addedCellSet
(
newMesh,
"addedCells",
findIndices(flaggedCells, true)
);
cellSet addedCellSet(newMesh," addedCells", nAddedCells);
forAll(cellNLayers, cellI)
{
if (cellNLayers[cellI] > 0)
{
addedCellSet.insert(cellI);
}
}
addedCellSet.instance() = meshRefiner_.timeName();
Info<< "Writing "
<< returnReduce(addedCellSet.size(), sumOp<label>())
<< " added cells to cellSet "
<< addedCellSet.name() << endl;
<< " added cells to cellSet " << addedCellSet.name() << endl;
addedCellSet.write();
faceSet layerFacesSet
(
newMesh,
"layerFaces",
findIndices(flaggedCells, true)
);
faceSet layerFacesSet(newMesh, "layerFaces", newMesh.nFaces()/100);
forAll(faceRealThickness, faceI)
{
if (faceRealThickness[faceI] > 0)
{
layerFacesSet.insert(faceI);
}
}
layerFacesSet.instance() = meshRefiner_.timeName();
Info<< "Writing "
<< returnReduce(layerFacesSet.size(), sumOp<label>())
@ -3140,26 +3362,17 @@ void Foam::autoLayerDriver::addLayers
extrudeStatus
);
label nExtruded = countExtrusion(pp, extrudeStatus);
label nTotExtruded = countExtrusion(pp, extrudeStatus);
label nTotFaces = returnReduce(pp().size(), sumOp<label>());
label nAddedCells = 0;
{
forAll(flaggedCells, cellI)
{
if (flaggedCells[cellI])
{
nAddedCells++;
}
}
reduce(nAddedCells, sumOp<label>());
}
Info<< "Extruding " << nExtruded
label nTotAddedCells = returnReduce(nAddedCells, sumOp<label>());
Info<< "Extruding " << nTotExtruded
<< " out of " << nTotFaces
<< " faces (" << 100.0*nExtruded/nTotFaces << "%)."
<< " faces (" << 100.0*nTotExtruded/nTotFaces << "%)."
<< " Removed extrusion at " << nTotChanged << " faces."
<< endl
<< "Added " << nAddedCells << " out of " << nIdealAddedCells
<< " cells (" << 100.0*nAddedCells/nIdealAddedCells << "%)."
<< "Added " << nTotAddedCells << " out of " << nIdealTotAddedCells
<< " cells (" << 100.0*nTotAddedCells/nIdealTotAddedCells << "%)."
<< endl;
if (nTotChanged == 0)
@ -3217,6 +3430,8 @@ void Foam::autoLayerDriver::addLayers
meshRefiner_.updateMesh(map, labelList(0));
// Update numbering of faceWantedThickness
meshRefinement::updateList(map().faceMap(), 0.0, faceWantedThickness);
// Update numbering on baffles
forAll(baffles, i)
@ -3237,8 +3452,9 @@ void Foam::autoLayerDriver::addLayers
autoPtr<mapPolyMesh> map = meshRefiner_.mergeBaffles(baffles);
inplaceReorder(map().reverseCellMap(), flaggedCells);
inplaceReorder(map().reverseFaceMap(), flaggedFaces);
inplaceReorder(map().reverseCellMap(), cellNLayers);
inplaceReorder(map().reverseFaceMap(), faceWantedThickness);
inplaceReorder(map().reverseFaceMap(), faceRealThickness);
Info<< "Converted baffles in = "
<< meshRefiner_.mesh().time().cpuTimeIncrement()
@ -3272,29 +3488,23 @@ void Foam::autoLayerDriver::addLayers
);
// Re-distribute flag of layer faces and cells
map().distributeCellData(flaggedCells);
map().distributeFaceData(flaggedFaces);
map().distributeCellData(cellNLayers);
map().distributeFaceData(faceWantedThickness);
map().distributeFaceData(faceRealThickness);
}
// Write mesh
// ~~~~~~~~~~
// Write mesh data
// ~~~~~~~~~~~~~~~
cellSet addedCellSet(mesh, "addedCells", findIndices(flaggedCells, true));
addedCellSet.instance() = meshRefiner_.timeName();
Info<< "Writing "
<< returnReduce(addedCellSet.size(), sumOp<label>())
<< " added cells to cellSet "
<< addedCellSet.name() << endl;
addedCellSet.write();
faceSet layerFacesSet(mesh, "layerFaces", findIndices(flaggedFaces, true));
layerFacesSet.instance() = meshRefiner_.timeName();
Info<< "Writing "
<< returnReduce(layerFacesSet.size(), sumOp<label>())
<< " faces inside added layer to faceSet "
<< layerFacesSet.name() << endl;
layerFacesSet.write();
writeLayerData
(
mesh,
patchIDs,
cellNLayers,
faceWantedThickness,
faceRealThickness
);
}

24
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
View File

@ -120,6 +120,13 @@ class autoLayerDriver
const List<extrudeMode>&
);
//- Average point wise data to face wise
static tmp<scalarField> avgPointData
(
const indirectPrimitivePatch&,
const scalarField& pointFld
);
//- Check that primitivePatch is not multiply connected.
// Collect non-manifold points in pointSet.
static void checkManifold
@ -367,10 +374,23 @@ class autoLayerDriver
(
const polyMesh&,
const addPatchCellLayer&,
boolList&,
boolList&
const scalarField& oldRealThickness,
labelList& cellStatus,
scalarField& faceRealThickness
);
//- Write cellSet,faceSet for layers
bool writeLayerData
(
const fvMesh& mesh,
const labelList& patchIDs,
const labelList& cellNLayers,
const scalarField& faceWantedThickness,
const scalarField& faceRealThickness
) const;
// Mesh shrinking (to create space for layers)
//- Average field (over all subset of mesh points) by

37
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
View File

@ -235,6 +235,12 @@ Foam::layerParameters::layerParameters
Info<< "Layer thickness specified as final layer and expansion ratio."
<< endl;
}
else if (haveTotal && haveExp)
{
layerSpec_ = TOTAL_AND_EXPANSION;
Info<< "Layer thickness specified as overall thickness"
<< " and expansion ratio." << endl;
}
if (layerSpec_ == ILLEGAL || nSpec != 2)
@ -253,7 +259,9 @@ Foam::layerParameters::layerParameters
<< " final layer thickness ('finalLayerThickness')"
<< " and expansion ratio ('expansionRatio') or" << nl
<< " final layer thickness ('finalLayerThickness')"
<< " and overall thickness ('thickness')"
<< " and overall thickness ('thickness') or" << nl
<< " overall thickness ('thickness')"
<< " and expansion ratio ('expansionRatio'"
<< exit(FatalIOError);
}
@ -354,6 +362,19 @@ Foam::layerParameters::layerParameters
);
break;
case TOTAL_AND_EXPANSION:
layerDict.readIfPresent
(
"thickness",
thickness_[patchI]
);
layerDict.readIfPresent
(
"expansionRatio",
expansionRatio_[patchI]
);
break;
default:
FatalIOErrorIn
(
@ -390,6 +411,7 @@ Foam::scalar Foam::layerParameters::layerThickness
{
case FIRST_AND_TOTAL:
case FINAL_AND_TOTAL:
case TOTAL_AND_EXPANSION:
{
return totalThickness;
}
@ -450,6 +472,7 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
{
case FIRST_AND_EXPANSION:
case FINAL_AND_EXPANSION:
case TOTAL_AND_EXPANSION:
{
return expansionRatio;
}
@ -524,6 +547,18 @@ Foam::scalar Foam::layerParameters::firstLayerThickness
}
break;
case TOTAL_AND_EXPANSION:
{
scalar r = finalLayerThicknessRatio
(
nLayers,
expansionRatio
);
scalar finalThickness = r*totalThickness;
return finalThickness/pow(expansionRatio, nLayers-1);
}
break;
default:
{
FatalErrorIn("layerParameters::layerThickness(..)")

4
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
View File

@ -64,13 +64,15 @@ public:
// - first and expansion ratio specified
// - final and total thickness specified
// - final and expansion ratio specified
// - total thickness and expansion ratio specified
enum layerSpecification
{
ILLEGAL,
FIRST_AND_TOTAL,
FIRST_AND_EXPANSION,
FINAL_AND_TOTAL,
FINAL_AND_EXPANSION
FINAL_AND_EXPANSION,
TOTAL_AND_EXPANSION
};

11
src/parallel/reconstruct/reconstruct/processorMeshes.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -31,6 +31,15 @@ License
void Foam::processorMeshes::read()
{
forAll(databases_, procI)
{
meshes_.set(procI, NULL);
pointProcAddressing_.set(procI, NULL);
faceProcAddressing_.set(procI, NULL);
cellProcAddressing_.set(procI, NULL);
boundaryProcAddressing_.set(procI, NULL);
}
forAll(databases_, procI)
{
meshes_.set

28
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -53,7 +53,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
rhoName_("rho"),
lookupGravity_(-1),
g_(vector::zero),
relaxationFactor_(1)
curTimeIndex_(-1)
{}
@ -71,7 +71,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
rhoName_(dict.lookupOrDefault<word>("rhoName", "rho")),
lookupGravity_(-1),
g_(vector::zero),
relaxationFactor_(dict.lookupOrDefault<scalar>("relaxationFactor", 1))
curTimeIndex_(-1)
{
if (rhoName_ == "rhoInf")
{
@ -115,7 +115,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
rhoName_(ptf.rhoName_),
lookupGravity_(ptf.lookupGravity_),
g_(ptf.g_),
relaxationFactor_(ptf.relaxationFactor_)
curTimeIndex_(-1)
{}
@ -133,7 +133,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
rhoName_(ptf.rhoName_),
lookupGravity_(ptf.lookupGravity_),
g_(ptf.g_),
relaxationFactor_(ptf.relaxationFactor_)
curTimeIndex_(-1)
{}
@ -203,6 +203,13 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
const Time& t = mesh.time();
const pointPatch& ptPatch = this->patch();
// Store the motion state at the beginning of the time-step
if (curTimeIndex_ != this->db().time().timeIndex())
{
motion_.newTime();
curTimeIndex_ = this->db().time().timeIndex();
}
// Patch force data is valid for the current positions, so
// calculate the forces on the motion object from this data, then
// update the positions
@ -231,7 +238,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
g_ = g.value();
}
motion_.updateForce
motion_.updateAcceleration
(
f.forceEff() + g_*motion_.mass(),
f.momentEff(),
@ -240,11 +247,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
Field<vector>::operator=
(
relaxationFactor_
*(
motion_.currentPosition(initialPoints_)
- initialPoints_
)
motion_.currentPosition(initialPoints_) - initialPoints_
);
fixedValuePointPatchField<vector>::updateCoeffs();
@ -267,9 +270,6 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::write(Ostream& os) const
os.writeKeyword("g") << g_ << token::END_STATEMENT << nl;
}
os.writeKeyword("relaxationFactor")
<< relaxationFactor_ << token::END_STATEMENT << nl;
motion_.write(os);
initialPoints_.writeEntry("initialPoints", os);

6
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -84,8 +84,8 @@ class sixDoFRigidBodyDisplacementPointPatchVectorField
//- Gravity vector to store when not available from the db
vector g_;
//- Optional under-relaxation factor for the motion
scalar relaxationFactor_;
//- Current time index (used for updating)
label curTimeIndex_;
public:

147
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
View File

@ -164,6 +164,7 @@ void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
:
motionState_(),
motionState0_(),
restraints_(),
restraintNames_(),
constraints_(),
@ -173,8 +174,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
initialQ_(I),
momentOfInertia_(diagTensor::one*VSMALL),
mass_(VSMALL),
cDamp_(0.0),
aLim_(VGREAT),
aRelax_(1.0),
report_(false)
{}
@ -191,8 +191,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
const point& initialCentreOfMass,
const tensor& initialQ,
const diagTensor& momentOfInertia,
scalar cDamp,
scalar aLim,
scalar aRelax,
bool report
)
:
@ -205,6 +204,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
pi,
tau
),
motionState0_(motionState_),
restraints_(),
restraintNames_(),
constraints_(),
@ -214,8 +214,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
initialQ_(initialQ),
momentOfInertia_(momentOfInertia),
mass_(mass),
cDamp_(cDamp),
aLim_(aLim),
aRelax_(aRelax),
report_(report)
{}
@ -223,6 +222,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
:
motionState_(dict),
motionState0_(motionState_),
restraints_(),
restraintNames_(),
constraints_(),
@ -238,8 +238,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
),
momentOfInertia_(dict.lookup("momentOfInertia")),
mass_(readScalar(dict.lookup("mass"))),
cDamp_(dict.lookupOrDefault<scalar>("accelerationDampingCoeff", 0.0)),
aLim_(dict.lookupOrDefault<scalar>("accelerationLimit", VGREAT)),
aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
report_(dict.lookupOrDefault<Switch>("report", false))
{
addRestraints(dict);
@ -254,6 +253,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
)
:
motionState_(sDoFRBM.motionState_),
motionState0_(sDoFRBM.motionState0_),
restraints_(sDoFRBM.restraints_),
restraintNames_(sDoFRBM.restraintNames_),
constraints_(sDoFRBM.constraints_),
@ -263,8 +263,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
initialQ_(sDoFRBM.initialQ_),
momentOfInertia_(sDoFRBM.momentOfInertia_),
mass_(sDoFRBM.mass_),
cDamp_(sDoFRBM.cDamp_),
aLim_(sDoFRBM.aLim_),
aRelax_(sDoFRBM.aRelax_),
report_(sDoFRBM.report_)
{}
@ -376,91 +375,59 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
if (Pstream::master())
{
vector aClip = a();
scalar aMag = mag(aClip);
if (aMag > SMALL)
{
aClip /= aMag;
}
if (aMag > aLim_)
{
WarningIn
(
"void Foam::sixDoFRigidBodyMotion::updatePosition"
"("
"scalar, "
"scalar"
")"
)
<< "Limited acceleration " << a()
<< " to " << aClip*aLim_
<< endl;
}
v() += 0.5*(1 - cDamp_)*deltaT0*aClip*min(aMag, aLim_);
pi() += 0.5*(1 - cDamp_)*deltaT0*tau();
v() = v0() + 0.5*deltaT0*a();
pi() = pi0() + 0.5*deltaT0*tau();
// Leapfrog move part
centreOfMass() += deltaT*v();
centreOfMass() = centreOfMass0() + deltaT*v();
// Leapfrog orientation adjustment
rotate(Q(), pi(), deltaT);
Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
Q() = Qpi.first();
pi() = Qpi.second();
}
Pstream::scatter(motionState_);
}
void Foam::sixDoFRigidBodyMotion::updateForce
void Foam::sixDoFRigidBodyMotion::updateAcceleration
(
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT
)
{
static bool first = true;
// Second leapfrog velocity adjust part, required after motion and
// force calculation
// acceleration calculation
if (Pstream::master())
{
// Save the previous iteration accelerations for relaxation
vector aPrevIter = a();
vector tauPrevIter = tau();
// Calculate new accelerations
a() = fGlobal/mass_;
tau() = (Q().T() & tauGlobal);
applyRestraints();
// Relax accelerations on all but first iteration
if (!first)
{
a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
}
first = false;
// Apply constraints after relaxation
applyConstraints(deltaT);
vector aClip = a();
scalar aMag = mag(aClip);
if (aMag > SMALL)
{
aClip /= aMag;
}
if (aMag > aLim_)
{
WarningIn
(
"void Foam::sixDoFRigidBodyMotion::updateForce"
"("
"const vector&, "
"const vector&, "
"scalar"
")"
)
<< "Limited acceleration " << a()
<< " to " << aClip*aLim_
<< endl;
}
v() += 0.5*(1 - cDamp_)*deltaT*aClip*min(aMag, aLim_);
pi() += 0.5*(1 - cDamp_)*deltaT*tau();
// Correct velocities
v() += 0.5*deltaT*a();
pi() += 0.5*deltaT*tau();
if (report_)
{
@ -472,7 +439,27 @@ void Foam::sixDoFRigidBodyMotion::updateForce
}
void Foam::sixDoFRigidBodyMotion::updateForce
void Foam::sixDoFRigidBodyMotion::updateVelocity(scalar deltaT)
{
// Second leapfrog velocity adjust part, required after motion and
// acceleration calculation
if (Pstream::master())
{
v() += 0.5*deltaT*a();
pi() += 0.5*deltaT*tau();
if (report_)
{
status();
}
}
Pstream::scatter(motionState_);
}
void Foam::sixDoFRigidBodyMotion::updateAcceleration
(
const pointField& positions,
const vectorField& forces,
@ -493,7 +480,7 @@ void Foam::sixDoFRigidBodyMotion::updateForce
}
}
updateForce(fGlobal, tauGlobal, deltaT);
updateAcceleration(fGlobal, tauGlobal, deltaT);
}
@ -506,19 +493,14 @@ Foam::point Foam::sixDoFRigidBodyMotion::predictedPosition
) const
{
vector vTemp = v() + deltaT*(a() + deltaForce/mass_);
vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
point centreOfMassTemp = centreOfMass() + deltaT*vTemp;
tensor QTemp = Q();
rotate(QTemp, piTemp, deltaT);
point centreOfMassTemp = centreOfMass0() + deltaT*vTemp;
Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
return
(
centreOfMassTemp
+ (QTemp & initialQ_.T() & (pInitial - initialCentreOfMass_))
+ (QpiTemp.first() & initialQ_.T() & (pInitial - initialCentreOfMass_))
);
}
@ -531,12 +513,9 @@ Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation
) const
{
vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
tensor QTemp = Q();
rotate(QTemp, piTemp, deltaT);
return (QTemp & initialQ_.T() & vInitial);
return (QpiTemp.first() & initialQ_.T() & vInitial);
}

64
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -59,7 +59,7 @@ SourceFiles
#include "pointField.H"
#include "sixDoFRigidBodyMotionRestraint.H"
#include "sixDoFRigidBodyMotionConstraint.H"
#include "Tuple2.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -87,6 +87,9 @@ class sixDoFRigidBodyMotion
//- Motion state data object
sixDoFRigidBodyMotionState motionState_;
//- Motion state data object for previous time-step
sixDoFRigidBodyMotionState motionState0_;
//- Restraints on the motion
PtrList<sixDoFRigidBodyMotionRestraint> restraints_;
@ -116,14 +119,8 @@ class sixDoFRigidBodyMotion
//- Mass of the body
scalar mass_;
//- Acceleration damping coefficient. Modify applied acceleration:
// v1 = v0 + a*dt - cDamp*a*dt
// = v0 + dt*f*(1 - cDamp)/m
// Increases effective mass by 1/(1 - cDamp).
scalar cDamp_;
//- Acceleration magnitude limit - clips large accelerations
scalar aLim_;
//- Acceleration relaxation coefficient
scalar aRelax_;
//- Switch to turn reporting of motion data on and off
Switch report_;
@ -143,8 +140,14 @@ class sixDoFRigidBodyMotion
// frame z-axis by the given angle
inline tensor rotationTensorZ(scalar deltaT) const;
//- Apply rotation tensors to Q for the given torque (pi) and deltaT
inline void rotate(tensor& Q, vector& pi, scalar deltaT) const;
//- Apply rotation tensors to Q0 for the given torque (pi) and deltaT
// and return the rotated Q and pi as a tuple
inline Tuple2<tensor, vector> rotate
(
const tensor& Q0,
const vector& pi,
const scalar deltaT
) const;
//- Apply the restraints to the object
void applyRestraints();
@ -152,6 +155,7 @@ class sixDoFRigidBodyMotion
//- Apply the constraints to the object
void applyConstraints(scalar deltaT);
// Access functions retained as private because of the risk of
// confusion over what is a body local frame vector and what is global
@ -199,6 +203,21 @@ class sixDoFRigidBodyMotion
//- Return access to torque
inline const vector& tau() const;
//- Return access to the orientation at previous time-step
inline const tensor& Q0() const;
//- Return access to velocity at previous time-step
inline const vector& v0() const;
//- Return access to acceleration at previous time-step
inline const vector& a0() const;
//- Return access to angular momentum at previous time-step
inline const vector& pi0() const;
//- Return access to torque at previous time-step
inline const vector& tau0() const;
// Edit
@ -244,8 +263,7 @@ public:
const point& initialCentreOfMass,
const tensor& initialQ,
const diagTensor& momentOfInertia,
scalar cDamp = 0.0,
scalar aLim = VGREAT,
scalar aRelax = 1.0,
bool report = false
);
@ -279,18 +297,22 @@ public:
scalar deltaT0
);
//- Second leapfrog velocity adjust part, required after motion and
// force calculation
void updateForce
//- Second leapfrog velocity adjust part
// required after motion and force calculation
void updateAcceleration
(
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT
);
//- Second leapfrog velocity adjust part
// required after motion and force calculation
void updateVelocity(scalar deltaT);
//- Global forces supplied at locations, calculating net force
// and moment
void updateForce
void updateAcceleration
(
const pointField& positions,
const vectorField& forces,
@ -354,6 +376,9 @@ public:
//- Return const access to the centre of mass
inline const point& centreOfMass() const;
//- Return const access to the centre of mass at previous time-step
inline const point& centreOfMass0() const;
//- Return access to the inertia tensor
inline const diagTensor& momentOfInertia() const;
@ -366,6 +391,9 @@ public:
// Edit
//- Store the motion state at the beginning of the time-step
inline void newTime();
//- Return non-const access to the centre of mass
inline point& centreOfMass();

61
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -61,16 +61,19 @@ Foam::sixDoFRigidBodyMotion::rotationTensorZ(scalar phi) const
}
inline void Foam::sixDoFRigidBodyMotion::rotate
inline Foam::Tuple2<Foam::tensor, Foam::vector>
Foam::sixDoFRigidBodyMotion::rotate
(
tensor& Q,
vector& pi,
scalar deltaT
const tensor& Q0,
const vector& pi0,
const scalar deltaT
) const
{
tensor R;
Tuple2<tensor, vector> Qpi(Q0, pi0);
tensor& Q = Qpi.first();
vector& pi = Qpi.second();
R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
tensor R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
pi = pi & R;
Q = Q & R;
@ -89,6 +92,8 @@ inline void Foam::sixDoFRigidBodyMotion::rotate
R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
pi = pi & R;
Q = Q & R;
return Qpi;
}
@ -176,6 +181,36 @@ inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau() const
}
inline const Foam::tensor& Foam::sixDoFRigidBodyMotion::Q0() const
{
return motionState0_.Q();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::v0() const
{
return motionState0_.v();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::a0() const
{
return motionState0_.a();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::pi0() const
{
return motionState0_.pi();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau0() const
{
return motionState0_.tau();
}
inline Foam::point& Foam::sixDoFRigidBodyMotion::initialCentreOfMass()
{
return initialCentreOfMass_;
@ -281,6 +316,12 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const
}
inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass0() const
{
return motionState0_.centreOfMass();
}
inline const Foam::diagTensor&
Foam::sixDoFRigidBodyMotion::momentOfInertia() const
{
@ -300,6 +341,12 @@ inline bool Foam::sixDoFRigidBodyMotion::report() const
}
inline void Foam::sixDoFRigidBodyMotion::newTime()
{
motionState0_ = motionState_;
}
inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass()
{
return motionState_.centreOfMass();

8
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -40,10 +40,8 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
<< momentOfInertia_ << token::END_STATEMENT << nl;
os.writeKeyword("mass")
<< mass_ << token::END_STATEMENT << nl;
os.writeKeyword("accelerationDampingCoeff")
<< cDamp_ << token::END_STATEMENT << nl;
os.writeKeyword("accelerationLimit")
<< aLim_ << token::END_STATEMENT << nl;
os.writeKeyword("accelerationRelaxation")
<< aRelax_ << token::END_STATEMENT << nl;
os.writeKeyword("report")
<< report_ << token::END_STATEMENT << nl;

8
src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -160,7 +160,11 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
}
// Do not modify the accelerations, except with gravity, where available
motion_.updateForce(gravity*motion_.mass(), vector::zero, t.deltaTValue());
motion_.updateAcceleration
(
gravity*motion_.mass(),
vector::zero, t.deltaTValue()
);
Field<vector>::operator=
(

2
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
View File

@ -36,7 +36,7 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
)
:
CompType(mesh),
ODE(),
ODESystem(),
Y_(this->thermo().composition().Y()),
reactions_
(

4
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
View File

@ -39,7 +39,7 @@ SourceFiles
#define chemistryModel_H
#include "Reaction.H"
#include "ODE.H"
#include "ODESystem.H"
#include "volFieldsFwd.H"
#include "simpleMatrix.H"
#include "DimensionedField.H"
@ -60,7 +60,7 @@ template<class CompType, class ThermoType>
class chemistryModel
:
public CompType,
public ODE
public ODESystem
{
// Private Member Functions

2
src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
View File

@ -37,7 +37,7 @@ solidChemistryModel
)
:
CompType(mesh),
ODE(),
ODESystem(),
Ys_(this->solidThermo().composition().Y()),
reactions_
(

4
src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
View File

@ -40,7 +40,7 @@ SourceFiles
#define solidChemistryModel_H
#include "Reaction.H"
#include "ODE.H"
#include "ODESystem.H"
#include "volFieldsFwd.H"
#include "DimensionedField.H"
#include "simpleMatrix.H"
@ -61,7 +61,7 @@ template<class CompType, class SolidThermo>
class solidChemistryModel
:
public CompType,
public ODE
public ODESystem
{
// Private Member Functions

10
tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict
View File

@ -17,7 +17,7 @@ FoamFile
application interDyMFoam;
startFrom latestTime;
startFrom startTime;
startTime 0;
@ -29,7 +29,7 @@ deltaT 0.01;
writeControl adjustableRunTime;
writeInterval 0.1;
writeInterval 0.2;
purgeWrite 0;
@ -47,9 +47,9 @@ runTimeModifiable yes;
adjustTimeStep yes;
maxCo 0.5;
maxAlphaCo 0.5;
maxDeltaT 0.01;
maxCo 5;
maxAlphaCo 5;
maxDeltaT 1;
libs
(

4
tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
View File

@ -27,7 +27,7 @@ gradSchemes
divSchemes
{
div(rho*phi,U) Gauss vanLeerV;
div(rhoPhi,U) Gauss vanLeerV;
div(phi,alpha) Gauss vanLeer;
div(phirb,alpha) Gauss interfaceCompression;
div(phi,k) Gauss upwind;
@ -55,7 +55,7 @@ fluxRequired
default no;
p_rgh;
pcorr;
alpha;
alpha.water;
}

24
tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
View File

@ -17,7 +17,7 @@ FoamFile
solvers
{
cellDisplacement
"cellDisplacement.*"
{
solver GAMG;
tolerance 1e-5;
@ -29,14 +29,24 @@ solvers
mergeLevels 1;
}
alpha.water
"alpha.water.*"
{
nAlphaCorr 1;
nAlphaSubCycles 3;
nAlphaSubCycles 1;
cAlpha 1;
alphaOuterCorrectors yes;
MULESCorr yes;
nLimiterIter 5;
solver PBiCG;
preconditioner DILU;
tolerance 1e-8;
relTol 0;
}
pcorr
"pcorr.*"
{
solver PCG;
preconditioner
@ -120,9 +130,11 @@ solvers
PIMPLE
{
momentumPredictor no;
nCorrectors 2;
nOuterCorrectors 6;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
correctPhi yes;
moveMeshOuterCorrectors yes;
}
relaxationFactors
@ -132,7 +144,7 @@ relaxationFactors
}
equations
{
"(U|k|epsilon|omega|R|nuTilda).*" 1;
".*" 1;
}
}