zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2010-02-04 11:51:50 +00:00
parent 76c0a7bfb9 10d4c94d82
commit 8505ac08dd
93 changed files with 6966 additions and 1284 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -434,6 +434,19 @@ int main(int argc, char *argv[])
// Layer addition parameters
layerParameters layerParams(layerDict, mesh.boundaryMesh());
//!!! Temporary hack to get access to maxLocalCells
bool preBalance;
{
refinementParameters refineParams(refineDict);
preBalance = returnReduce
(
(mesh.nCells() >= refineParams.maxLocalCells()),
orOp<bool>()
);
}
if (!overwrite)
{
const_cast<Time&>(mesh.time())++;
@ -444,6 +457,7 @@ int main(int argc, char *argv[])
layerDict,
motionDict,
layerParams,
preBalance,
decomposer,
distributor
);

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

@ -10,7 +10,7 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
object autoHexMeshDict;
object snappyHexMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -40,7 +40,7 @@ geometry
{
type triSurfaceMesh;
//tolerance 1E-6; // optional:non-default tolerance on intersections
//tolerance 1E-5; // optional:non-default tolerance on intersections
//maxTreeDepth 10; // optional:depth of octree. Decrease only in case
// of memory limitations.
@ -71,9 +71,9 @@ castellatedMeshControls
// Refinement parameters
// ~~~~~~~~~~~~~~~~~~~~~
// While refining maximum number of cells per processor. This is basically
// the number of cells that fit on a processor. If you choose this too small
// it will do just more refinement iterations to obtain a similar mesh.
// If local number of cells is >= maxLocalCells on any processor
// switches from from refinement followed by balancing
// (current method) to (weighted) balancing before refinement.
maxLocalCells 1000000;
// Overall cell limit (approximately). Refinement will stop immediately
@ -89,6 +89,13 @@ castellatedMeshControls
// (unless the number of cells to refine is 0)
minRefinementCells 0;
// Allow a certain level of imbalance during refining
// (since balancing is quite expensive)
// Expressed as fraction of perfect balance (= overall number of cells /
// nProcs). 0=balance always.
maxLoadUnbalance 0.10;
// Number of buffer layers between different levels.
// 1 means normal 2:1 refinement restriction, larger means slower
// refinement.

248
applications/utilities/mesh/manipulation/setSet/setSet.C
View File

@ -45,6 +45,7 @@ Description
#include "cellZoneSet.H"
#include "faceZoneSet.H"
#include "pointZoneSet.H"
#include "timeSelector.H"
#include <stdio.h>
@ -82,45 +83,6 @@ Istream& selectStream(Istream* is0Ptr, Istream* is1Ptr)
}
}
// Copy set
void backup
(
const word& setType,
const polyMesh& mesh,
const word& fromName,
const topoSet& fromSet,
const word& toName
)
{
if (fromSet.size())
{
Info<< " Backing up " << fromName << " into " << toName << endl;
topoSet::New(setType, mesh, toName, fromSet)().write();
}
}
// Read and copy set
void backup
(
const word& setType,
const polyMesh& mesh,
const word& fromName,
const word& toName
)
{
autoPtr<topoSet> fromSet = topoSet::New
(
setType,
mesh,
fromName,
IOobject::READ_IF_PRESENT
);
backup(setType, mesh, fromName, fromSet(), toName);
}
// Write set to VTK readable files
void writeVTK
@ -304,7 +266,13 @@ void printAllSets(const polyMesh& mesh, Ostream& os)
IOobjectList objects
(
mesh,
mesh.pointsInstance(),
mesh.time().findInstance
(
polyMesh::meshSubDir/"sets",
word::null,
IOobject::READ_IF_PRESENT,
mesh.facesInstance()
),
polyMesh::meshSubDir/"sets"
);
IOobjectList cellSets(objects.lookupClass(cellSet::typeName));
@ -417,7 +385,13 @@ void removeSet
IOobjectList objects
(
mesh,
mesh.pointsInstance(),
mesh.time().findInstance
(
polyMesh::meshSubDir/"sets",
word::null,
IOobject::READ_IF_PRESENT,
mesh.facesInstance()
),
polyMesh::meshSubDir/"sets"
);
@ -465,6 +439,7 @@ bool doCommand
const word& setName,
const word& actionName,
const bool writeVTKFile,
const bool writeCurrentTime,
Istream& is
)
{
@ -628,8 +603,6 @@ bool doCommand
/currentSet.name()
<< " and to vtk file " << vtkName << endl << endl;
currentSet.write();
writeVTK(mesh, currentSet, vtkName);
}
else
@ -638,9 +611,13 @@ bool doCommand
<< " (size " << currentSet.size() << ") to "
<< currentSet.instance()/currentSet.local()
/currentSet.name() << endl << endl;
currentSet.write();
}
if (writeCurrentTime)
{
currentSet.instance() = mesh.time().timeName();
}
currentSet.write();
}
}
}
@ -692,6 +669,48 @@ void printMesh(const Time& runTime, const polyMesh& mesh)
}
polyMesh::readUpdateState meshReadUpdate(polyMesh& mesh)
{
polyMesh::readUpdateState stat = mesh.readUpdate();
switch(stat)
{
case polyMesh::UNCHANGED:
{
Info<< " mesh not changed." << endl;
break;
}
case polyMesh::POINTS_MOVED:
{
Info<< " points moved; topology unchanged." << endl;
break;
}
case polyMesh::TOPO_CHANGE:
{
Info<< " topology changed; patches unchanged." << nl
<< " ";
printMesh(mesh.time(), mesh);
break;
}
case polyMesh::TOPO_PATCH_CHANGE:
{
Info<< " topology changed and patches changed." << nl
<< " ";
printMesh(mesh.time(), mesh);
break;
}
default:
{
FatalErrorIn("meshReadUpdate(polyMesh&)")
<< "Illegal mesh update state "
<< stat << abort(FatalError);
break;
}
}
return stat;
}
commandStatus parseType
(
@ -730,44 +749,10 @@ commandStatus parseType
<< " to " << Times[nearestIndex].name()
<< endl;
// Set time
runTime.setTime(Times[nearestIndex], nearestIndex);
polyMesh::readUpdateState stat = mesh.readUpdate();
switch(stat)
{
case polyMesh::UNCHANGED:
{
Info<< " mesh not changed." << endl;
break;
}
case polyMesh::POINTS_MOVED:
{
Info<< " points moved; topology unchanged." << endl;
break;
}
case polyMesh::TOPO_CHANGE:
{
Info<< " topology changed; patches unchanged." << nl
<< " ";
printMesh(runTime, mesh);
break;
}
case polyMesh::TOPO_PATCH_CHANGE:
{
Info<< " topology changed and patches changed." << nl
<< " ";
printMesh(runTime, mesh);
break;
}
default:
{
FatalErrorIn("parseType")
<< "Illegal mesh update state "
<< stat << abort(FatalError);
break;
}
}
// Optionally re-read mesh
meshReadUpdate(mesh);
return INVALID;
}
@ -840,23 +825,28 @@ commandStatus parseAction(const word& actionName)
int main(int argc, char *argv[])
{
timeSelector::addOptions(true, false);
# include "addRegionOption.H"
# include "addTimeOptions.H"
argList::addBoolOption("noVTK");
argList::addBoolOption("noVTK", "do not write VTK files");
argList::addBoolOption("loop", "execute batch commands for all timesteps");
argList::addOption("batch", "file");
# include "setRootCase.H"
# include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
bool writeVTK = !args.optionFound("noVTK");
bool loop = args.optionFound("loop");
bool batch = args.optionFound("batch");
// Get times list
instantList Times = runTime.times();
# include "checkTimeOptions.H"
if (loop && !batch)
{
FatalErrorIn(args.executable())
<< "Can only loop when in batch mode."
<< exit(FatalError);
}
runTime.setTime(Times[startTime], startTime);
# include "createNamedPolyMesh.H"
@ -866,11 +856,34 @@ int main(int argc, char *argv[])
// Print current sets
printAllSets(mesh, Info);
// Read history if interactive
# if READLINE != 0
if (!batch && !read_history(historyFile))
{
Info<< "Successfully read history from " << historyFile << endl;
}
# endif
std::ifstream* fileStreamPtr(NULL);
// Exit status
int status = 0;
if (args.optionFound("batch"))
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Time = " << runTime.timeName() << endl;
// Handle geometry/topology changes
meshReadUpdate(mesh);
// Main command read & execute loop
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
autoPtr<IFstream> fileStreamPtr(NULL);
if (batch)
{
fileName batchFile(args.option("batch"));
@ -883,18 +896,14 @@ int main(int argc, char *argv[])
<< "Cannot open file " << batchFile << exit(FatalError);
}
fileStreamPtr = new std::ifstream(batchFile.c_str());
fileStreamPtr.reset(new IFstream(batchFile));
}
#if READLINE != 0
else if (!read_history(historyFile))
{
Info<< "Successfully read history from " << historyFile << endl;
}
#endif
Info<< "Please type 'help', 'quit' or a set command after prompt." << endl;
Info<< "Please type 'help', 'quit' or a set command after prompt."
<< endl;
bool ok = true;
// Whether to quit
bool quit = false;
FatalError.throwExceptions();
FatalIOError.throwExceptions();
@ -912,15 +921,16 @@ int main(int argc, char *argv[])
commandStatus stat = INVALID;
if (fileStreamPtr)
if (fileStreamPtr.valid())
{
if (!fileStreamPtr->good())
if (!fileStreamPtr().good())
{
Info<< "End of batch file" << endl;
// No error.
break;
}
std::getline(*fileStreamPtr, rawLine);
fileStreamPtr().getLine(rawLine);
if (rawLine.size())
{
@ -973,28 +983,44 @@ int main(int argc, char *argv[])
}
}
}
ok = true;
if (stat == QUIT)
{
break;
// Make sure to quit
quit = true;
}
else if (stat == VALIDSETCMD || stat == VALIDZONECMD)
{
ok = doCommand(mesh, setType, setName, actionName, writeVTK, is);
bool ok = doCommand
(
mesh,
setType,
setName,
actionName,
writeVTK,
loop, // if in looping mode dump sets to time directory
is
);
if (!ok)
{
// Exit with error.
quit = true;
status = 1;
}
}
} while (ok);
} while (!quit);
if (fileStreamPtr)
if (quit)
{
delete fileStreamPtr;
break;
}
}
Info<< "\nEnd\n" << endl;
return 0;
return status;
}

13
applications/utilities/mesh/manipulation/setsToZones/setsToZones.C
View File

@ -81,14 +81,17 @@ int main(int argc, char *argv[])
# include "createNamedPolyMesh.H"
// Search for list of objects for the time of the mesh
IOobjectList objects
word setsInstance = runTime.findInstance
(
mesh,
mesh.pointsInstance(),
polyMesh::meshSubDir/"sets"
polyMesh::meshSubDir/"sets",
word::null,
IOobject::MUST_READ,
mesh.facesInstance()
);
Info<< "Searched : " << mesh.pointsInstance()/polyMesh::meshSubDir/"sets"
IOobjectList objects(mesh, setsInstance, polyMesh::meshSubDir/"sets");
Info<< "Searched : " << setsInstance/polyMesh::meshSubDir/"sets"
<< nl
<< "Found : " << objects.names() << nl
<< endl;

33
applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C
View File

@ -400,25 +400,30 @@ void Foam::vtkPV3Foam::updateInfoSets
Info<< "<beg> Foam::vtkPV3Foam::updateInfoSets" << endl;
}
// Add names of sets
IOobjectList objects
(
dbPtr_(),
dbPtr_().findInstance(meshDir_, "faces", IOobject::READ_IF_PRESENT),
meshDir_/"sets"
);
// Add names of sets. Search for last time directory with a sets
// subdirectory. Take care not to search beyond the last mesh.
if (debug)
{
Info<< " Foam::vtkPV3Foam::updateInfoSets read "
<< objects.names() << " from "
<< dbPtr_().findInstance
word facesInstance = dbPtr_().findInstance
(
meshDir_,
"faces",
IOobject::READ_IF_PRESENT
)
<< endl;
);
word setsInstance = dbPtr_().findInstance
(
meshDir_/"sets",
word::null,
IOobject::READ_IF_PRESENT,
facesInstance
);
IOobjectList objects(dbPtr_(), setsInstance, meshDir_/"sets");
if (debug)
{
Info<< " Foam::vtkPV3Foam::updateInfoSets read "
<< objects.names() << " from " << setsInstance << endl;
}

26
applications/utilities/postProcessing/sampling/sample/sampleDict
View File

@ -118,9 +118,7 @@ sets
);
// Surface sampling definition: choice of
// plane : values on plane defined by point, normal.
// patch : values on patch.
// Surface sampling definition
//
// 1] patches are not triangulated by default
// 2] planes are always triangulated
@ -209,6 +207,28 @@ surfaces
// regularise false; // Optional: do not simplify
}
distance
{
// Isosurface from signed/unsigned distance to surface
type distanceSurface;
signed true;
// Definition of surface
surfaceType triSurfaceMesh;
surfaceName integrationPlane.stl;
// Distance to surface
distance 0.0;
interpolate false;
}
triSurfaceSampling
{
// Sampling on triSurface
type sampledTriSurfaceMesh;
surface integrationPlane.stl;
interpolate true;
}
);

2
applications/utilities/surface/surfaceInertia/surfaceInertia.C
View File

@ -296,7 +296,7 @@ int main(int argc, char *argv[])
(
"density",
"scalar",
"Specify density,"
"Specify density, "
"kg/m3 for solid properties, kg/m2 for shell properties"
);

6
src/OpenFOAM/db/Time/Time.H
View File

@ -283,12 +283,14 @@ public:
//- Return the location of "dir" containing the file "name".
// (eg, used in reading mesh data)
// If name is null, search for the directory "dir" only
// If name is null, search for the directory "dir" only.
// Does not search beyond stopInstance (if set) or constant.
word findInstance
(
const fileName& dir,
const word& name = word::null,
const IOobject::readOption rOpt = IOobject::MUST_READ
const IOobject::readOption rOpt = IOobject::MUST_READ,
const word& stopInstance = word::null
) const;
//- Search the case for valid time directories

62
src/OpenFOAM/db/Time/findInstance.C
View File

@ -39,19 +39,24 @@ Foam::word Foam::Time::findInstance
(
const fileName& dir,
const word& name,
const IOobject::readOption rOpt
const IOobject::readOption rOpt,
const word& stopInstance
) const
{
// Note: if name is empty, just check the directory itself
const fileName tPath(path());
const fileName dirPath(tPath/timeName()/dir);
// check the current time directory
if
(
name.empty()
? isDir(path()/timeName()/dir)
? isDir(dirPath)
:
(
isFile(path()/timeName()/dir/name)
isFile(dirPath/name)
&& IOobject(name, timeName(), dir, *this).headerOk()
)
)
@ -59,7 +64,8 @@ Foam::word Foam::Time::findInstance
if (debug)
{
Info<< "Time::findInstance"
"(const fileName&, const word&, const IOobject::readOption)"
"(const fileName&, const word&"
", const IOobject::readOption, const word&)"
<< " : found \"" << name
<< "\" in " << timeName()/dir
<< endl;
@ -88,10 +94,10 @@ Foam::word Foam::Time::findInstance
if
(
name.empty()
? isDir(path()/ts[instanceI].name()/dir)
? isDir(tPath/ts[instanceI].name()/dir)
:
(
isFile(path()/ts[instanceI].name()/dir/name)
isFile(tPath/ts[instanceI].name()/dir/name)
&& IOobject(name, ts[instanceI].name(), dir, *this).headerOk()
)
)
@ -99,7 +105,8 @@ Foam::word Foam::Time::findInstance
if (debug)
{
Info<< "Time::findInstance"
"(const fileName&, const word&, const IOobject::readOption)"
"(const fileName&, const word&"
", const IOobject::readOption, const word&)"
<< " : found \"" << name
<< "\" in " << ts[instanceI].name()/dir
<< endl;
@ -107,6 +114,34 @@ Foam::word Foam::Time::findInstance
return ts[instanceI].name();
}
// Check if hit minimum instance
if (ts[instanceI].name() == stopInstance)
{
if (debug)
{
Info<< "Time::findInstance"
"(const fileName&, const word&"
", const IOobject::readOption, const word&)"
<< " : hit stopInstance " << stopInstance
<< endl;
}
if (rOpt == IOobject::MUST_READ)
{
FatalErrorIn
(
"Time::findInstance"
"(const fileName&, const word&"
", const IOobject::readOption, const word&)"
) << "Cannot find file \"" << name << "\" in directory "
<< dir << " in times " << timeName()
<< " down to " << stopInstance
<< exit(FatalError);
}
return ts[instanceI].name();
}
}
@ -119,10 +154,10 @@ Foam::word Foam::Time::findInstance
if
(
name.empty()
? isDir(path()/constant()/dir)
? isDir(tPath/constant()/dir)
:
(
isFile(path()/constant()/dir/name)
isFile(tPath/constant()/dir/name)
&& IOobject(name, constant(), dir, *this).headerOk()
)
)
@ -130,7 +165,8 @@ Foam::word Foam::Time::findInstance
if (debug)
{
Info<< "Time::findInstance"
"(const fileName&, const word&, const IOobject::readOption)"
"(const fileName&, const word&"
", const IOobject::readOption, const word&)"
<< " : found \"" << name
<< "\" in " << constant()/dir
<< endl;
@ -144,9 +180,11 @@ Foam::word Foam::Time::findInstance
FatalErrorIn
(
"Time::findInstance"
"(const fileName&, const word&, const IOobject::readOption)"
"(const fileName&, const word&"
", const IOobject::readOption, const word&)"
) << "Cannot find file \"" << name << "\" in directory "
<< constant()/dir
<< dir << " in times " << timeName()
<< " down to " << constant()
<< exit(FatalError);
}

24
src/OpenFOAM/dimensionSet/dimensionSets.C
View File

@ -38,19 +38,21 @@ const Foam::dimensionSet Foam::dimMoles(0, 0, 0, 0, 1, 0, 0);
const Foam::dimensionSet Foam::dimCurrent(0, 0, 0, 0, 0, 1, 0);
const Foam::dimensionSet Foam::dimLuminousIntensity(0, 0, 0, 0, 0, 0, 1);
const Foam::dimensionSet Foam::dimArea(0, 2, 0, 0, 0, 0, 0);
const Foam::dimensionSet Foam::dimVolume(0, 3, 0, 0, 0, 0, 0);
const Foam::dimensionSet Foam::dimVol(0, 3, 0, 0, 0, 0, 0);
const Foam::dimensionSet Foam::dimArea(sqr(dimLength));
const Foam::dimensionSet Foam::dimVolume(pow3(dimLength));
const Foam::dimensionSet Foam::dimVol(dimVolume);
const Foam::dimensionSet Foam::dimDensity(1, -3, 0, 0, 0, 0, 0);
const Foam::dimensionSet Foam::dimForce(1, 1, -2, 0, 0, 0, 0);
const Foam::dimensionSet Foam::dimEnergy(1, 2, -2, 0, 0, 0, 0);
const Foam::dimensionSet Foam::dimVelocity(dimLength/dimTime);
const Foam::dimensionSet Foam::dimAcceleration(dimVelocity/dimTime);
const Foam::dimensionSet Foam::dimVelocity(0, 1, -1, 0, 0, 0, 0);
const Foam::dimensionSet Foam::dimAcceleration(0, 1, -2, 0, 0, 0, 0);
const Foam::dimensionSet Foam::dimPressure(1, -1, -2, 0, 0, 0, 0);
const Foam::dimensionSet Foam::dimGasConstant(0, 2, -2, -1, 0, 0, 0);
const Foam::dimensionSet Foam::dimSpecificHeatCapacity(0, 2, -2, -1, 0, 0, 0);
const Foam::dimensionSet Foam::dimDensity(dimMass/dimVolume);
const Foam::dimensionSet Foam::dimForce(dimMass*dimAcceleration);
const Foam::dimensionSet Foam::dimEnergy(dimForce*dimLength);
const Foam::dimensionSet Foam::dimPower(dimEnergy/dimTime);
const Foam::dimensionSet Foam::dimPressure(dimForce/dimArea);
const Foam::dimensionSet Foam::dimGasConstant(dimEnergy/dimMass/dimTemperature);
const Foam::dimensionSet Foam::dimSpecificHeatCapacity(dimGasConstant);
// ************************************************************************* //

1
src/OpenFOAM/dimensionSet/dimensionSets.H
View File

@ -61,6 +61,7 @@ extern const dimensionSet dimVol;
extern const dimensionSet dimDensity;
extern const dimensionSet dimForce;
extern const dimensionSet dimEnergy;
extern const dimensionSet dimPower;
extern const dimensionSet dimVelocity;
extern const dimensionSet dimAcceleration;

19
src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
View File

@ -128,7 +128,20 @@ Foam::plane::plane(const vector& normalVector)
:
unitVector_(normalVector),
basePoint_(vector::zero)
{}
{
scalar magUnitVector(mag(unitVector_));
if (magUnitVector > VSMALL)
{
unitVector_ /= magUnitVector;
}
else
{
FatalErrorIn("plane::plane(const point&, const vector&)")
<< "plane normal has zero length"
<< abort(FatalError);
}
}
// Construct from point and normal vector
@ -146,7 +159,7 @@ Foam::plane::plane(const point& basePoint, const vector& normalVector)
else
{
FatalErrorIn("plane::plane(const point&, const vector&)")
<< "plane normal has got zero length"
<< "plane normal has zero length"
<< abort(FatalError);
}
}
@ -238,7 +251,7 @@ Foam::plane::plane(Istream& is)
else
{
FatalErrorIn("plane::plane(Istream& is)")
<< "plane normal has got zero length"
<< "plane normal has zero length"
<< abort(FatalError);
}
}

6
src/dynamicMesh/setUpdater/setUpdaterTemplates.C
View File

@ -77,11 +77,7 @@ void setUpdater::updateSets(const mapPolyMesh& morphMap) const
IOobjectList Objects
(
morphMap.mesh().time(),
morphMap.mesh().time().findInstance
(
morphMap.mesh().meshDir(),
"faces"
),
morphMap.mesh().facesInstance(),
"polyMesh/sets"
);

262
src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchField.C
View File

@ -1,262 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "expDirectionMixedFvPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
expDirectionMixedFvPatchField<Type>::expDirectionMixedFvPatchField
(
const fvPatch& p,
const DimensionedField<Type, volMesh>& iF
)
:
fvPatchField<Type>(p, iF),
refValue_(p.size()),
refGrad_(p.size()),
valueFraction_(p.size())
{}
template<class Type>
expDirectionMixedFvPatchField<Type>::expDirectionMixedFvPatchField
(
const expDirectionMixedFvPatchField<Type>& ptf,
const fvPatch& p,
const DimensionedField<Type, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fvPatchField<Type>(ptf, p, iF, mapper),
refValue_(ptf.refValue_, mapper),
refGrad_(ptf.refGrad_, mapper),
valueFraction_(ptf.valueFraction_, mapper)
{}
template<class Type>
expDirectionMixedFvPatchField<Type>::expDirectionMixedFvPatchField
(
const fvPatch& p,
const DimensionedField<Type, volMesh>& iF,
const dictionary& dict
)
:
fvPatchField<Type>(p, iF, dict),
refValue_("refValue", dict, p.size()),
refGrad_("refGradient", dict, p.size()),
valueFraction_("valueFraction", dict, p.size())
{
evaluate();
}
template<class Type>
expDirectionMixedFvPatchField<Type>::expDirectionMixedFvPatchField
(
const expDirectionMixedFvPatchField<Type>& ptf
)
:
fvPatchField<Type>(ptf),
refValue_(ptf.refValue_),
refGrad_(ptf.refGrad_),
valueFraction_(ptf.valueFraction_)
{}
template<class Type>
expDirectionMixedFvPatchField<Type>::expDirectionMixedFvPatchField
(
const expDirectionMixedFvPatchField<Type>& ptf,
const DimensionedField<Type, volMesh>& iF
)
:
fvPatchField<Type>(ptf, iF),
refValue_(ptf.refValue_),
refGrad_(ptf.refGrad_),
valueFraction_(ptf.valueFraction_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void expDirectionMixedFvPatchField<Type>::autoMap
(
const fvPatchFieldMapper& m
)
{
fvPatchField<Type>::autoMap(m);
refValue_.autoMap(m);
refGrad_.autoMap(m);
valueFraction_.autoMap(m);
}
// Reverse-map the given fvPatchField onto this fvPatchField
template<class Type>
void expDirectionMixedFvPatchField<Type>::rmap
(
const fvPatchField<Type>& ptf,
const labelList& addr
)
{
fvPatchField<Type>::rmap(ptf, addr);
const expDirectionMixedFvPatchField<Type>& edmptf =
refCast<const expDirectionMixedFvPatchField<Type> >(ptf);
refValue_.rmap(edmptf.refValue_, addr);
refGrad_.rmap(edmptf.refGrad_, addr);
valueFraction_.rmap(edmptf.valueFraction_, addr);
}
template<class Type>
tmp<Field<Type> > expDirectionMixedFvPatchField<Type>::snGrad() const
{
const vectorField& nHat = patch().faceNormals();
Field<Type> gradValue =
patchInternalField() + refGrad_/patch().deltaCoeffs();
Field<Type> mixedValue =
nHat*(nHat & refValue_)
+ gradValue - nHat*(nHat & gradValue);
return
valueFraction_*
(mixedValue - patchInternalField())*patch().deltaCoeffs()
+ (1.0 - valueFraction_)*refGrad_;
}
template<class Type>
void expDirectionMixedFvPatchField<Type>::evaluate(const Pstream::commsTypes)
{
if (!updated())
{
updateCoeffs();
}
const vectorField& nHat = patch().faceNormals();
Field<Type> gradValue =
patchInternalField() + refGrad_/patch().deltaCoeffs();
Field<Type> mixedValue =
nHat*(nHat & refValue_)
+ gradValue - nHat*(nHat & gradValue);
Field<Type>::operator=
(
valueFraction_*mixedValue + (1.0 - valueFraction_)*gradValue
);
fvPatchField<Type>::evaluate();
}
template<class Type>
tmp<Field<Type> > expDirectionMixedFvPatchField<Type>::valueInternalCoeffs
(
const tmp<scalarField>&
) const
{
return Type(pTraits<Type>::one)*(1.0 - valueFraction_);
}
template<class Type>
tmp<Field<Type> > expDirectionMixedFvPatchField<Type>::valueBoundaryCoeffs
(
const tmp<scalarField>&
) const
{
const vectorField& nHat = patch().faceNormals();
Field<Type> gradValue =
patchInternalField() + refGrad_/patch().deltaCoeffs();
Field<Type> mixedValue =
nHat*(nHat & refValue_)
+ gradValue - nHat*(nHat & gradValue);
return
valueFraction_*mixedValue
+ (1.0 - valueFraction_)*refGrad_/patch().deltaCoeffs();
}
template<class Type>
tmp<Field<Type> > expDirectionMixedFvPatchField<Type>::
gradientInternalCoeffs() const
{
return -Type(pTraits<Type>::one)*valueFraction_*patch().deltaCoeffs();
}
template<class Type>
tmp<Field<Type> > expDirectionMixedFvPatchField<Type>::
gradientBoundaryCoeffs() const
{
const vectorField& nHat = patch().faceNormals();
Field<Type> gradValue =
patchInternalField() + refGrad_/patch().deltaCoeffs();
Field<Type> mixedValue =
nHat*(nHat & refValue_)
+ gradValue - nHat*(nHat & gradValue);
return
valueFraction_*patch().deltaCoeffs()*mixedValue
+ (1.0 - valueFraction_)*refGrad_;
}
template<class Type>
void expDirectionMixedFvPatchField<Type>::write(Ostream& os) const
{
fvPatchField<Type>::write(os);
refValue_.writeEntry("refValue", os);
refGrad_.writeEntry("refGradient", os);
valueFraction_.writeEntry("valueFraction", os);
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

271
src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchField.H
View File

@ -1,271 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::expDirectionMixedFvPatchField
Description
Foam::expDirectionMixedFvPatchField
SourceFiles
expDirectionMixedFvPatchField.C
\*---------------------------------------------------------------------------*/
#ifndef expDirectionMixedFvPatchField_H
#define expDirectionMixedFvPatchField_H
#include "fvPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class expDirectionMixedFvPatch Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class expDirectionMixedFvPatchField
:
public fvPatchField<Type>
{
// Private data
//- Value field
Field<Type> refValue_;
//- Normal gradient field
Field<Type> refGrad_;
//- Fraction (0-1) of value used for boundary condition
scalarField valueFraction_;
public:
//- Runtime type information
TypeName("expDirectionMixedMixed");
// Constructors
//- Construct from patch and internal field
expDirectionMixedFvPatchField
(
const fvPatch&,
const DimensionedField<Type, volMesh>&
);
//- Construct from patch, internal field and dictionary
expDirectionMixedFvPatchField
(
const fvPatch&,
const DimensionedField<Type, volMesh>&,
const dictionary&
);
//- Construct by mapping given expDirectionMixedFvPatchField
// onto a new patch
expDirectionMixedFvPatchField
(
const expDirectionMixedFvPatchField<Type>&,
const fvPatch&,
const DimensionedField<Type, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
expDirectionMixedFvPatchField
(
const expDirectionMixedFvPatchField<Type>&
);
//- Construct and return a clone
virtual tmp<fvPatchField<Type> > clone() const
{
return tmp<fvPatchField<Type> >
(
new expDirectionMixedFvPatchField<Type>(*this)
);
}
//- Construct as copy setting internal field reference
expDirectionMixedFvPatchField
(
const expDirectionMixedFvPatchField<Type>&,
const DimensionedField<Type, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchField<Type> > clone
(
const DimensionedField<Type, volMesh>& iF
) const
{
return tmp<fvPatchField<Type> >
(
new expDirectionMixedFvPatchField<Type>(*this, iF)
);
}
// Member functions
// Access
//- Return true if this patch field fixes a value.
// Needed to check if a level has to be specified while solving
// Poissons equations.
virtual bool fixesrefValue() const
{
return true;
}
// Return defining fields
virtual Field<Type>& refValue()
{
return refValue_;
}
virtual const Field<Type>& refValue() const
{
return refValue_;
}
virtual Field<Type>& refGrad()
{
return refGrad_;
}
virtual const Field<Type>& refGrad() const
{
return refGrad_;
}
virtual scalarField& valueFraction()
{
return valueFraction_;
}
virtual const scalarField& valueFraction() const
{
return valueFraction_;
}
// Mapping functions
//- Map (and resize as needed) from self given a mapping object
virtual void autoMap
(
const fvPatchFieldMapper&
);
//- Reverse map the given fvPatchField onto this fvPatchField
virtual void rmap
(
const fvPatchField<Type>&,
const labelList&
);
// Evaluation functions
//- Return gradient at boundary
virtual tmp<Field<Type> > snGrad() const;
//- Evaluate the patch field
virtual void evaluate
(
const Pstream::commsTypes commsType=Pstream::blocking
);
//- Return the matrix diagonal coefficients corresponding to the
// evaluation of the value of this patchField with given weights
virtual tmp<Field<Type> > valueInternalCoeffs
(
const tmp<scalarField>&
) const;
//- Return the matrix source coefficients corresponding to the
// evaluation of the value of this patchField with given weights
virtual tmp<Field<Type> > valueBoundaryCoeffs
(
const tmp<scalarField>&
) const;
//- Return the matrix diagonal coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type> > gradientInternalCoeffs() const;
//- Return the matrix source coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type> > gradientBoundaryCoeffs() const;
//- Write
virtual void write(Ostream&) const;
// Member operators
virtual void operator=(const fvPatchField<Type>&) {}
virtual void operator+=(const fvPatchField<Type>&) {}
virtual void operator-=(const fvPatchField<Type>&) {}
virtual void operator*=(const fvPatchField<scalar>&) {}
virtual void operator/=(const fvPatchField<scalar>&) {}
virtual void operator=(const Field<Type>&) {}
virtual void operator+=(const Field<Type>&) {}
virtual void operator-=(const Field<Type>&) {}
virtual void operator*=(const Field<scalar>&) {}
virtual void operator/=(const Field<scalar>&) {}
virtual void operator=(const Type&) {}
virtual void operator+=(const Type&) {}
virtual void operator-=(const Type&) {}
virtual void operator*=(const scalar) {}
virtual void operator/=(const scalar) {}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "expDirectionMixedFvPatchField.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/AverageIOField.C
View File

@ -51,7 +51,7 @@ Foam::AverageIOField<Type>::AverageIOField
:
regIOobject(io),
Field<Type>(size),
average_(0)
average_(pTraits<Type>::zero)
{}

1
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
View File

@ -539,6 +539,7 @@ void Foam::autoHexMeshDriver::doMesh()
shrinkDict,
motionDict,
layerParams,
true, // pre-balance
decomposerPtr_(),
distributorPtr_()
);

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

@ -124,6 +124,7 @@ Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
}
Pout<< "Writing all faces to be merged to set "
<< allSets.objectPath() << endl;
allSets.instance() = mesh.time().timeName();
allSets.write();
}
@ -221,6 +222,7 @@ Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
if (debug)
{
errorFaces.instance() = mesh.time().timeName();
Pout<< "Writing all faces in error to faceSet "
<< errorFaces.objectPath() << nl << endl;
errorFaces.write();
@ -266,6 +268,7 @@ Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
if (debug)
{
faceSet restoreSet(mesh, "mastersToRestore", mastersToRestore);
restoreSet.instance() = mesh.time().timeName();
Pout<< "Writing all " << mastersToRestore.size()
<< " masterfaces to be restored to set "
<< restoreSet.objectPath() << endl;
@ -624,6 +627,7 @@ Foam::label Foam::autoLayerDriver::mergeEdgesUndo
if (debug)
{
errorFaces.instance() = mesh.time().timeName();
Pout<< "**Writing all faces in error to faceSet "
<< errorFaces.objectPath() << nl << endl;
errorFaces.write();
@ -806,6 +810,7 @@ void Foam::autoLayerDriver::checkMeshManifold() const
<< " points where this happens to pointSet "
<< nonManifoldPoints.name() << endl;
nonManifoldPoints.instance() = mesh.time().timeName();
nonManifoldPoints.write();
}
Info<< endl;
@ -945,6 +950,7 @@ void Foam::autoLayerDriver::handleNonManifolds
<< "and setting layer thickness to zero on these points."
<< endl;
nonManifoldPoints.instance() = mesh.time().timeName();
nonManifoldPoints.write();
forAll(meshPoints, patchPointI)
@ -1207,6 +1213,7 @@ void Foam::autoLayerDriver::handleWarpedFaces
//
// if (nMultiPatchCells > 0)
// {
// multiPatchCells.instance() = mesh.time().timeName();
// Info<< "Writing " << nMultiPatchCells
// << " cells with multiple (connected) patch faces to cellSet "
// << multiPatchCells.objectPath() << endl;
@ -2490,18 +2497,13 @@ void Foam::autoLayerDriver::mergePatchFacesUndo
<< "---------------------------" << nl
<< " - which are on the same patch" << nl
<< " - which make an angle < " << layerParams.featureAngle()
<< "- which are on the same patch" << nl
<< "- which make an angle < " << layerParams.featureAngle()
<< " degrees"
<< nl
<< " (cos:" << minCos << ')' << nl
<< " - as long as the resulting face doesn't become concave"
<< " (cos:" << minCos << ')' << nl
<< "- as long as the resulting face doesn't become concave"
<< " by more than "
<< layerParams.concaveAngle() << " degrees" << nl
<< " (0=straight, 180=fully concave)" << nl
<< " (0=straight, 180=fully concave)" << nl
<< endl;
label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
@ -3234,6 +3236,7 @@ void Foam::autoLayerDriver::doLayers
const dictionary& shrinkDict,
const dictionary& motionDict,
const layerParameters& layerParams,
const bool preBalance,
decompositionMethod& decomposer,
fvMeshDistribute& distributor
)
@ -3292,6 +3295,7 @@ void Foam::autoLayerDriver::doLayers
// Balance
if (Pstream::parRun())
if (Pstream::parRun() && preBalance)
{
Info<< nl
<< "Doing initial balancing" << nl

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

@ -535,6 +535,7 @@ public:
const dictionary& shrinkDict,
const dictionary& motionDict,
const layerParameters& layerParams,
const bool preBalance, // balance before adding?
decompositionMethod& decomposer,
fvMeshDistribute& distributor
);

69
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
View File

@ -205,6 +205,27 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
const_cast<Time&>(mesh.time())++;
}
if
(
returnReduce
(
(mesh.nCells() >= refineParams.maxLocalCells()),
orOp<bool>()
)
)
{
meshRefiner_.balanceAndRefine
(
"feature refinement iteration " + name(iter),
decomposer_,
distributor_,
cellsToRefine,
refineParams.maxLoadUnbalance()
);
}
else
{
meshRefiner_.refineAndBalance
(
"feature refinement iteration " + name(iter),
@ -215,6 +236,7 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
);
}
}
}
return iter;
}
@ -306,6 +328,27 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
const_cast<Time&>(mesh.time())++;
}
if
(
returnReduce
(
(mesh.nCells() >= refineParams.maxLocalCells()),
orOp<bool>()
)
)
{
meshRefiner_.balanceAndRefine
(
"surface refinement iteration " + name(iter),
decomposer_,
distributor_,
cellsToRefine,
refineParams.maxLoadUnbalance()
);
}
else
{
meshRefiner_.refineAndBalance
(
"surface refinement iteration " + name(iter),
@ -315,6 +358,7 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
refineParams.maxLoadUnbalance()
);
}
}
return iter;
}
@ -415,7 +459,9 @@ Foam::label Foam::autoRefineDriver::shellRefine
Pout<< "Dumping " << candidateCells.size()
<< " cells to cellSet candidateCellsFromShells." << endl;
cellSet(mesh, "candidateCellsFromShells", candidateCells).write();
cellSet c(mesh, "candidateCellsFromShells", candidateCells);
c.instance() = mesh.time().timeName();
c.write();
}
// Problem choosing starting faces for bufferlayers (bFaces)
@ -489,6 +535,26 @@ Foam::label Foam::autoRefineDriver::shellRefine
const_cast<Time&>(mesh.time())++;
}
if
(
returnReduce
(
(mesh.nCells() >= refineParams.maxLocalCells()),
orOp<bool>()
)
)
{
meshRefiner_.balanceAndRefine
(
"shell refinement iteration " + name(iter),
decomposer_,
distributor_,
cellsToRefine,
refineParams.maxLoadUnbalance()
);
}
else
{
meshRefiner_.refineAndBalance
(
"shell refinement iteration " + name(iter),
@ -498,6 +564,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
refineParams.maxLoadUnbalance()
);
}
}
meshRefiner_.userFaceData().clear();
return iter;

10
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
View File

@ -656,6 +656,16 @@ public:
const scalar maxLoadUnbalance
);
//- Balance before refining some cells
autoPtr<mapDistributePolyMesh> balanceAndRefine
(
const string& msg,
decompositionMethod& decomposer,
fvMeshDistribute& distributor,
const labelList& cellsToRefine,
const scalar maxLoadUnbalance
);
// Baffle handling

2
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
View File

@ -598,6 +598,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::getDuplicateFaces
}
Pout<< "Writing duplicate faces (baffles) to faceSet "
<< duplicateFaceSet.name() << nl << endl;
duplicateFaceSet.instance() = mesh_.time().timeName();
duplicateFaceSet.write();
}
@ -1559,6 +1560,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
problemTopo.insert(faceI);
}
}
problemTopo.instance() = mesh_.time().timeName();
Pout<< "Dumping " << problemTopo.size()
<< " problem faces to " << problemTopo.objectPath() << endl;
problemTopo.write();

3
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
View File

@ -203,6 +203,7 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
if (debug)
{
faceSet fSet(mesh_, "edgeConnectedFaces", candidateFaces);
fSet.instance() = mesh_.time().timeName();
Pout<< "Writing " << fSet.size()
<< " with problematic topology to faceSet "
<< fSet.objectPath() << endl;
@ -265,6 +266,7 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
if (debug)
{
perpFaces.instance() = mesh_.time().timeName();
Pout<< "Writing " << perpFaces.size()
<< " faces that are perpendicular to the surface to set "
<< perpFaces.objectPath() << endl;
@ -482,6 +484,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
if (debug)
{
cellSet problemCellSet(mesh_, "problemCells", problemCells.toc());
problemCellSet.instance() = mesh_.time().timeName();
Pout<< "Writing " << problemCellSet.size()
<< " cells that are edge connected to coarser cell to set "
<< problemCellSet.objectPath() << endl;

184
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
View File

@ -1245,90 +1245,110 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
}
//// Do refinement of consistent set of cells followed by truncation and
//// load balancing.
//Foam::autoPtr<Foam::mapDistributePolyMesh>
//Foam::meshRefinement::refineAndBalance
//(
// const string& msg,
// decompositionMethod& decomposer,
// fvMeshDistribute& distributor,
// const labelList& cellsToRefine
//)
//{
// // Do all refinement
// refine(cellsToRefine);
//
// if (debug)
// {
// Pout<< "Writing refined but unbalanced " << msg
// << " mesh to time " << timeName() << endl;
// write
// (
// debug,
// mesh_.time().path()
// /timeName()
// );
// Pout<< "Dumped debug data in = "
// << mesh_.time().cpuTimeIncrement() << " s" << endl;
//
// // test all is still synced across proc patches
// checkData();
// }
//
// Info<< "Refined mesh in = "
// << mesh_.time().cpuTimeIncrement() << " s" << endl;
// printMeshInfo(debug, "After refinement " + msg);
//
//
// // Load balancing
// // ~~~~~~~~~~~~~~
//
// autoPtr<mapDistributePolyMesh> distMap;
//
// if (Pstream::nProcs() > 1)
// {
// scalarField cellWeights(mesh_.nCells(), 1);
//
// distMap = balance
// (
// false, //keepZoneFaces
// false, //keepBaffles
// cellWeights,
// decomposer,
// distributor
// );
//
// Info<< "Balanced mesh in = "
// << mesh_.time().cpuTimeIncrement() << " s" << endl;
//
// printMeshInfo(debug, "After balancing " + msg);
//
//
// if (debug)
// {
// Pout<< "Writing balanced " << msg
// << " mesh to time " << timeName() << endl;
// write
// (
// debug,
// mesh_.time().path()/timeName()
// );
// Pout<< "Dumped debug data in = "
// << mesh_.time().cpuTimeIncrement() << " s" << endl;
//
// // test all is still synced across proc patches
// checkData();
// }
// }
//
// return distMap;
//}
// Do refinement of consistent set of cells followed by truncation and
// load balancing.
Foam::autoPtr<Foam::mapDistributePolyMesh>
Foam::meshRefinement::refineAndBalance
(
const string& msg,
decompositionMethod& decomposer,
fvMeshDistribute& distributor,
const labelList& cellsToRefine,
const scalar maxLoadUnbalance
)
{
// Do all refinement
refine(cellsToRefine);
if (debug)
{
Pout<< "Writing refined but unbalanced " << msg
<< " mesh to time " << timeName() << endl;
write
(
debug,
mesh_.time().path()
/timeName()
);
Pout<< "Dumped debug data in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
// test all is still synced across proc patches
checkData();
}
Info<< "Refined mesh in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
printMeshInfo(debug, "After refinement " + msg);
// Load balancing
// ~~~~~~~~~~~~~~
autoPtr<mapDistributePolyMesh> distMap;
if (Pstream::nProcs() > 1)
{
scalar nIdealCells =
mesh_.globalData().nTotalCells()
/ Pstream::nProcs();
scalar unbalance = returnReduce
(
mag(1.0-mesh_.nCells()/nIdealCells),
maxOp<scalar>()
);
if (unbalance <= maxLoadUnbalance)
{
Info<< "Skipping balancing since max unbalance " << unbalance
<< " is less than allowable " << maxLoadUnbalance
<< endl;
}
else
{
scalarField cellWeights(mesh_.nCells(), 1);
distMap = balance
(
false, //keepZoneFaces
false, //keepBaffles
cellWeights,
decomposer,
distributor
);
Info<< "Balanced mesh in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
printMeshInfo(debug, "After balancing " + msg);
if (debug)
{
Pout<< "Writing balanced " << msg
<< " mesh to time " << timeName() << endl;
write
(
debug,
mesh_.time().path()/timeName()
);
Pout<< "Dumped debug data in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
// test all is still synced across proc patches
checkData();
}
}
}
return distMap;
}
// Do load balancing followed by refinement of consistent set of cells.
Foam::autoPtr<Foam::mapDistributePolyMesh>
Foam::meshRefinement::refineAndBalance
Foam::meshRefinement::balanceAndRefine
(
const string& msg,
decompositionMethod& decomposer,
@ -1385,8 +1405,8 @@ Foam::meshRefinement::refineAndBalance
if (unbalance <= maxLoadUnbalance)
{
Info<< "Skipping balancing since max unbalance " << unbalance
<< " in = "
<< mesh_.time().cpuTimeIncrement() << " s" << endl;
<< " is less than allowable " << maxLoadUnbalance
<< endl;
}
else
{

403
src/meshTools/momentOfInertia/momentOfInertia.C Normal file
View File

@ -0,0 +1,403 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
momentOfInertia
Description
Reimplementation of volInt.c by Brian Mirtich.
* mirtich@cs.berkeley.edu *
* http://www.cs.berkeley.edu/~mirtich *
-------------------------------------------------------------------------------
*/
#include "momentOfInertia.H"
//#include "pyramidPointFaceRef.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
//Foam::tensor Foam::momentOfInertia
//(
// const pointField& points,
// const faceList& faces,
// const cell& cFaces,
// const point& cc
//)
//{
// tensor t(tensor::zero);
//
// forAll(cFaces, i)
// {
// const face& f = faces[cFaces[i]];
//
// scalar pyrVol = pyramidPointFaceRef(f, cc).mag(points);
//
// vector pyrCentre = pyramidPointFaceRef(f, cc).centre(points);
//
// vector d = pyrCentre - cc;
//
// t.xx() += pyrVol*(sqr(d.y()) + sqr(d.z()));
// t.yy() += pyrVol*(sqr(d.x()) + sqr(d.z()));
// t.zz() += pyrVol*(sqr(d.x()) + sqr(d.y()));
//
// t.xy() -= pyrVol*d.x()*d.y();
// t.xz() -= pyrVol*d.x()*d.z();
// t.yz() -= pyrVol*d.y()*d.z();
// }
//
// // Symmetric
// t.yx() = t.xy();
// t.zx() = t.xz();
// t.zy() = t.yz();
//
// return t;
//}
#define sqr(x) ((x)*(x))
#define pow3(x) ((x)*(x)*(x))
// compute various integrations over projection of face
void Foam::compProjectionIntegrals
(
const pointField& points,
const face& f,
const direction A,
const direction B,
scalar& P1,
scalar& Pa,
scalar& Pb,
scalar& Paa,
scalar& Pab,
scalar& Pbb,
scalar& Paaa,
scalar& Paab,
scalar& Pabb,
scalar& Pbbb
)
{
P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0;
forAll(f, i)
{
scalar a0 = points[f[i]][A];
scalar b0 = points[f[i]][B];
scalar a1 = points[f[(i+1) % f.size()]][A];
scalar b1 = points[f[(i+1) % f.size()]][B];
scalar da = a1 - a0;
scalar db = b1 - b0;
scalar a0_2 = a0 * a0;
scalar a0_3 = a0_2 * a0;
scalar a0_4 = a0_3 * a0;
scalar b0_2 = b0 * b0;
scalar b0_3 = b0_2 * b0;
scalar b0_4 = b0_3 * b0;
scalar a1_2 = a1 * a1;
scalar a1_3 = a1_2 * a1;
scalar b1_2 = b1 * b1;
scalar b1_3 = b1_2 * b1;
scalar C1 = a1 + a0;
scalar Ca = a1*C1 + a0_2;
scalar Caa = a1*Ca + a0_3;
scalar Caaa = a1*Caa + a0_4;
scalar Cb = b1*(b1 + b0) + b0_2;
scalar Cbb = b1*Cb + b0_3;
scalar Cbbb = b1*Cbb + b0_4;
scalar Cab = 3*a1_2 + 2*a1*a0 + a0_2;
scalar Kab = a1_2 + 2*a1*a0 + 3*a0_2;
scalar Caab = a0*Cab + 4*a1_3;
scalar Kaab = a1*Kab + 4*a0_3;
scalar Cabb = 4*b1_3 + 3*b1_2*b0 + 2*b1*b0_2 + b0_3;
scalar Kabb = b1_3 + 2*b1_2*b0 + 3*b1*b0_2 + 4*b0_3;
P1 += db*C1;
Pa += db*Ca;
Paa += db*Caa;
Paaa += db*Caaa;
Pb += da*Cb;
Pbb += da*Cbb;
Pbbb += da*Cbbb;
Pab += db*(b1*Cab + b0*Kab);
Paab += db*(b1*Caab + b0*Kaab);
Pabb += da*(a1*Cabb + a0*Kabb);
}
P1 /= 2.0;
Pa /= 6.0;
Paa /= 12.0;
Paaa /= 20.0;
Pb /= -6.0;
Pbb /= -12.0;
Pbbb /= -20.0;
Pab /= 24.0;
Paab /= 60.0;
Pabb /= -60.0;
}
void Foam::compFaceIntegrals
(
const pointField& points,
const face& f,
const vector& n,
const scalar w,
const direction A,
const direction B,
const direction C,
scalar& Fa,
scalar& Fb,
scalar& Fc,
scalar& Faa,
scalar& Fbb,
scalar& Fcc,
scalar& Faaa,
scalar& Fbbb,
scalar& Fccc,
scalar& Faab,
scalar& Fbbc,
scalar& Fcca
)
{
scalar P1, Pa, Pb, Paa, Pab, Pbb, Paaa, Paab, Pabb, Pbbb;
compProjectionIntegrals
(
points,
f,
A,
B,
P1,
Pa,
Pb,
Paa,
Pab,
Pbb,
Paaa,
Paab,
Pabb,
Pbbb
);
scalar k1 = 1 / n[C];
scalar k2 = k1 * k1;
scalar k3 = k2 * k1;
scalar k4 = k3 * k1;
Fa = k1 * Pa;
Fb = k1 * Pb;
Fc = -k2 * (n[A]*Pa + n[B]*Pb + w*P1);
Faa = k1 * Paa;
Fbb = k1 * Pbb;
Fcc = k3 * (sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb
+ w*(2*(n[A]*Pa + n[B]*Pb) + w*P1));
Faaa = k1 * Paaa;
Fbbb = k1 * Pbbb;
Fccc = -k4 * (pow3(n[A])*Paaa + 3*sqr(n[A])*n[B]*Paab
+ 3*n[A]*sqr(n[B])*Pabb + pow3(n[B])*Pbbb
+ 3*w*(sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb)
+ w*w*(3*(n[A]*Pa + n[B]*Pb) + w*P1));
Faab = k1 * Paab;
Fbbc = -k2 * (n[A]*Pabb + n[B]*Pbbb + w*Pbb);
Fcca = k3 * (sqr(n[A])*Paaa + 2*n[A]*n[B]*Paab + sqr(n[B])*Pabb
+ w*(2*(n[A]*Paa + n[B]*Pab) + w*Pa));
}
void Foam::compVolumeIntegrals
(
const pointField& points,
const faceList& faces,
const cell& cFaces,
const vectorField& fNorm,
const scalarField& fW,
scalar& T0,
vector& T1,
vector& T2,
vector& TP
)
{
T0 = 0;
T1 = vector::zero;
T2 = vector::zero;
TP = vector::zero;
forAll(cFaces, i)
{
const vector& n = fNorm[i];
scalar nx = mag(n[0]);
scalar ny = mag(n[1]);
scalar nz = mag(n[2]);
direction A, B, C;
if (nx > ny && nx > nz)
{
C = 0;
}
else
{
C = (ny > nz) ? 1 : 2;
}
A = (C + 1) % 3;
B = (A + 1) % 3;
scalar Fa, Fb, Fc, Faa, Fbb, Fcc, Faaa, Fbbb, Fccc, Faab, Fbbc, Fcca;
compFaceIntegrals
(
points,
faces[cFaces[i]],
n,
fW[i],
A,
B,
C,
Fa,
Fb,
Fc,
Faa,
Fbb,
Fcc,
Faaa,
Fbbb,
Fccc,
Faab,
Fbbc,
Fcca
);
T0 += n[0] * ((A == 0) ? Fa : ((B == 0) ? Fb : Fc));
T1[A] += n[A] * Faa;
T1[B] += n[B] * Fbb;
T1[C] += n[C] * Fcc;
T2[A] += n[A] * Faaa;
T2[B] += n[B] * Fbbb;
T2[C] += n[C] * Fccc;
TP[A] += n[A] * Faab;
TP[B] += n[B] * Fbbc;
TP[C] += n[C] * Fcca;
}
T1 /= 2;
T2 /= 3;
TP /= 2;
}
// Calculate
// - r: centre of mass
// - J: inertia around origin (point 0,0,0)
void Foam::momentOfIntertia
(
const pointField& points,
const faceList& faces,
const cell& cFaces,
point& r,
tensor& J
)
{
// Face normals
vectorField fNorm(cFaces.size());
scalarField fW(cFaces.size());
forAll(cFaces, i)
{
label faceI = cFaces[i];
const face& f = faces[faceI];
fNorm[i] = f.normal(points);
fNorm[i] /= mag(fNorm[i]) + VSMALL;
fW[i] = - (fNorm[i] & points[f[0]]);
}
scalar T0;
vector T1, T2, TP;
compVolumeIntegrals
(
points,
faces,
cFaces,
fNorm,
fW,
T0,
T1,
T2,
TP
);
const scalar density = 1.0; /* assume unit density */
scalar mass = density * T0;
/* compute center of mass */
r = T1 / T0;
/* compute inertia tensor */
J.xx() = density * (T2[1] + T2[2]);
J.yy() = density * (T2[2] + T2[0]);
J.zz() = density * (T2[0] + T2[1]);
J.xy() = J.yx() = - density * TP[0];
J.yz() = J.zy() = - density * TP[1];
J.zx() = J.xz() = - density * TP[2];
///* translate inertia tensor to center of mass */
//J[XX] -= mass * (r[1]*r[1] + r[2]*r[2]);
//J[YY] -= mass * (r[2]*r[2] + r[0]*r[0]);
//J[ZZ] -= mass * (r[0]*r[0] + r[1]*r[1]);
//J[XY] = J[YX] += mass * r[0] * r[1];
//J[YZ] = J[ZY] += mass * r[1] * r[2];
//J[ZX] = J[XZ] += mass * r[2] * r[0];
}
// ************************************************************************* //

40
src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFieldsFwd.H → src/meshTools/momentOfInertia/momentOfInertia.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2007-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,23 +22,47 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
momentOfInertia
Description
SourceFiles
momentOfInertia.H
\*---------------------------------------------------------------------------*/
#ifndef expDirectionMixedFvPatchFieldsFwd_H
#define expDirectionMixedFvPatchFieldsFwd_H
#ifndef momentOfInertia_H
#define momentOfInertia_H
#include "fieldTypes.H"
#include "tensor.H"
#include "primitiveMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
////- Moment of inertia around cell centre for single cell.
//tensor momentOfInertia
//(
// const pointField&,
// const faceList&,
// const cell&,
// const point& cc
//);
template<class Type> class expDirectionMixedFvPatchField;
makePatchTypeFieldTypedefs(expDirectionMixed)
// Calculate
// - centre of mass
// - inertia tensor around (0,0,0)
void momentOfIntertia
(
const pointField&,
const faceList&,
const cell&,
point& r,
tensor& Jorigin
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

BIN
src/meshTools/momentOfInertia/volInt.ps.gz Normal file
View File

Binary file not shown.

146
src/meshTools/momentOfInertia/volumeIntegration/README Normal file
View File

@ -0,0 +1,146 @@
1. OVERVIEW
This code accompanies the paper:
Brian Mirtich, "Fast and Accurate Computation of
Polyhedral Mass Properties," journal of graphics
tools, volume 1, number 2, 1996.
It computes the ten volume integrals needed for
determining the center of mass, moments of
inertia, and products of inertia for a uniform
density polyhedron. From this information, a
body frame can be computed.
To compile the program, use an ANSI compiler, and
type something like
% cc volInt.c -O2 -lm -o volInt
Revision history
26 Jan 1996 Program creation.
3 Aug 1996 Corrected bug arising when polyhedron density
is not 1.0. Changes confined to function main().
Thanks to Zoran Popovic for catching this one.
2. POLYHEDRON GEOMETRY FILES
The program reads a data file specified on the
command line. This data file describes the
geometry of a polyhedron, and has the following
format:
N
x_0 y_0 z_0
x_1 y_1 z_1
.
.
.
x_{N-1} y_{N-1} z_{N-1}
M
k1 v_{1,1} v_{1,2} ... v_{1,k1}
k2 v_{2,1} v_{2,2} ... v_{2,k2}
.
.
.
kM v_{M,1} v_{M,2} ... v_{M,kM}
where:
N number of vertices on polyhedron
x_i y_i z_i x, y, and z coordinates of ith vertex
M number of faces on polyhedron
ki number of vertices on ith face
v_{i,j} jth vertex on ith face
In English:
First the number of vertices are specified. Next
the vertices are defined by listing the 3
coordinates of each one. Next the number of faces
are specified. Finally, the faces themselves are
specified. A face is specified by first giving
the number of vertices around the polygonal face,
followed by the (integer) indices of these
vertices. When specifying indices, note that
they must be given in counter-clockwise order
(when looking at the face from outside the
polyhedron), and the vertices are indexed from 0
to N-1 for a polyhedron with N faces.
White space can be inserted (or not) as desired.
Three example polyhedron geometry files are included:
cube A cube, 20 units on a side, centered at
the origin and aligned with the coordinate axes.
tetra A tetrahedron formed by taking the convex
hull of the origin, and the points (5,0,0),
(0,4,0), and (0,0,3).
icosa An icosahedron with vertices lying on the unit
sphere, centered at the origin.
3. RUNNING THE PROGRAM
Simply type,
% volInt <polyhedron geometry filename>
The program will read in the geometry of the
polyhedron, and the print out the ten volume
integrals.
The program also computes some of the mass
properties which may be inferred from the volume
integrals. A density of 1.0 is assumed, although
this may be changed in function main(). The
center of mass is computed as well as the inertia
tensor relative to a frame with origin at the
center of mass. Note, however, that this matrix
may not be diagonal. If not, a diagonalization
routine must be performed, to determine the
orientation of the true body frame relative to
the original model frame. The Jacobi method
works quite well (see Numerical Recipes in C by
Press, et. al.).
4. DISCLAIMERS
1. The volume integration code has been written
to match the development and algorithms presented
in the paper, and not with maximum optimization
in mind. While inherently very efficient, a few
more cycles can be squeezed out of the algorithm.
This is left as an exercise. :)
2. Don't like global variables? The three
procedures which evaluate the volume integrals
can be combined into a single procedure with two
nested loops. In addition to providing some
speedup, all of the global variables can then be
made local.
3. The polyhedron data structure used by the
program is admittedly lame; much better schemes
are possible. The idea here is just to give the
basic integral evaluation code, which will have
to be adjusted for other polyhedron data
structures.
4. There is no error checking for the input
files. Be careful. Note the hard limits
#defined for the number of vertices, number of
faces, and number of vertices per faces.

21
src/meshTools/momentOfInertia/volumeIntegration/cube Normal file
View File

@ -0,0 +1,21 @@
8
-10 -10 -10
+10 -10 -10
+10 +10 -10
-10 +10 -10
-10 -10 +10
+10 -10 +10
+10 +10 +10
-10 +10 +10
6
4 0 3 2 1
4 4 5 6 7
4 0 1 5 4
4 6 2 3 7
4 1 2 6 5
4 0 4 7 3

38
src/meshTools/momentOfInertia/volumeIntegration/icosa Normal file
View File

@ -0,0 +1,38 @@
12
+0.00000000 +0.00000000 +1.00000000
+0.00000000 +0.00000000 -1.00000000
+0.89442719 +0.00000000 +0.44721360
+0.27639320 +0.85065081 +0.44721360
-0.72360680 +0.52573111 +0.44721360
-0.72360680 -0.52573111 +0.44721360
+0.27639320 -0.85065081 +0.44721360
+0.72360680 +0.52573111 -0.44721360
-0.27639320 +0.85065081 -0.44721360
-0.89442719 +0.00000000 -0.44721360
-0.27639320 -0.85065081 -0.44721360
+0.72360680 -0.52573111 -0.44721360
20
3 6 11 2
3 3 2 7
3 7 2 11
3 0 2 3
3 0 6 2
3 10 1 11
3 1 7 11
3 10 11 6
3 8 7 1
3 8 3 7
3 5 10 6
3 5 6 0
3 4 3 8
3 4 0 3
3 9 8 1
3 9 1 10
3 4 5 0
3 9 10 5
3 9 5 4
3 9 4 8

16
src/meshTools/momentOfInertia/volumeIntegration/tetra Normal file
View File

@ -0,0 +1,16 @@
4
0 0 0
5 0 0
0 4 0
0 0 3
4
3 0 3 2
3 3 0 1
3 2 1 0
3 1 2 3

380
src/meshTools/momentOfInertia/volumeIntegration/volInt.c Normal file
View File

@ -0,0 +1,380 @@
/*******************************************************
* *
* volInt.c *
* *
* This code computes volume integrals needed for *
* determining mass properties of polyhedral bodies. *
* *
* For more information, see the accompanying README *
* file, and the paper *
* *
* Brian Mirtich, "Fast and Accurate Computation of *
* Polyhedral Mass Properties," journal of graphics *
* tools, volume 1, number 1, 1996. *
* *
* This source code is public domain, and may be used *
* in any way, shape or form, free of charge. *
* *
* Copyright 1995 by Brian Mirtich *
* *
* mirtich@cs.berkeley.edu *
* http://www.cs.berkeley.edu/~mirtich *
* *
*******************************************************/
/*
Revision history
26 Jan 1996 Program creation.
3 Aug 1996 Corrected bug arising when polyhedron density
is not 1.0. Changes confined to function main().
Thanks to Zoran Popovic for catching this one.
27 May 1997 Corrected sign error in translation of inertia
product terms to center of mass frame. Changes
confined to function main(). Thanks to
Chris Hecker.
*/
#include <stdio.h>
#include <string.h>
#include <math.h>
/*
============================================================================
constants
============================================================================
*/
#define MAX_VERTS 100 /* maximum number of polyhedral vertices */
#define MAX_FACES 100 /* maximum number of polyhedral faces */
#define MAX_POLYGON_SZ 10 /* maximum number of verts per polygonal face */
#define X 0
#define Y 1
#define Z 2
/*
============================================================================
macros
============================================================================
*/
#define SQR(x) ((x)*(x))
#define CUBE(x) ((x)*(x)*(x))
/*
============================================================================
data structures
============================================================================
*/
typedef struct {
int numVerts;
double norm[3];
double w;
int verts[MAX_POLYGON_SZ];
struct polyhedron *poly;
} FACE;
typedef struct polyhedron {
int numVerts, numFaces;
double verts[MAX_VERTS][3];
FACE faces[MAX_FACES];
} POLYHEDRON;
/*
============================================================================
globals
============================================================================
*/
static int A; /* alpha */
static int B; /* beta */
static int C; /* gamma */
/* projection integrals */
static double P1, Pa, Pb, Paa, Pab, Pbb, Paaa, Paab, Pabb, Pbbb;
/* face integrals */
static double Fa, Fb, Fc, Faa, Fbb, Fcc, Faaa, Fbbb, Fccc, Faab, Fbbc, Fcca;
/* volume integrals */
static double T0, T1[3], T2[3], TP[3];
/*
============================================================================
read in a polyhedron
============================================================================
*/
void readPolyhedron(char *name, POLYHEDRON *p)
{
FILE *fp;
char line[200], *c;
int i, j, n;
double dx1, dy1, dz1, dx2, dy2, dz2, nx, ny, nz, len;
FACE *f;
if (!(fp = fopen(name, "r"))) {
printf("i/o error\n");
exit(1);
}
fscanf(fp, "%d", &p->numVerts);
printf("Reading in %d vertices\n", p->numVerts);
for (i = 0; i < p->numVerts; i++)
fscanf(fp, "%lf %lf %lf",
&p->verts[i][X], &p->verts[i][Y], &p->verts[i][Z]);
fscanf(fp, "%d", &p->numFaces);
printf("Reading in %d faces\n", p->numFaces);
for (i = 0; i < p->numFaces; i++) {
f = &p->faces[i];
f->poly = p;
fscanf(fp, "%d", &f->numVerts);
for (j = 0; j < f->numVerts; j++) fscanf(fp, "%d", &f->verts[j]);
/* compute face normal and offset w from first 3 vertices */
dx1 = p->verts[f->verts[1]][X] - p->verts[f->verts[0]][X];
dy1 = p->verts[f->verts[1]][Y] - p->verts[f->verts[0]][Y];
dz1 = p->verts[f->verts[1]][Z] - p->verts[f->verts[0]][Z];
dx2 = p->verts[f->verts[2]][X] - p->verts[f->verts[1]][X];
dy2 = p->verts[f->verts[2]][Y] - p->verts[f->verts[1]][Y];
dz2 = p->verts[f->verts[2]][Z] - p->verts[f->verts[1]][Z];
nx = dy1 * dz2 - dy2 * dz1;
ny = dz1 * dx2 - dz2 * dx1;
nz = dx1 * dy2 - dx2 * dy1;
len = sqrt(nx * nx + ny * ny + nz * nz);
f->norm[X] = nx / len;
f->norm[Y] = ny / len;
f->norm[Z] = nz / len;
f->w = - f->norm[X] * p->verts[f->verts[0]][X]
- f->norm[Y] * p->verts[f->verts[0]][Y]
- f->norm[Z] * p->verts[f->verts[0]][Z];
}
fclose(fp);
}
/*
============================================================================
compute mass properties
============================================================================
*/
/* compute various integrations over projection of face */
void compProjectionIntegrals(FACE *f)
{
double a0, a1, da;
double b0, b1, db;
double a0_2, a0_3, a0_4, b0_2, b0_3, b0_4;
double a1_2, a1_3, b1_2, b1_3;
double C1, Ca, Caa, Caaa, Cb, Cbb, Cbbb;
double Cab, Kab, Caab, Kaab, Cabb, Kabb;
int i;
P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0;
for (i = 0; i < f->numVerts; i++) {
a0 = f->poly->verts[f->verts[i]][A];
b0 = f->poly->verts[f->verts[i]][B];
a1 = f->poly->verts[f->verts[(i+1) % f->numVerts]][A];
b1 = f->poly->verts[f->verts[(i+1) % f->numVerts]][B];
da = a1 - a0;
db = b1 - b0;
a0_2 = a0 * a0; a0_3 = a0_2 * a0; a0_4 = a0_3 * a0;
b0_2 = b0 * b0; b0_3 = b0_2 * b0; b0_4 = b0_3 * b0;
a1_2 = a1 * a1; a1_3 = a1_2 * a1;
b1_2 = b1 * b1; b1_3 = b1_2 * b1;
C1 = a1 + a0;
Ca = a1*C1 + a0_2; Caa = a1*Ca + a0_3; Caaa = a1*Caa + a0_4;
Cb = b1*(b1 + b0) + b0_2; Cbb = b1*Cb + b0_3; Cbbb = b1*Cbb + b0_4;
Cab = 3*a1_2 + 2*a1*a0 + a0_2; Kab = a1_2 + 2*a1*a0 + 3*a0_2;
Caab = a0*Cab + 4*a1_3; Kaab = a1*Kab + 4*a0_3;
Cabb = 4*b1_3 + 3*b1_2*b0 + 2*b1*b0_2 + b0_3;
Kabb = b1_3 + 2*b1_2*b0 + 3*b1*b0_2 + 4*b0_3;
P1 += db*C1;
Pa += db*Ca;
Paa += db*Caa;
Paaa += db*Caaa;
Pb += da*Cb;
Pbb += da*Cbb;
Pbbb += da*Cbbb;
Pab += db*(b1*Cab + b0*Kab);
Paab += db*(b1*Caab + b0*Kaab);
Pabb += da*(a1*Cabb + a0*Kabb);
}
P1 /= 2.0;
Pa /= 6.0;
Paa /= 12.0;
Paaa /= 20.0;
Pb /= -6.0;
Pbb /= -12.0;
Pbbb /= -20.0;
Pab /= 24.0;
Paab /= 60.0;
Pabb /= -60.0;
}
compFaceIntegrals(FACE *f)
{
double *n, w;
double k1, k2, k3, k4;
compProjectionIntegrals(f);
w = f->w;
n = f->norm;
k1 = 1 / n[C]; k2 = k1 * k1; k3 = k2 * k1; k4 = k3 * k1;
Fa = k1 * Pa;
Fb = k1 * Pb;
Fc = -k2 * (n[A]*Pa + n[B]*Pb + w*P1);
Faa = k1 * Paa;
Fbb = k1 * Pbb;
Fcc = k3 * (SQR(n[A])*Paa + 2*n[A]*n[B]*Pab + SQR(n[B])*Pbb
+ w*(2*(n[A]*Pa + n[B]*Pb) + w*P1));
Faaa = k1 * Paaa;
Fbbb = k1 * Pbbb;
Fccc = -k4 * (CUBE(n[A])*Paaa + 3*SQR(n[A])*n[B]*Paab
+ 3*n[A]*SQR(n[B])*Pabb + CUBE(n[B])*Pbbb
+ 3*w*(SQR(n[A])*Paa + 2*n[A]*n[B]*Pab + SQR(n[B])*Pbb)
+ w*w*(3*(n[A]*Pa + n[B]*Pb) + w*P1));
Faab = k1 * Paab;
Fbbc = -k2 * (n[A]*Pabb + n[B]*Pbbb + w*Pbb);
Fcca = k3 * (SQR(n[A])*Paaa + 2*n[A]*n[B]*Paab + SQR(n[B])*Pabb
+ w*(2*(n[A]*Paa + n[B]*Pab) + w*Pa));
}
void compVolumeIntegrals(POLYHEDRON *p)
{
FACE *f;
double nx, ny, nz;
int i;
T0 = T1[X] = T1[Y] = T1[Z]
= T2[X] = T2[Y] = T2[Z]
= TP[X] = TP[Y] = TP[Z] = 0;
for (i = 0; i < p->numFaces; i++) {
f = &p->faces[i];
nx = fabs(f->norm[X]);
ny = fabs(f->norm[Y]);
nz = fabs(f->norm[Z]);
if (nx > ny && nx > nz) C = X;
else C = (ny > nz) ? Y : Z;
A = (C + 1) % 3;
B = (A + 1) % 3;
compFaceIntegrals(f);
T0 += f->norm[X] * ((A == X) ? Fa : ((B == X) ? Fb : Fc));
T1[A] += f->norm[A] * Faa;
T1[B] += f->norm[B] * Fbb;
T1[C] += f->norm[C] * Fcc;
T2[A] += f->norm[A] * Faaa;
T2[B] += f->norm[B] * Fbbb;
T2[C] += f->norm[C] * Fccc;
TP[A] += f->norm[A] * Faab;
TP[B] += f->norm[B] * Fbbc;
TP[C] += f->norm[C] * Fcca;
}
T1[X] /= 2; T1[Y] /= 2; T1[Z] /= 2;
T2[X] /= 3; T2[Y] /= 3; T2[Z] /= 3;
TP[X] /= 2; TP[Y] /= 2; TP[Z] /= 2;
}
/*
============================================================================
main
============================================================================
*/
int main(int argc, char *argv[])
{
POLYHEDRON p;
double density, mass;
double r[3]; /* center of mass */
double J[3][3]; /* inertia tensor */
if (argc != 2) {
printf("usage: %s <polyhedron geometry filename>\n", argv[0]);
exit(0);
}
readPolyhedron(argv[1], &p);
compVolumeIntegrals(&p);
printf("\nT1 = %+20.6f\n\n", T0);
printf("Tx = %+20.6f\n", T1[X]);
printf("Ty = %+20.6f\n", T1[Y]);
printf("Tz = %+20.6f\n\n", T1[Z]);
printf("Txx = %+20.6f\n", T2[X]);
printf("Tyy = %+20.6f\n", T2[Y]);
printf("Tzz = %+20.6f\n\n", T2[Z]);
printf("Txy = %+20.6f\n", TP[X]);
printf("Tyz = %+20.6f\n", TP[Y]);
printf("Tzx = %+20.6f\n\n", TP[Z]);
density = 1.0; /* assume unit density */
mass = density * T0;
/* compute center of mass */
r[X] = T1[X] / T0;
r[Y] = T1[Y] / T0;
r[Z] = T1[Z] / T0;
/* compute inertia tensor */
J[X][X] = density * (T2[Y] + T2[Z]);
J[Y][Y] = density * (T2[Z] + T2[X]);
J[Z][Z] = density * (T2[X] + T2[Y]);
J[X][Y] = J[Y][X] = - density * TP[X];
J[Y][Z] = J[Z][Y] = - density * TP[Y];
J[Z][X] = J[X][Z] = - density * TP[Z];
/* translate inertia tensor to center of mass */
J[X][X] -= mass * (r[Y]*r[Y] + r[Z]*r[Z]);
J[Y][Y] -= mass * (r[Z]*r[Z] + r[X]*r[X]);
J[Z][Z] -= mass * (r[X]*r[X] + r[Y]*r[Y]);
J[X][Y] = J[Y][X] += mass * r[X] * r[Y];
J[Y][Z] = J[Z][Y] += mass * r[Y] * r[Z];
J[Z][X] = J[X][Z] += mass * r[Z] * r[X];
printf("center of mass: (%+12.6f,%+12.6f,%+12.6f)\n\n", r[X], r[Y], r[Z]);
printf("inertia tensor with origin at c.o.m. :\n");
printf("%+15.6f %+15.6f %+15.6f\n", J[X][X], J[X][Y], J[X][Z]);
printf("%+15.6f %+15.6f %+15.6f\n", J[Y][X], J[Y][Y], J[Y][Z]);
printf("%+15.6f %+15.6f %+15.6f\n\n", J[Z][X], J[Z][Y], J[Z][Z]);
}

9
src/meshTools/searchableSurface/searchableSurfaceCollection.C
View File

@ -681,15 +681,17 @@ void Foam::searchableSurfaceCollection::getField
List<List<label> > infoMap;
sortHits(info, surfInfo, infoMap);
values.setSize(info.size());
//?Misses do not get set? values = 0;
// Do surface tests
forAll(surfInfo, surfI)
{
labelList surfValues;
subGeom_[surfI].getField(surfInfo[surfI], surfValues);
if (surfValues.size())
{
// Size values only when we have a surface that supports it.
values.setSize(info.size());
const labelList& map = infoMap[surfI];
forAll(map, i)
{
@ -697,6 +699,7 @@ void Foam::searchableSurfaceCollection::getField
}
}
}
}
}

43
src/meshTools/sets/topoSets/cellSet.C
View File

@ -117,7 +117,18 @@ cellSet::cellSet
IOobject
(
name,
runTime.findInstance(polyMesh::meshSubDir, "faces"),
runTime.findInstance
(
polyMesh::meshSubDir/"sets", //polyMesh::meshSubDir,
word::null, //"faces"
IOobject::MUST_READ,
runTime.findInstance
(
polyMesh::meshSubDir,
"faces",
IOobject::READ_IF_PRESENT
)
),
polyMesh::meshSubDir/"sets",
runTime,
r,
@ -141,10 +152,21 @@ cellSet::cellSet
IOobject
(
name,
runTime.findInstance(polyMesh::meshSubDir, "faces"),
runTime.findInstance
(
polyMesh::meshSubDir/"sets", //polyMesh::meshSubDir,
word::null, //"faces"
IOobject::NO_READ,
runTime.findInstance
(
polyMesh::meshSubDir,
"faces",
IOobject::READ_IF_PRESENT
)
),
polyMesh::meshSubDir/"sets",
runTime,
NO_READ,
IOobject::NO_READ,
w
),
size
@ -165,10 +187,21 @@ cellSet::cellSet
IOobject
(
name,
runTime.findInstance(polyMesh::meshSubDir, "faces"),
runTime.findInstance
(
polyMesh::meshSubDir/"sets", //polyMesh::meshSubDir,
word::null, //"faces"
IOobject::NO_READ,
runTime.findInstance
(
polyMesh::meshSubDir,
"faces",
IOobject::READ_IF_PRESENT
)
),
polyMesh::meshSubDir/"sets",
runTime,
NO_READ,
IOobject::NO_READ,
w
),
set

29
src/meshTools/sets/topoSets/topoSet.C
View File

@ -28,6 +28,7 @@ License
#include "mapPolyMesh.H"
#include "polyMesh.H"
#include "boundBox.H"
#include "Time.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -362,7 +363,13 @@ topoSet::topoSet
IOobject
(
name,
mesh.facesInstance(),
mesh.time().findInstance
(
polyMesh::meshSubDir/"sets",
word::null,
IOobject::MUST_READ,
mesh.facesInstance()
),
polyMesh::meshSubDir/"sets",
mesh,
r,
@ -402,10 +409,16 @@ topoSet::topoSet
IOobject
(
name,
mesh.facesInstance(),
mesh.time().findInstance
(
polyMesh::meshSubDir/"sets",
word::null,
IOobject::NO_READ,
mesh.facesInstance()
),
polyMesh::meshSubDir/"sets",
mesh,
NO_READ,
IOobject::NO_READ,
w
)
),
@ -426,10 +439,16 @@ topoSet::topoSet
IOobject
(
name,
mesh.facesInstance(),
mesh.time().findInstance
(
polyMesh::meshSubDir/"sets",
word::null,
IOobject::NO_READ,
mesh.facesInstance()
),
polyMesh::meshSubDir/"sets",
mesh,
NO_READ,
IOobject::NO_READ,
w
)
),

13
src/meshTools/sets/topoSets/topoSet.H
View File

@ -176,7 +176,8 @@ public:
// Can't use typeName info here since subclasses not yet instantiated
topoSet(const IOobject&, const word& wantedType);
//- Construct from polyMesh and name
//- Construct from polyMesh and name. Searches for a polyMesh/sets
// directory but not beyond mesh.facesInstance.
topoSet
(
const polyMesh& mesh,
@ -186,7 +187,9 @@ public:
writeOption w=NO_WRITE
);
//- Construct empty from additional size of labelHashSet
//- Construct empty from additional size of labelHashSet.
// Searches for a polyMesh/sets
// directory but not beyond mesh.facesInstance.
topoSet
(
const polyMesh& mesh,
@ -196,6 +199,8 @@ public:
);
//- Construct empty from additional labelHashSet
// Searches for a polyMesh/sets
// directory but not beyond mesh.facesInstance.
topoSet
(
const polyMesh& mesh,
@ -204,10 +209,10 @@ public:
writeOption w=NO_WRITE
);
//- Construct empty from IOobject and size
//- Construct empty from IOobject and size.
topoSet(const IOobject&, const label size);
//- Construct from IOobject and labelHashSet
//- Construct from IOobject and labelHashSet.
topoSet(const IOobject&, const labelHashSet&);

15
src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C
View File

@ -253,9 +253,19 @@ Foam::label Foam::scotchDecomp::decompose
// Check for externally provided cellweights and if so initialise weights
scalar maxWeights = gMax(cWeights);
scalar minWeights = gMin(cWeights);
if (cWeights.size() > 0)
{
if (minWeights <= 0)
{
WarningIn
(
"scotchDecomp::decompose"
"(const pointField&, const scalarField&)"
) << "Illegal minimum weight " << minWeights
<< endl;
}
if (cWeights.size() != xadj.size()-1)
{
FatalErrorIn
@ -266,11 +276,12 @@ Foam::label Foam::scotchDecomp::decompose
<< " does not equal number of cells " << xadj.size()-1
<< exit(FatalError);
}
// Convert to integers.
velotab.setSize(cWeights.size());
forAll(velotab, i)
{
velotab[i] = int(1000000*cWeights[i]/maxWeights);
velotab[i] = int(cWeights[i]/minWeights);
}
}

5
src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.H
View File

@ -104,6 +104,9 @@ public:
//- Return for every coordinate the wanted processor number. Use the
// mesh connectivity (if needed)
// Weights get truncated to convert into integer
// so e.g. 3.5 is seen as 3. The overall sum of weights
// might otherwise overflow.
virtual labelList decompose
(
const pointField& points,
@ -114,6 +117,7 @@ public:
// passed agglomeration map (from fine to coarse cells) and coarse cell
// location. Can be overridden by decomposers that provide this
// functionality natively.
// See note on weights above.
virtual labelList decompose
(
const labelList& agglom,
@ -128,6 +132,7 @@ public:
// from 0 at processor0 and then incrementing all through the
// processors)
// - the connections are across coupled patches
// See note on weights above.
virtual labelList decompose
(
const labelListList& globalCellCells,

14
src/parallel/metisDecomp/metisDecomp.C
View File

@ -89,9 +89,19 @@ Foam::label Foam::metisDecomp::decompose
// Check for externally provided cellweights and if so initialise weights
scalar maxWeights = gMax(cWeights);
scalar minWeights = gMin(cWeights);
if (cWeights.size() > 0)
{
if (minWeights <= 0)
{
WarningIn
(
"metisDecomp::decompose"
"(const pointField&, const scalarField&)"
) << "Illegal minimum weight " << minWeights
<< endl;
}
if (cWeights.size() != numCells)
{
FatalErrorIn
@ -106,7 +116,7 @@ Foam::label Foam::metisDecomp::decompose
cellWeights.setSize(cWeights.size());
forAll(cellWeights, i)
{
cellWeights[i] = int(1000000*cWeights[i]/maxWeights);
cellWeights[i] = int(cWeights[i]/minWeights);
}
}

5
src/parallel/metisDecomp/metisDecomp.H
View File

@ -101,6 +101,9 @@ public:
//- Return for every coordinate the wanted processor number. Use the
// mesh connectivity (if needed)
// Weights get normalised so the minimum value is 1 before truncation
// to an integer so the weights should be multiples of the minimum
// value. The overall sum of weights might otherwise overflow.
virtual labelList decompose
(
const pointField& points,
@ -111,6 +114,7 @@ public:
// passed agglomeration map (from fine to coarse cells) and coarse cell
// location. Can be overridden by decomposers that provide this
// functionality natively.
// See note on weights above.
virtual labelList decompose
(
const labelList& agglom,
@ -125,6 +129,7 @@ public:
// from 0 at processor0 and then incrementing all through the
// processors)
// - the connections are across coupled patches
// See note on weights above.
virtual labelList decompose
(
const labelListList& globalCellCells,

30
src/parallel/parMetisDecomp/parMetisDecomp.C
View File

@ -428,9 +428,19 @@ Foam::labelList Foam::parMetisDecomp::decompose
// Check for externally provided cellweights and if so initialise weights
scalar maxWeights = gMax(cWeights);
scalar minWeights = gMin(cWeights);
if (cWeights.size() > 0)
{
if (minWeights <= 0)
{
WarningIn
(
"metisDecomp::decompose"
"(const pointField&, const scalarField&)"
) << "Illegal minimum weight " << minWeights
<< endl;
}
if (cWeights.size() != mesh_.nCells())
{
FatalErrorIn
@ -441,11 +451,12 @@ Foam::labelList Foam::parMetisDecomp::decompose
<< " does not equal number of cells " << mesh_.nCells()
<< exit(FatalError);
}
// Convert to integers.
cellWeights.setSize(cWeights.size());
forAll(cellWeights, i)
{
cellWeights[i] = int(1000000*cWeights[i]/maxWeights);
cellWeights[i] = int(cWeights[i]/minWeights);
}
}
@ -779,9 +790,19 @@ Foam::labelList Foam::parMetisDecomp::decompose
// Check for externally provided cellweights and if so initialise weights
scalar maxWeights = gMax(cWeights);
scalar minWeights = gMin(cWeights);
if (cWeights.size() > 0)
{
if (minWeights <= 0)
{
WarningIn
(
"parMetisDecomp::decompose(const labelListList&"
", const pointField&, const scalarField&)"
) << "Illegal minimum weight " << minWeights
<< endl;
}
if (cWeights.size() != globalCellCells.size())
{
FatalErrorIn
@ -792,11 +813,12 @@ Foam::labelList Foam::parMetisDecomp::decompose
<< " does not equal number of cells " << globalCellCells.size()
<< exit(FatalError);
}
// Convert to integers.
cellWeights.setSize(cWeights.size());
forAll(cellWeights, i)
{
cellWeights[i] = int(1000000*cWeights[i]/maxWeights);
cellWeights[i] = int(cWeights[i]/minWeights);
}
}

5
src/parallel/parMetisDecomp/parMetisDecomp.H
View File

@ -112,6 +112,9 @@ public:
//- Return for every coordinate the wanted processor number. Use the
// mesh connectivity (if needed)
// Weights get normalised so the minimum value is 1 before truncation
// to an integer so the weights should be multiples of the minimum
// value. The overall sum of weights might otherwise overflow.
virtual labelList decompose
(
const pointField& points,
@ -122,6 +125,7 @@ public:
// passed agglomeration map (from fine to coarse cells) and coarse cell
// location. Can be overridden by decomposers that provide this
// functionality natively.
// See note on weights above.
virtual labelList decompose
(
const labelList& cellToRegion,
@ -136,6 +140,7 @@ public:
// from 0 at processor0 and then incrementing all through the
// processors)
// - the connections are across coupled patches
// See note on weights above.
virtual labelList decompose
(
const labelListList& globalCellCells,

31
src/postProcessing/functionObjects/forces/Make/files
View File

@ -4,10 +4,33 @@ forces/forcesFunctionObject.C
forceCoeffs/forceCoeffs.C
forceCoeffs/forceCoeffsFunctionObject.C
sDoFRBM = pointPatchFields/derived/sixDoFRigidBodyMotion
$(sDoFRBM)/sixDoFRigidBodyMotion.C
$(sDoFRBM)/sixDoFRigidBodyMotionIO.C
$(sDoFRBM)/sixDoFRigidBodyMotionState.C
$(sDoFRBM)/sixDoFRigidBodyMotionStateIO.C
sDoFRBMR = $(sDoFRBM)/sixDoFRigidBodyMotionRestraint
$(sDoFRBMR)/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C
$(sDoFRBMR)/sixDoFRigidBodyMotionRestraint/newSixDoFRigidBodyMotionRestraint.C
$(sDoFRBMR)/linearAxialAngularSpring/linearAxialAngularSpring.C
$(sDoFRBMR)/linearSpring/linearSpring.C
$(sDoFRBMR)/sphericalAngularSpring/sphericalAngularSpring.C
$(sDoFRBMR)/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C
sDoFRBMC = $(sDoFRBM)/sixDoFRigidBodyMotionConstraint
$(sDoFRBMC)/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C
$(sDoFRBMC)/sixDoFRigidBodyMotionConstraint/newSixDoFRigidBodyMotionConstraint.C
$(sDoFRBMC)/fixedAxis/fixedAxis.C
$(sDoFRBMC)/fixedLine/fixedLine.C
$(sDoFRBMC)/fixedOrientation/fixedOrientation.C
$(sDoFRBMC)/fixedPlane/fixedPlane.C
$(sDoFRBMC)/fixedPoint/fixedPoint.C
pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C
pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C
pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
LIB = $(FOAM_LIBBIN)/libforces

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

@ -160,7 +160,27 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
t.deltaTValue()
);
Field<vector>::operator=(motion_.generatePositions(p0_) - p0_);
// ----------------------------------------
// vector rotationAxis(0, 1, 0);
// vector torque
// (
// (
// (fm.second().first() + fm.second().second())
// & rotationAxis
// )
// *rotationAxis
// );
// motion_.updateForce
// (
// vector::zero, // Force no centre of mass motion
// torque, // Only rotation allowed around the unit rotationAxis
// t.deltaTValue()
// );
// ----------------------------------------
Field<vector>::operator=(motion_.currentPosition(p0_) - p0_);
fixedValuePointPatchField<vector>::updateCoeffs();
}

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

@ -1,210 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2009-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "sixDoFRigidBodyMotion.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
:
motionState_(),
refCentreOfMass_(vector::zero),
momentOfInertia_(diagTensor::one*VSMALL),
mass_(VSMALL)
{}
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
(
const point& centreOfMass,
const tensor& Q,
const vector& v,
const vector& a,
const vector& pi,
const vector& tau,
scalar mass,
const point& refCentreOfMass,
const diagTensor& momentOfInertia
)
:
motionState_
(
centreOfMass,
Q,
v,
a,
pi,
tau
),
refCentreOfMass_(refCentreOfMass),
momentOfInertia_(momentOfInertia),
mass_(mass)
{}
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
:
motionState_(dict),
refCentreOfMass_(dict.lookupOrDefault("refCentreOfMass", centreOfMass())),
momentOfInertia_(dict.lookup("momentOfInertia")),
mass_(readScalar(dict.lookup("mass")))
{}
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
(
const sixDoFRigidBodyMotion& sDoFRBM
)
:
motionState_(sDoFRBM.motionState()),
refCentreOfMass_(sDoFRBM.refCentreOfMass()),
momentOfInertia_(sDoFRBM.momentOfInertia()),
mass_(sDoFRBM.mass())
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotion::~sixDoFRigidBodyMotion()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::sixDoFRigidBodyMotion::updatePosition
(
scalar deltaT
)
{
// First leapfrog velocity adjust and motion part, required before
// force calculation
if (Pstream::master())
{
v() += 0.5*deltaT*a();
pi() += 0.5*deltaT*tau();
// Leapfrog move part
centreOfMass() += deltaT*v();
// Leapfrog orientation adjustment
tensor R;
R = rotationTensorX(0.5*deltaT*pi().x()/momentOfInertia_.xx());
pi() = pi() & R;
Q() = Q() & R;
R = rotationTensorY(0.5*deltaT*pi().y()/momentOfInertia_.yy());
pi() = pi() & R;
Q() = Q() & R;
R = rotationTensorZ(deltaT*pi().z()/momentOfInertia_.zz());
pi() = pi() & R;
Q() = Q() & R;
R = rotationTensorY(0.5*deltaT*pi().y()/momentOfInertia_.yy());
pi() = pi() & R;
Q() = Q() & R;
R = rotationTensorX(0.5*deltaT*pi().x()/momentOfInertia_.xx());
pi() = pi() & R;
Q() = Q() & R;
}
Pstream::scatter(motionState_);
}
void Foam::sixDoFRigidBodyMotion::updateForce
(
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT
)
{
// Second leapfrog velocity adjust part, required after motion and
// force calculation part
if (Pstream::master())
{
a() = fGlobal/mass_;
tau() = (Q().T() & tauGlobal);
v() += 0.5*deltaT*a();
pi() += 0.5*deltaT*tau();
}
Pstream::scatter(motionState_);
}
void Foam::sixDoFRigidBodyMotion::updateForce
(
const pointField& positions,
const vectorField& forces,
scalar deltaT
)
{
// Second leapfrog velocity adjust part, required after motion and
// force calculation part
if (Pstream::master())
{
a() = vector::zero;
tau() = vector::zero;
forAll(positions, i)
{
const vector& f = forces[i];
a() += f/mass_;
tau() += (positions[i] ^ (Q().T() & f));
}
v() += 0.5*deltaT*a();
pi() += 0.5*deltaT*tau();
}
Pstream::scatter(motionState_);
}
Foam::tmp<Foam::pointField>
Foam::sixDoFRigidBodyMotion::generatePositions(const pointField& pts) const
{
return (centreOfMass() + (Q() & (pts - refCentreOfMass_)));
}
// ************************************************************************* //

438
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C Normal file
View File

@ -0,0 +1,438 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2009-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "sixDoFRigidBodyMotion.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::sixDoFRigidBodyMotion::applyRestraints()
{
if (Pstream::master())
{
forAll(restraints_, rI)
{
// restraint position
point rP = vector::zero;
// restraint force
vector rF = vector::zero;
// restraint moment
vector rM = vector::zero;
restraints_[rI].restrain(*this, rP, rF, rM);
a() += rF/mass_;
// Moments are returned in global axes, transforming to
// body local to add to torque.
tau() += Q().T() & (rM + ((rP - centreOfMass()) ^ rF));
}
}
}
void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
{
if (Pstream::master())
{
label iter = 0;
bool allConverged = true;
// constraint force accumulator
vector cFA = vector::zero;
// constraint moment accumulator
vector cMA = vector::zero;
do
{
allConverged = true;
forAll(constraints_, cI)
{
// constraint position
point cP = vector::zero;
// constraint force
vector cF = vector::zero;
// constraint moment
vector cM = vector::zero;
bool constraintConverged = constraints_[cI].constrain
(
*this,
cFA,
cMA,
deltaT,
cP,
cF,
cM
);
allConverged = allConverged && constraintConverged;
// Accumulate constraint force
cFA += cF;
// Accumulate constraint moment
cMA += cM + ((cP - centreOfMass()) ^ cF);
}
} while(++iter < maxConstraintIters_ && !allConverged);
if (iter >= maxConstraintIters_)
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotion::applyConstraints"
"(scalar deltaT)"
)
<< nl << "Maximum number of sixDoFRigidBodyMotion constraint "
<< "iterations (" << maxConstraintIters_ << ") exceeded." << nl
<< exit(FatalError);
}
else if (report_)
{
Info<< "sixDoFRigidBodyMotion constraints converged in "
<< iter << " iterations"
<< nl << "Constraint force: " << cFA << nl
<< "Constraint moment: " << cMA
<< endl;
}
// Add the constraint forces and moments to the motion state variables
a() += cFA/mass_;
// The moment of constraint forces has already been added
// during accumulation. Moments are returned in global axes,
// transforming to body local
tau() += Q().T() & cMA;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
:
motionState_(),
restraints_(),
constraints_(),
maxConstraintIters_(0),
refCentreOfMass_(vector::zero),
momentOfInertia_(diagTensor::one*VSMALL),
mass_(VSMALL),
report_(false)
{}
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
(
const point& centreOfMass,
const tensor& Q,
const vector& v,
const vector& a,
const vector& pi,
const vector& tau,
scalar mass,
const point& refCentreOfMass,
const diagTensor& momentOfInertia,
bool report
)
:
motionState_
(
centreOfMass,
Q,
v,
a,
pi,
tau
),
restraints_(),
constraints_(),
maxConstraintIters_(0),
refCentreOfMass_(refCentreOfMass),
momentOfInertia_(momentOfInertia),
mass_(mass),
report_(report)
{}
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
:
motionState_(dict),
restraints_(),
constraints_(),
maxConstraintIters_(0),
refCentreOfMass_(dict.lookupOrDefault("refCentreOfMass", centreOfMass())),
momentOfInertia_(dict.lookup("momentOfInertia")),
mass_(readScalar(dict.lookup("mass"))),
report_(dict.lookupOrDefault<Switch>("report", false))
{
addRestraints(dict);
addConstraints(dict);
}
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
(
const sixDoFRigidBodyMotion& sDoFRBM
)
:
motionState_(sDoFRBM.motionState()),
restraints_(sDoFRBM.restraints()),
constraints_(sDoFRBM.constraints()),
maxConstraintIters_(sDoFRBM.maxConstraintIters()),
refCentreOfMass_(sDoFRBM.refCentreOfMass()),
momentOfInertia_(sDoFRBM.momentOfInertia()),
mass_(sDoFRBM.mass()),
report_(sDoFRBM.report())
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotion::~sixDoFRigidBodyMotion()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::sixDoFRigidBodyMotion::addRestraints
(
const dictionary& dict
)
{
if (dict.found("restraints"))
{
const dictionary& restraintDict = dict.subDict("restraints");
label i = 0;
restraints_.setSize(restraintDict.size());
forAllConstIter(IDLList<entry>, restraintDict, iter)
{
if (iter().isDict())
{
Info<< "Adding restraint: " << iter().keyword() << endl;
restraints_.set
(
i,
sixDoFRigidBodyMotionRestraint::New(iter().dict())
);
}
i++;
}
restraints_.setSize(i);
}
}
void Foam::sixDoFRigidBodyMotion::addConstraints
(
const dictionary& dict
)
{
if (dict.found("constraints"))
{
const dictionary& constraintDict = dict.subDict("constraints");
label i = 0;
constraints_.setSize(constraintDict.size());
forAllConstIter(IDLList<entry>, constraintDict, iter)
{
if (iter().isDict())
{
Info<< "Adding constraint: " << iter().keyword() << endl;
constraints_.set
(
i++,
sixDoFRigidBodyMotionConstraint::New(iter().dict())
);
}
}
constraints_.setSize(i);
if (constraints_.size())
{
maxConstraintIters_ = readLabel
(
constraintDict.lookup("maxIterations")
);
}
}
}
void Foam::sixDoFRigidBodyMotion::updatePosition
(
scalar deltaT
)
{
// First leapfrog velocity adjust and motion part, required before
// force calculation
if (Pstream::master())
{
v() += 0.5*deltaT*a();
pi() += 0.5*deltaT*tau();
// Leapfrog move part
centreOfMass() += deltaT*v();
// Leapfrog orientation adjustment
rotate(Q(), pi(), deltaT);
}
Pstream::scatter(motionState_);
}
void Foam::sixDoFRigidBodyMotion::updateForce
(
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT
)
{
// Second leapfrog velocity adjust part, required after motion and
// force calculation
if (Pstream::master())
{
a() = fGlobal/mass_;
tau() = (Q().T() & tauGlobal);
applyRestraints();
applyConstraints(deltaT);
v() += 0.5*deltaT*a();
pi() += 0.5*deltaT*tau();
if(report_)
{
status();
}
}
Pstream::scatter(motionState_);
}
void Foam::sixDoFRigidBodyMotion::updateForce
(
const pointField& positions,
const vectorField& forces,
scalar deltaT
)
{
vector a = vector::zero;
vector tau = vector::zero;
if (Pstream::master())
{
forAll(positions, i)
{
const vector& f = forces[i];
a += f/mass_;
tau += Q().T() & ((positions[i] - centreOfMass()) ^ f);
}
}
updateForce(a, tau, deltaT);
}
Foam::point Foam::sixDoFRigidBodyMotion::predictedPosition
(
const point& pt,
const vector& deltaForce,
const vector& deltaMoment,
scalar deltaT
) 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);
return (centreOfMassTemp + (QTemp & (pt - refCentreOfMass_)));
}
Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation
(
const vector& v,
const vector& deltaMoment,
scalar deltaT
) const
{
vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
tensor QTemp = Q();
rotate(QTemp, piTemp, deltaT);
return (QTemp & v);
}
void Foam::sixDoFRigidBodyMotion::status() const
{
Info<< "Centre of mass: " << centreOfMass() << nl
<< "Linear velocity: " << v() << nl
<< "Angular velocity: " << omega()
<< endl;
}
// ************************************************************************* //

207
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H → src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
View File

@ -26,10 +26,10 @@ Class
Foam::sixDoFRigidBodyMotion
Description
Six degree of freedom motion for a rigid body. Angular momentum stored in
body fixed reference frame. Reference orientation of the body must align
with the cartesian axes such that the Inertia tensor is in principle
component form.
Six degree of freedom motion for a rigid body. Angular momentum
stored in body fixed reference frame. Reference orientation of
the body (where Q = I) must align with the cartesian axes such
that the Inertia tensor is in principle component form.
Symplectic motion as per:
@ -43,6 +43,9 @@ Description
url = {http://link.aip.org/link/?JCP/107/5840/1},
doi = {10.1063/1.474310}
Can add restraints (i.e. a spring) and constraints (i.e. motion
may only be on a plane).
SourceFiles
sixDoFRigidBodyMotionI.H
sixDoFRigidBodyMotion.C
@ -55,6 +58,9 @@ SourceFiles
#include "sixDoFRigidBodyMotionState.H"
#include "pointField.H"
#include "sixDoFRigidBodyMotionRestraint.H"
#include "sixDoFRigidBodyMotionConstraint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -79,9 +85,19 @@ class sixDoFRigidBodyMotion
{
// Private data
// state data object
//- Motion state data object
sixDoFRigidBodyMotionState motionState_;
//- Restraints on the motion
PtrList<sixDoFRigidBodyMotionRestraint> restraints_;
//- Constaints on the motion
PtrList<sixDoFRigidBodyMotionConstraint> constraints_;
//- Maximum number of iterations allowed to attempt to obey
// constraints
label maxConstraintIters_;
//- Centre of mass of reference state
point refCentreOfMass_;
@ -91,6 +107,9 @@ class sixDoFRigidBodyMotion
//- Mass of the body
scalar mass_;
//- Switch to turn reporting of motion data on and off
Switch report_;
// Private Member Functions
@ -106,6 +125,74 @@ 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 the restraints to the object
void applyRestraints();
//- 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
// Access
//- Return access to the motion state
inline const sixDoFRigidBodyMotionState& motionState() const;
//- Return access to the restraints
inline const PtrList<sixDoFRigidBodyMotionRestraint>&
restraints() const;
//- Return access to the constraints
inline const PtrList<sixDoFRigidBodyMotionConstraint>&
constraints() const;
//- Return access to the maximum allowed number of
// constraint iterations
inline label maxConstraintIters() const;
//- Return access to the centre of mass
inline const point& refCentreOfMass() const;
//- Return access to the orientation
inline const tensor& Q() const;
//- Return access to velocity
inline const vector& v() const;
//- Return access to acceleration
inline const vector& a() const;
//- Return access to angular momentum
inline const vector& pi() const;
//- Return access to torque
inline const vector& tau() const;
// Edit
//- Return access to the centre of mass
inline point& refCentreOfMass();
//- Return non-const access to the orientation
inline tensor& Q();
//- Return non-const access to vector
inline vector& v();
//- Return non-const access to acceleration
inline vector& a();
//- Return non-const access to angular momentum
inline vector& pi();
//- Return non-const access to torque
inline vector& tau();
public:
@ -125,7 +212,8 @@ public:
const vector& tau,
scalar mass,
const point& refCentreOfMass,
const diagTensor& momentOfInertia
const diagTensor& momentOfInertia,
bool report = false
);
//- Construct from dictionary
@ -141,11 +229,23 @@ public:
// Member Functions
//- Add restraints to the motion, public to allow external
// addition of restraints after construction
void addRestraints(const dictionary& dict);
//- Add restraints to the motion, public to allow external
// addition of restraints after construction
void addConstraints(const dictionary& dict);
//- First leapfrog velocity adjust and motion part, required
// before force calculation
void updatePosition
(
scalar deltaT
);
//- Second leapfrog velocity adjust part, required after motion and
// force calculation
void updateForce
(
const vector& fGlobal,
@ -153,6 +253,8 @@ public:
scalar deltaT
);
//- Global forces supplied at locations, calculating net force
// and moment
void updateForce
(
const pointField& positions,
@ -160,39 +262,68 @@ public:
scalar deltaT
);
tmp<pointField> generatePositions(const pointField& pts) const;
//- Transform the given reference state pointField by the
// current motion state
inline tmp<pointField> currentPosition(const pointField& refPts) const;
//- Transform the given reference state point by the current
// motion state
inline point currentPosition(const point& refPt) const;
//- Transform the given reference state direction by the current
// motion state
inline vector currentOrientation(const vector& refDir) const;
//- Access the orientation tensor, Q.
// globalVector = Q & bodyLocalVector
// bodyLocalVector = Q.T() & globalVector
inline const tensor& currentOrientation() const;
//- Predict the position of the supplied point after deltaT
// given the current motion state and the additional supplied
// force and moment
point predictedPosition
(
const point& pt,
const vector& deltaForce,
const vector& deltaMoment,
scalar deltaT
) const;
//- Predict the orientation of the supplied vector after deltaT
// given the current motion state and the additional supplied
// moment
vector predictedOrientation
(
const vector& v,
const vector& deltaMoment,
scalar deltaT
) const;
//- Return the angular velocity in the global frame
inline vector omega() const;
//- Return the velocity of a position given by the current
// motion state
inline point currentVelocity(const point& pt) const;
//- Report the status of the motion
void status() const;
// Access
//- Return access to the motion state
inline const sixDoFRigidBodyMotionState& motionState() const;
//- Return access to the centre of mass
//- Return const access to the centre of mass
inline const point& centreOfMass() const;
//- Return access to the centre of mass
inline const point& refCentreOfMass() const;
//- Return access to the inertia tensor
inline const diagTensor& momentOfInertia() const;
//- Return access to the mass
//- Return const access to the mass
inline scalar mass() const;
//- Return access to the orientation
inline const tensor& Q() const;
//- Return access to velocity
inline const vector& v() const;
//- Return access to acceleration
inline const vector& a() const;
//- Return access to angular momentum
inline const vector& pi() const;
//- Return access to torque
inline const vector& tau() const;
//- Return the report Switch
inline bool report() const;
// Edit
@ -200,30 +331,12 @@ public:
//- Return non-const access to the centre of mass
inline point& centreOfMass();
//- Return access to the centre of mass
inline point& refCentreOfMass();
//- Return non-const access to the inertia tensor
inline diagTensor& momentOfInertia();
//- Return non-const access to the mass
inline scalar& mass();
//- Return non-const access to the orientation
inline tensor& Q();
//- Return non-const access to vector
inline vector& v();
//- Return non-const access to acceleration
inline vector& a();
//- Return non-const access to angular momentum
inline vector& pi();
//- Return non-const access to torque
inline vector& tau();
//- Write
void write(Ostream&) const;

178
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.C Normal file
View File

@ -0,0 +1,178 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "fixedAxis.H"
#include "addToRunTimeSelectionTable.H"
#include "sixDoFRigidBodyMotion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionConstraints
{
defineTypeNameAndDebug(fixedAxis, 0);
addToRunTimeSelectionTable
(
sixDoFRigidBodyMotionConstraint,
fixedAxis,
dictionary
);
};
};
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::fixedAxis
(
const dictionary& sDoFRBMCDict
)
:
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
fixedAxis_()
{
read(sDoFRBMCDict);
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::~fixedAxis()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const
{
constraintMomentIncrement = vector::zero;
vector predictedDir = motion.predictedOrientation
(
fixedAxis_,
existingConstraintMoment,
deltaT
);
scalar theta = acos(min(predictedDir & fixedAxis_, 1.0));
vector rotationAxis = fixedAxis_ ^ predictedDir;
scalar magRotationAxis = mag(rotationAxis);
if (magRotationAxis > VSMALL)
{
rotationAxis /= magRotationAxis;
const tensor& Q = motion.currentOrientation();
// Transform rotationAxis to body local system
rotationAxis = Q.T() & rotationAxis;
constraintMomentIncrement =
-relaxationFactor_
*(motion.momentOfInertia() & rotationAxis)
*theta/sqr(deltaT);
// Transform moment increment to global system
constraintMomentIncrement = Q & constraintMomentIncrement;
// Remove any moment that is around the fixedAxis_
constraintMomentIncrement -=
(constraintMomentIncrement & fixedAxis_)*fixedAxis_;
}
constraintPosition = motion.centreOfMass();
constraintForceIncrement = vector::zero;
bool converged(mag(theta) < tolerance_);
if (motion.report())
{
Info<< "Constraint " << this->name()
<< " angle " << theta
<< " force " << constraintForceIncrement
<< " moment " << constraintMomentIncrement;
if (converged)
{
Info<< " converged";
}
else
{
Info<< " not converged";
}
Info<< endl;
}
return converged;
}
bool Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::read
(
const dictionary& sDoFRBMCDict
)
{
sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
sDoFRBMCCoeffs_.lookup("axis") >> fixedAxis_;
scalar magFixedAxis(mag(fixedAxis_));
if (magFixedAxis > VSMALL)
{
fixedAxis_ /= magFixedAxis;
}
else
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::read"
"("
"const dictionary& sDoFRBMCDict"
")"
)
<< "axis has zero length"
<< abort(FatalError);
}
return true;
}
// ************************************************************************* //

126
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.H Normal file
View File

@ -0,0 +1,126 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::sixDoFRigidBodyMotionConstraints::fixedAxis
Description
sixDoFRigidBodyMotionConstraint. Axis of body fixed global
space.
SourceFiles
fixedAxis.C
\*---------------------------------------------------------------------------*/
#ifndef fixedAxis_H
#define fixedAxis_H
#include "sixDoFRigidBodyMotionConstraint.H"
#include "point.H"
#include "tensor.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionConstraints
{
/*---------------------------------------------------------------------------*\
Class fixedAxis Declaration
\*---------------------------------------------------------------------------*/
class fixedAxis
:
public sixDoFRigidBodyMotionConstraint
{
// Private data
//- Reference axis in global space
vector fixedAxis_;
public:
//- Runtime type information
TypeName("fixedAxis");
// Constructors
//- Construct from components
fixedAxis
(
const dictionary& sDoFRBMCDict
);
//- Construct and return a clone
virtual autoPtr<sixDoFRigidBodyMotionConstraint> clone() const
{
return autoPtr<sixDoFRigidBodyMotionConstraint>
(
new fixedAxis(*this)
);
}
// Destructor
virtual ~fixedAxis();
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCCoeff);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solidBodyMotionFunctions
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

166
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.C Normal file
View File

@ -0,0 +1,166 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "fixedLine.H"
#include "addToRunTimeSelectionTable.H"
#include "sixDoFRigidBodyMotion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionConstraints
{
defineTypeNameAndDebug(fixedLine, 0);
addToRunTimeSelectionTable
(
sixDoFRigidBodyMotionConstraint,
fixedLine,
dictionary
);
};
};
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::fixedLine::fixedLine
(
const dictionary& sDoFRBMCDict
)
:
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
refPt_(),
dir_()
{
read(sDoFRBMCDict);
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::fixedLine::~fixedLine()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionConstraints::fixedLine::constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const
{
point predictedPosition = motion.predictedPosition
(
refPt_,
existingConstraintForce,
existingConstraintMoment,
deltaT
);
constraintPosition = motion.currentPosition(refPt_);
// Info<< "current position " << constraintPosition << nl
// << "next predictedPosition " << predictedPosition
// << endl;
// Vector from reference point to predicted point
vector rC = predictedPosition - refPt_;
vector error = rC - ((rC) & dir_)*dir_;
// Info<< "error " << error << endl;
constraintForceIncrement =
-relaxationFactor_*error*motion.mass()/sqr(deltaT);
constraintMomentIncrement = vector::zero;
bool converged(mag(error) < tolerance_);
if (motion.report())
{
Info<< "Constraint " << this->name()
<< " error << " << error
<< " force " << constraintForceIncrement
<< " moment " << constraintMomentIncrement;
if (converged)
{
Info<< " converged";
}
else
{
Info<< " not converged";
}
Info<< endl;
}
return converged;
}
bool Foam::sixDoFRigidBodyMotionConstraints::fixedLine::read
(
const dictionary& sDoFRBMCDict
)
{
sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
sDoFRBMCCoeffs_.lookup("refPoint") >> refPt_;
sDoFRBMCCoeffs_.lookup("direction") >> dir_;
scalar magDir(mag(dir_));
if (magDir > VSMALL)
{
dir_ /= magDir;
}
else
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotionConstraints::fixedLine::read"
"("
"const dictionary& sDoFRBMCDict"
")"
)
<< "line direction has zero length"
<< abort(FatalError);
}
return true;
}
// ************************************************************************* //

127
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.H Normal file
View File

@ -0,0 +1,127 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::sixDoFRigidBodyMotionConstraints::fixedLine
Description
sixDoFRigidBodyMotionConstraint. Reference point may only move
along a line.
SourceFiles
fixedLine.C
\*---------------------------------------------------------------------------*/
#ifndef fixedLine_H
#define fixedLine_H
#include "sixDoFRigidBodyMotionConstraint.H"
#include "point.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionConstraints
{
/*---------------------------------------------------------------------------*\
Class fixedLine Declaration
\*---------------------------------------------------------------------------*/
class fixedLine
:
public sixDoFRigidBodyMotionConstraint
{
// Private data
//- Reference point for the constraining line
point refPt_;
//- Direction of the constraining line
vector dir_;
public:
//- Runtime type information
TypeName("fixedLine");
// Constructors
//- Construct from components
fixedLine
(
const dictionary& sDoFRBMCDict
);
//- Construct and return a clone
virtual autoPtr<sixDoFRigidBodyMotionConstraint> clone() const
{
return autoPtr<sixDoFRigidBodyMotionConstraint>
(
new fixedLine(*this)
);
}
// Destructor
virtual ~fixedLine();
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCCoeff);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solidBodyMotionFunctions
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

201
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.C Normal file
View File

@ -0,0 +1,201 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "fixedOrientation.H"
#include "addToRunTimeSelectionTable.H"
#include "sixDoFRigidBodyMotion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionConstraints
{
defineTypeNameAndDebug(fixedOrientation, 0);
addToRunTimeSelectionTable
(
sixDoFRigidBodyMotionConstraint,
fixedOrientation,
dictionary
);
};
};
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::fixedOrientation
(
const dictionary& sDoFRBMCDict
)
:
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
fixedOrientation_()
{
read(sDoFRBMCDict);
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::~fixedOrientation()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const
{
constraintMomentIncrement = vector::zero;
scalar maxTheta = -SMALL;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
vector axis = vector::zero;
axis[cmpt] = 1;
vector refDir = vector::zero;
refDir[(cmpt + 1) % 3] = 1;
vector predictedDir = motion.predictedOrientation
(
refDir,
existingConstraintMoment,
deltaT
);
axis = (fixedOrientation_ & axis);
refDir = (fixedOrientation_ & refDir);
// Removing any axis component from predictedDir
predictedDir -= (axis & predictedDir)*axis;
scalar theta = GREAT;
scalar magPredictedDir = mag(predictedDir);
if (magPredictedDir > VSMALL)
{
predictedDir /= magPredictedDir;
theta = acos(min(predictedDir & refDir, 1.0));
// Recalculating axis to give correct sign to angle
axis = (refDir ^ predictedDir);
scalar magAxis = mag(axis);
if (magAxis > VSMALL)
{
axis /= magAxis;
}
else
{
axis = vector::zero;
}
}
if (theta > maxTheta)
{
maxTheta = theta;
}
constraintMomentIncrement +=
-relaxationFactor_
*theta*axis
*motion.momentOfInertia()[cmpt]/sqr(deltaT);
}
constraintPosition = motion.centreOfMass();
constraintForceIncrement = vector::zero;
bool converged(mag(maxTheta) < tolerance_);
if (motion.report())
{
Info<< "Constraint " << this->name()
<< " max angle " << maxTheta
<< " force " << constraintForceIncrement
<< " moment " << constraintMomentIncrement;
if (converged)
{
Info<< " converged";
}
else
{
Info<< " not converged";
}
Info<< endl;
}
return converged;
}
bool Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::read
(
const dictionary& sDoFRBMCDict
)
{
sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
fixedOrientation_ =
sDoFRBMCCoeffs_.lookupOrDefault<tensor>("fixedOrientation", I);
if (mag(mag(fixedOrientation_) - sqrt(3.0)) > 1e-9)
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::read"
"("
"const dictionary& sDoFRBMCDict"
")"
)
<< "fixedOrientation " << fixedOrientation_
<< " is not a rotation tensor."
<< exit(FatalError);
}
return true;
}
// ************************************************************************* //

127
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.H Normal file
View File

@ -0,0 +1,127 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation
Description
sixDoFRigidBodyMotionConstraint. Orientation of body fixed global
space. Only valid where the predicted deviation from alignment is
< 90 degrees.
SourceFiles
fixedOrientation.C
\*---------------------------------------------------------------------------*/
#ifndef fixedOrientation_H
#define fixedOrientation_H
#include "sixDoFRigidBodyMotionConstraint.H"
#include "point.H"
#include "tensor.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionConstraints
{
/*---------------------------------------------------------------------------*\
Class fixedOrientation Declaration
\*---------------------------------------------------------------------------*/
class fixedOrientation
:
public sixDoFRigidBodyMotionConstraint
{
// Private data
//- Reference orientation where there is no moment
tensor fixedOrientation_;
public:
//- Runtime type information
TypeName("fixedOrientation");
// Constructors
//- Construct from components
fixedOrientation
(
const dictionary& sDoFRBMCDict
);
//- Construct and return a clone
virtual autoPtr<sixDoFRigidBodyMotionConstraint> clone() const
{
return autoPtr<sixDoFRigidBodyMotionConstraint>
(
new fixedOrientation(*this)
);
}
// Destructor
virtual ~fixedOrientation();
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCCoeff);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solidBodyMotionFunctions
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

149
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.C Normal file
View File

@ -0,0 +1,149 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "fixedPlane.H"
#include "addToRunTimeSelectionTable.H"
#include "sixDoFRigidBodyMotion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionConstraints
{
defineTypeNameAndDebug(fixedPlane, 0);
addToRunTimeSelectionTable
(
sixDoFRigidBodyMotionConstraint,
fixedPlane,
dictionary
);
};
};
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::fixedPlane
(
const dictionary& sDoFRBMCDict
)
:
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
fixedPlane_(vector::one)
{
read(sDoFRBMCDict);
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::~fixedPlane()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const
{
const point& refPt = fixedPlane_.refPoint();
const vector& n = fixedPlane_.normal();
point predictedPosition = motion.predictedPosition
(
refPt,
existingConstraintForce,
existingConstraintMoment,
deltaT
);
constraintPosition = motion.currentPosition(refPt);
// Info<< "current position " << constraintPosition << nl
// << "next predictedPosition " << predictedPosition
// << endl;
vector error = ((predictedPosition - refPt) & n)*n;
// Info<< "error " << error << endl;
constraintForceIncrement =
-relaxationFactor_*error*motion.mass()/sqr(deltaT);
constraintMomentIncrement = vector::zero;
bool converged(mag(error) < tolerance_);
if (motion.report())
{
Info<< "Constraint " << this->name()
<< " error " << error
<< " force " << constraintForceIncrement
<< " moment " << constraintMomentIncrement;
if (converged)
{
Info<< " converged";
}
else
{
Info<< " not converged";
}
Info<< endl;
}
return converged;
}
bool Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::read
(
const dictionary& sDoFRBMCDict
)
{
sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
point refPt = sDoFRBMCCoeffs_.lookup("refPoint");
vector normal = sDoFRBMCCoeffs_.lookup("normal");
fixedPlane_ = plane(refPt, normal);
return true;
}
// ************************************************************************* //

125
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.H Normal file
View File

@ -0,0 +1,125 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::sixDoFRigidBodyMotionConstraints::fixedPlane
Description
sixDoFRigidBodyMotionConstraint. Reference point may only move
along a plane.
SourceFiles
fixedPlane.C
\*---------------------------------------------------------------------------*/
#ifndef fixedPlane_H
#define fixedPlane_H
#include "sixDoFRigidBodyMotionConstraint.H"
#include "plane.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionConstraints
{
/*---------------------------------------------------------------------------*\
Class fixedPlane Declaration
\*---------------------------------------------------------------------------*/
class fixedPlane
:
public sixDoFRigidBodyMotionConstraint
{
// Private data
//- Plane which the transformed reference point is constrained
// to move along
plane fixedPlane_;
public:
//- Runtime type information
TypeName("fixedPlane");
// Constructors
//- Construct from components
fixedPlane
(
const dictionary& sDoFRBMCDict
);
//- Construct and return a clone
virtual autoPtr<sixDoFRigidBodyMotionConstraint> clone() const
{
return autoPtr<sixDoFRigidBodyMotionConstraint>
(
new fixedPlane(*this)
);
}
// Destructor
virtual ~fixedPlane();
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCCoeff);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solidBodyMotionFunctions
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

154
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.C Normal file
View File

@ -0,0 +1,154 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "fixedPoint.H"
#include "addToRunTimeSelectionTable.H"
#include "sixDoFRigidBodyMotion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionConstraints
{
defineTypeNameAndDebug(fixedPoint, 0);
addToRunTimeSelectionTable
(
sixDoFRigidBodyMotionConstraint,
fixedPoint,
dictionary
);
};
};
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::fixedPoint
(
const dictionary& sDoFRBMCDict
)
:
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
fixedPoint_()
{
read(sDoFRBMCDict);
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::~fixedPoint()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const
{
point predictedPosition = motion.predictedPosition
(
fixedPoint_,
existingConstraintForce,
existingConstraintMoment,
deltaT
);
constraintPosition = motion.currentPosition(fixedPoint_);
// Info<< "current position " << constraintPosition << nl
// << "next predictedPosition " << predictedPosition
// << endl;
vector error = predictedPosition - fixedPoint_;
// Info<< "error " << error << endl;
// Correction force derived from Lagrange multiplier:
// G = -lambda*grad(sigma)
// where
// sigma = mag(error) = 0
// so
// grad(sigma) = error/mag(error)
// Solving for lambda using the SHAKE methodology gives
// lambda = mass*mag(error)/sqr(deltaT)
// This is only strictly applicable (i.e. will converge in one
// iteration) to constraints at the centre of mass. Everything
// else will need to iterate, and may need under-relaxed to be
// stable.
constraintForceIncrement =
-relaxationFactor_*error*motion.mass()/sqr(deltaT);
constraintMomentIncrement = vector::zero;
bool converged(mag(error) < tolerance_);
if (motion.report())
{
Info<< "Constraint " << this->name()
<< " error " << error
<< " force " << constraintForceIncrement
<< " moment " << constraintMomentIncrement;
if (converged)
{
Info<< " converged";
}
else
{
Info<< " not converged";
}
Info<< endl;
}
return converged;
}
bool Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::read
(
const dictionary& sDoFRBMCDict
)
{
sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
sDoFRBMCCoeffs_.lookup("fixedPoint") >> fixedPoint_;
return true;
}
// ************************************************************************* //

126
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.H Normal file
View File

@ -0,0 +1,126 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::sixDoFRigidBodyMotionConstraints::fixedPoint
Description
sixDoFRigidBodyMotionConstraint. Point fixed in space.
SourceFiles
fixedPoint.C
\*---------------------------------------------------------------------------*/
#ifndef fixedPoint_H
#define fixedPoint_H
#include "sixDoFRigidBodyMotionConstraint.H"
#include "point.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionConstraints
{
/*---------------------------------------------------------------------------*\
Class fixedPoint Declaration
\*---------------------------------------------------------------------------*/
class fixedPoint
:
public sixDoFRigidBodyMotionConstraint
{
// Private data
//- Point which is rigidly attached to the body and constrains
// it so that this point remains fixed. This serves as the
// reference point for displacements as well as the target
// position.
point fixedPoint_;
public:
//- Runtime type information
TypeName("fixedPoint");
// Constructors
//- Construct from components
fixedPoint
(
const dictionary& sDoFRBMCDict
);
//- Construct and return a clone
virtual autoPtr<sixDoFRigidBodyMotionConstraint> clone() const
{
return autoPtr<sixDoFRigidBodyMotionConstraint>
(
new fixedPoint(*this)
);
}
// Destructor
virtual ~fixedPoint();
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCCoeff);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solidBodyMotionFunctions
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

65
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/newSixDoFRigidBodyMotionConstraint.C Normal file
View File

@ -0,0 +1,65 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "sixDoFRigidBodyMotionConstraint.H"
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::sixDoFRigidBodyMotionConstraint>
Foam::sixDoFRigidBodyMotionConstraint::New(const dictionary& sDoFRBMCDict)
{
word sixDoFRigidBodyMotionConstraintTypeName =
sDoFRBMCDict.lookup("sixDoFRigidBodyMotionConstraint");
Info<< "Selecting sixDoFRigidBodyMotionConstraint function "
<< sixDoFRigidBodyMotionConstraintTypeName << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find
(
sixDoFRigidBodyMotionConstraintTypeName
);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorIn
(
"sixDoFRigidBodyMotionConstraint::New"
"("
"const dictionary& sDoFRBMCDict"
")"
) << "Unknown sixDoFRigidBodyMotionConstraint type "
<< sixDoFRigidBodyMotionConstraintTypeName << endl << endl
<< "Valid sixDoFRigidBodyMotionConstraints are : " << endl
<< dictionaryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<sixDoFRigidBodyMotionConstraint>(cstrIter()(sDoFRBMCDict));
}
// ************************************************************************* //

86
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C Normal file
View File

@ -0,0 +1,86 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "sixDoFRigidBodyMotionConstraint.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::sixDoFRigidBodyMotionConstraint, 0);
defineRunTimeSelectionTable(Foam::sixDoFRigidBodyMotionConstraint, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraint::sixDoFRigidBodyMotionConstraint
(
const dictionary& sDoFRBMCDict
)
:
name_(fileName(sDoFRBMCDict.name().name()).components(token::COLON).last()),
sDoFRBMCCoeffs_
(
sDoFRBMCDict.subDict
(
word(sDoFRBMCDict.lookup("sixDoFRigidBodyMotionConstraint"))
+ "Coeffs"
)
),
tolerance_(readScalar(sDoFRBMCDict.lookup("tolerance"))),
relaxationFactor_
(
sDoFRBMCDict.lookupOrDefault<scalar>("relaxationFactor", 1)
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraint::~sixDoFRigidBodyMotionConstraint()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionConstraint::read
(
const dictionary& sDoFRBMCDict
)
{
tolerance_ = (readScalar(sDoFRBMCDict.lookup("tolerance")));
relaxationFactor_ = sDoFRBMCDict.lookupOrDefault<scalar>
(
"relaxationFactor",
1
);
sDoFRBMCCoeffs_ = sDoFRBMCDict.subDict(type() + "Coeffs");
return true;
}
// ************************************************************************* //

186
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H Normal file
View File

@ -0,0 +1,186 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Namespace
Foam::sixDoFRigidBodyMotionConstraints
Description
Namespace for six DoF motion constraints
Class
Foam::sixDoFRigidBodyMotionConstraint
Description
Base class for defining constraints for sixDoF motions
SourceFiles
sixDoFRigidBodyMotionConstraint.C
newDynamicFvMesh.C
\*---------------------------------------------------------------------------*/
#ifndef sixDoFRigidBodyMotionConstraint_H
#define sixDoFRigidBodyMotionConstraint_H
#include "Time.H"
#include "dictionary.H"
#include "autoPtr.H"
#include "vector.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class sixDoFRigidBodyMotion;
/*---------------------------------------------------------------------------*\
Class sixDoFRigidBodyMotionConstraint Declaration
\*---------------------------------------------------------------------------*/
class sixDoFRigidBodyMotionConstraint
{
protected:
// Protected data
//- Name of the constraint in dictionary
word name_;
//- Constraint model specific coefficient dictionary
dictionary sDoFRBMCCoeffs_;
//- Solution tolerance. Meaning depends on model, usually an
// absolute distance or angle.
scalar tolerance_;
//- Relaxation factor for solution, default to one
scalar relaxationFactor_;
public:
//- Runtime type information
TypeName("sixDoFRigidBodyMotionConstraint");
// Declare run-time constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
sixDoFRigidBodyMotionConstraint,
dictionary,
(const dictionary& sDoFRBMCDict),
(sDoFRBMCDict)
);
// Constructors
//- Construct from the sDoFRBMCDict dictionary and Time
sixDoFRigidBodyMotionConstraint
(
const dictionary& sDoFRBMCDict
);
//- Construct and return a clone
virtual autoPtr<sixDoFRigidBodyMotionConstraint> clone() const = 0;
// Selectors
//- Select constructed from the sDoFRBMCDict dictionary and Time
static autoPtr<sixDoFRigidBodyMotionConstraint> New
(
const dictionary& sDoFRBMCDict
);
// Destructor
virtual ~sixDoFRigidBodyMotionConstraint();
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const = 0;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCDict) = 0;
// Access
//- Return access to the name of the restraint
inline const word& name() const
{
return name_;
}
// Return access to sDoFRBMCCoeffs
inline const dictionary& coeffDict() const
{
return sDoFRBMCCoeffs_;
}
//- Return access to the tolerance
inline scalar tolerance() const
{
return tolerance_;
}
//- Return access to the relaxationFactor
inline scalar relaxationFactor() const
{
return relaxationFactor_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

173
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H → src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
View File

@ -24,8 +24,6 @@ License
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline Foam::tensor
@ -64,6 +62,37 @@ Foam::sixDoFRigidBodyMotion::rotationTensorZ(scalar phi) const
}
inline void Foam::sixDoFRigidBodyMotion::rotate
(
tensor& Q,
vector& pi,
scalar deltaT
) const
{
tensor R;
R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
pi = pi & R;
Q = Q & R;
R = rotationTensorY(0.5*deltaT*pi.y()/momentOfInertia_.yy());
pi = pi & R;
Q = Q & R;
R = rotationTensorZ(deltaT*pi.z()/momentOfInertia_.zz());
pi = pi & R;
Q = Q & R;
R = rotationTensorY(0.5*deltaT*pi.y()/momentOfInertia_.yy());
pi = pi & R;
Q = Q & R;
R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
pi = pi & R;
Q = Q & R;
}
inline const Foam::sixDoFRigidBodyMotionState&
Foam::sixDoFRigidBodyMotion::motionState() const
{
@ -71,9 +100,23 @@ Foam::sixDoFRigidBodyMotion::motionState() const
}
inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const
inline const Foam::PtrList<Foam::sixDoFRigidBodyMotionRestraint>&
Foam::sixDoFRigidBodyMotion::restraints() const
{
return motionState_.centreOfMass();
return restraints_;
}
inline const Foam::PtrList<Foam::sixDoFRigidBodyMotionConstraint>&
Foam::sixDoFRigidBodyMotion::constraints() const
{
return constraints_;
}
inline Foam::label Foam::sixDoFRigidBodyMotion::maxConstraintIters() const
{
return maxConstraintIters_;
}
@ -83,19 +126,6 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::refCentreOfMass() const
}
inline const Foam::diagTensor&
Foam::sixDoFRigidBodyMotion::momentOfInertia() const
{
return momentOfInertia_;
}
inline Foam::scalar Foam::sixDoFRigidBodyMotion::mass() const
{
return mass_;
}
inline const Foam::tensor& Foam::sixDoFRigidBodyMotion::Q() const
{
return motionState_.Q();
@ -126,30 +156,12 @@ inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau() const
}
inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass()
{
return motionState_.centreOfMass();
}
inline Foam::point& Foam::sixDoFRigidBodyMotion::refCentreOfMass()
{
return refCentreOfMass_;
}
inline Foam::diagTensor& Foam::sixDoFRigidBodyMotion::momentOfInertia()
{
return momentOfInertia_;
}
inline Foam::scalar& Foam::sixDoFRigidBodyMotion::mass()
{
return mass_;
}
inline Foam::tensor& Foam::sixDoFRigidBodyMotion::Q()
{
return motionState_.Q();
@ -180,4 +192,95 @@ inline Foam::vector& Foam::sixDoFRigidBodyMotion::tau()
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
inline Foam::tmp<Foam::pointField>
Foam::sixDoFRigidBodyMotion::currentPosition(const pointField& refPts) const
{
return (centreOfMass() + (Q() & (refPts - refCentreOfMass_)));
}
inline Foam::point Foam::sixDoFRigidBodyMotion::currentPosition
(
const point& refPt
) const
{
return (centreOfMass() + (Q() & (refPt - refCentreOfMass_)));
}
inline Foam::vector Foam::sixDoFRigidBodyMotion::currentOrientation
(
const vector& refDir
) const
{
return (Q() & refDir);
}
inline const Foam::tensor&
Foam::sixDoFRigidBodyMotion::currentOrientation() const
{
return Q();
}
inline Foam::vector Foam::sixDoFRigidBodyMotion::omega() const
{
return Q() & (inv(momentOfInertia_) & pi());
}
inline Foam::point Foam::sixDoFRigidBodyMotion::currentVelocity
(
const point& pt
) const
{
return (omega() ^ (pt - centreOfMass())) + v();
}
inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const
{
return motionState_.centreOfMass();
}
inline const Foam::diagTensor&
Foam::sixDoFRigidBodyMotion::momentOfInertia() const
{
return momentOfInertia_;
}
inline Foam::scalar Foam::sixDoFRigidBodyMotion::mass() const
{
return mass_;
}
inline bool Foam::sixDoFRigidBodyMotion::report() const
{
return report_;
}
inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass()
{
return motionState_.centreOfMass();
}
inline Foam::diagTensor& Foam::sixDoFRigidBodyMotion::momentOfInertia()
{
return momentOfInertia_;
}
inline Foam::scalar& Foam::sixDoFRigidBodyMotion::mass()
{
return mass_;
}
// ************************************************************************* //

68
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C → src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
View File

@ -39,6 +39,72 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
<< momentOfInertia_ << token::END_STATEMENT << nl;
os.writeKeyword("mass")
<< mass_ << token::END_STATEMENT << nl;
os.writeKeyword("report")
<< report_ << token::END_STATEMENT << nl;
if (restraints_.size())
{
dictionary restraintsDict;
forAll(restraints_, rI)
{
word restraintType = restraints_[rI].type();
dictionary restraintDict;
restraintDict.add("sixDoFRigidBodyMotionRestraint", restraintType);
restraintDict.add
(
word(restraintType + "Coeffs"), restraints_[rI].coeffDict()
);
restraintsDict.add(restraints_[rI].name(), restraintDict);
}
os.writeKeyword("restraints") << restraintsDict;
}
if (constraints_.size())
{
dictionary constraintsDict;
constraintsDict.add("maxIterations", maxConstraintIters_);
forAll(constraints_, rI)
{
word constraintType = constraints_[rI].type();
dictionary constraintDict;
constraintDict.add
(
"sixDoFRigidBodyMotionConstraint",
constraintType
);
constraintDict.add
(
"tolerance",
constraints_[rI].tolerance()
);
constraintDict.add
(
"relaxationFactor",
constraints_[rI].relaxationFactor()
);
constraintDict.add
(
word(constraintType + "Coeffs"), constraints_[rI].coeffDict()
);
constraintsDict.add(constraints_[rI].name(), constraintDict);
}
os.writeKeyword("constraints") << constraintsDict;
}
}
@ -71,7 +137,7 @@ Foam::Ostream& Foam::operator<<
os << sDoFRBM.motionState()
<< token::SPACE << sDoFRBM.refCentreOfMass()
<< token::SPACE << sDoFRBM.momentOfInertia()
<< token::SPACE << sDoFRBM.mass() ;
<< token::SPACE << sDoFRBM.mass();
// Check state of Ostream
os.check

202
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C Normal file
View File

@ -0,0 +1,202 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "linearAxialAngularSpring.H"
#include "addToRunTimeSelectionTable.H"
#include "sixDoFRigidBodyMotion.H"
#include "transform.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionRestraints
{
defineTypeNameAndDebug(linearAxialAngularSpring, 0);
addToRunTimeSelectionTable
(
sixDoFRigidBodyMotionRestraint,
linearAxialAngularSpring,
dictionary
);
};
};
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::
linearAxialAngularSpring
(
const dictionary& sDoFRBMRDict
)
:
sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
refQ_(),
axis_(),
stiffness_(),
damping_()
{
read(sDoFRBMRDict);
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::
~linearAxialAngularSpring()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void
Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::restrain
(
const sixDoFRigidBodyMotion& motion,
vector& restraintPosition,
vector& restraintForce,
vector& restraintMoment
) const
{
vector refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 1, 0);
vector oldDir = refQ_ & refDir;
vector newDir = motion.currentOrientation(refDir);
if (mag(oldDir & axis_) > 0.95 || mag(newDir & axis_) > 0.95)
{
// Directions getting close to the axis, change reference
refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 0, 1);
vector oldDir = refQ_ & refDir;
vector newDir = motion.currentOrientation(refDir);
}
// Removing any axis component from oldDir and newDir and normalising
oldDir -= (axis_ & oldDir)*axis_;
oldDir /= mag(oldDir);
newDir -= (axis_ & newDir)*axis_;
newDir /= mag(newDir);
scalar theta = mag(acos(min(oldDir & newDir, 1.0)));
// Temporary axis with sign information.
vector a = (oldDir ^ newDir);
// Remove any component that is not along axis that may creep in
a = (a & axis_)*axis_;
scalar magA = mag(a);
if (magA > VSMALL)
{
a /= magA;
}
else
{
a = vector::zero;
}
// Damping of along axis angular velocity only
restraintMoment = -stiffness_*theta*a - damping_*(motion.omega() & a)*a;
restraintForce = vector::zero;
// Not needed to be altered as restraintForce is zero, but set to
// centreOfMass to be sure of no spurious moment
restraintPosition = motion.centreOfMass();
if (motion.report())
{
Info<< "Restraint " << this->name()
<< " angle " << theta
<< " force " << restraintForce
<< " moment " << restraintMoment
<< endl;
}
}
bool Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::read
(
const dictionary& sDoFRBMRDict
)
{
sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict);
refQ_ = sDoFRBMRCoeffs_.lookupOrDefault<tensor>("referenceOrientation", I);
if (mag(mag(refQ_) - sqrt(3.0)) > 1e-9)
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::"
"read"
"("
"const dictionary& sDoFRBMRDict"
")"
)
<< "referenceOrientation " << refQ_ << " is not a rotation tensor. "
<< "mag(referenceOrientation) - sqrt(3) = "
<< mag(refQ_) - sqrt(3.0) << nl
<< exit(FatalError);
}
axis_ = sDoFRBMRCoeffs_.lookup("axis");
scalar magAxis(mag(axis_));
if (magAxis > VSMALL)
{
axis_ /= magAxis;
}
else
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::"
"read"
"("
"const dictionary& sDoFRBMCDict"
")"
)
<< "axis has zero length"
<< abort(FatalError);
}
sDoFRBMRCoeffs_.lookup("stiffness") >> stiffness_;
sDoFRBMRCoeffs_.lookup("damping") >> damping_;
return true;
}
// ************************************************************************* //

129
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.H Normal file
View File

@ -0,0 +1,129 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring
Description
sixDoFRigidBodyMotionRestraints model. Linear axial angular spring.
SourceFiles
linearAxialAngularSpring.C
\*---------------------------------------------------------------------------*/
#ifndef linearAxialAngularSpring_H
#define linearAxialAngularSpring_H
#include "sixDoFRigidBodyMotionRestraint.H"
#include "point.H"
#include "tensor.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionRestraints
{
/*---------------------------------------------------------------------------*\
Class linearAxialAngularSpring Declaration
\*---------------------------------------------------------------------------*/
class linearAxialAngularSpring
:
public sixDoFRigidBodyMotionRestraint
{
// Private data
//- Reference orientation where there is no moment
tensor refQ_;
//- Global unit axis around which the motion is sprung
vector axis_;
//- Spring stiffness coefficient (Nm/rad)
scalar stiffness_;
//- Damping coefficient (Nms/rad)
scalar damping_;
public:
//- Runtime type information
TypeName("linearAxialAngularSpring");
// Constructors
//- Construct from components
linearAxialAngularSpring
(
const dictionary& sDoFRBMRDict
);
//- Construct and return a clone
virtual autoPtr<sixDoFRigidBodyMotionRestraint> clone() const
{
return autoPtr<sixDoFRigidBodyMotionRestraint>
(
new linearAxialAngularSpring(*this)
);
}
// Destructor
virtual ~linearAxialAngularSpring();
// Member Functions
//- Calculate the restraint position, force and moment.
// Global reference frame vectors.
virtual void restrain
(
const sixDoFRigidBodyMotion& motion,
vector& restraintPosition,
vector& restraintForce,
vector& restraintMoment
) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMRCoeff);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solidBodyMotionFunctions
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

128
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.C Normal file
View File

@ -0,0 +1,128 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "linearSpring.H"
#include "addToRunTimeSelectionTable.H"
#include "sixDoFRigidBodyMotion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionRestraints
{
defineTypeNameAndDebug(linearSpring, 0);
addToRunTimeSelectionTable
(
sixDoFRigidBodyMotionRestraint,
linearSpring,
dictionary
);
};
};
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionRestraints::linearSpring::linearSpring
(
const dictionary& sDoFRBMRDict
)
:
sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
anchor_(),
refAttachmentPt_(),
stiffness_(),
damping_(),
restLength_()
{
read(sDoFRBMRDict);
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionRestraints::linearSpring::~linearSpring()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::sixDoFRigidBodyMotionRestraints::linearSpring::restrain
(
const sixDoFRigidBodyMotion& motion,
vector& restraintPosition,
vector& restraintForce,
vector& restraintMoment
) const
{
restraintPosition = motion.currentPosition(refAttachmentPt_);
vector r = restraintPosition - anchor_;
scalar magR = mag(r);
// r is now the r unit vector
r /= magR;
vector v = motion.currentVelocity(restraintPosition);
restraintForce = -stiffness_*(magR - restLength_)*r - damping_*(r & v)*r;
restraintMoment = vector::zero;
if (motion.report())
{
Info<< "Restraint " << this->name()
<< " spring length " << magR
<< " force " << restraintForce
<< " moment " << restraintMoment
<< endl;
}
}
bool Foam::sixDoFRigidBodyMotionRestraints::linearSpring::read
(
const dictionary& sDoFRBMRDict
)
{
sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict);
sDoFRBMRCoeffs_.lookup("anchor") >> anchor_;
sDoFRBMRCoeffs_.lookup("refAttachmentPt") >> refAttachmentPt_;
sDoFRBMRCoeffs_.lookup("stiffness") >> stiffness_;
sDoFRBMRCoeffs_.lookup("damping") >> damping_;
sDoFRBMRCoeffs_.lookup("restLength") >> restLength_;
return true;
}
// ************************************************************************* //

132
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.H Normal file
View File

@ -0,0 +1,132 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::sixDoFRigidBodyMotionRestraints::linearSpring
Description
sixDoFRigidBodyMotionRestraints model. Linear spring.
SourceFiles
linearSpring.C
\*---------------------------------------------------------------------------*/
#ifndef linearSpring_H
#define linearSpring_H
#include "sixDoFRigidBodyMotionRestraint.H"
#include "point.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionRestraints
{
/*---------------------------------------------------------------------------*\
Class linearSpring Declaration
\*---------------------------------------------------------------------------*/
class linearSpring
:
public sixDoFRigidBodyMotionRestraint
{
// Private data
//- Anchor point, where the spring is attached to an immovable
// object
point anchor_;
//- Reference point of attachment to the solid body
point refAttachmentPt_;
//- Spring stiffness coefficient (N/m)
scalar stiffness_;
//- Damping coefficient (Ns/m)
scalar damping_;
//- Rest length - length of spring when no forces are applied to it
scalar restLength_;
public:
//- Runtime type information
TypeName("linearSpring");
// Constructors
//- Construct from components
linearSpring
(
const dictionary& sDoFRBMRDict
);
//- Construct and return a clone
virtual autoPtr<sixDoFRigidBodyMotionRestraint> clone() const
{
return autoPtr<sixDoFRigidBodyMotionRestraint>
(
new linearSpring(*this)
);
}
// Destructor
virtual ~linearSpring();
// Member Functions
//- Calculate the restraint position, force and moment.
// Global reference frame vectors.
virtual void restrain
(
const sixDoFRigidBodyMotion& motion,
vector& restraintPosition,
vector& restraintForce,
vector& restraintMoment
) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMRCoeff);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solidBodyMotionFunctions
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

43
src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFields.H → src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/newSixDoFRigidBodyMotionRestraint.C
View File

@ -24,27 +24,42 @@ License
\*---------------------------------------------------------------------------*/
#ifndef expDirectionMixedFvPatchFields_H
#define expDirectionMixedFvPatchFields_H
#include "sixDoFRigidBodyMotionRestraint.H"
#include "expDirectionMixedFvPatchField.H"
#include "fieldTypes.H"
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
Foam::autoPtr<Foam::sixDoFRigidBodyMotionRestraint>
Foam::sixDoFRigidBodyMotionRestraint::New(const dictionary& sDoFRBMRDict)
{
word sixDoFRigidBodyMotionRestraintTypeName =
sDoFRBMRDict.lookup("sixDoFRigidBodyMotionRestraint");
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "Selecting sixDoFRigidBodyMotionRestraint function "
<< sixDoFRigidBodyMotionRestraintTypeName << endl;
makePatchTypeFieldTypedefs(expDirectionMixed)
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find
(
sixDoFRigidBodyMotionRestraintTypeName
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorIn
(
"sixDoFRigidBodyMotionRestraint::New"
"("
"const dictionary& sDoFRBMRDict"
")"
) << "Unknown sixDoFRigidBodyMotionRestraint type "
<< sixDoFRigidBodyMotionRestraintTypeName << endl << endl
<< "Valid sixDoFRigidBodyMotionRestraints are : " << endl
<< dictionaryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
} // End namespace Foam
return autoPtr<sixDoFRigidBodyMotionRestraint>(cstrIter()(sDoFRBMRDict));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

53
src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFields.C → src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C
View File

@ -24,25 +24,50 @@ License
\*---------------------------------------------------------------------------*/
#include "expDirectionMixedFvPatchFields.H"
#include "addToRunTimeSelectionTable.H"
#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
#include "sixDoFRigidBodyMotionRestraint.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template<>
makePatchTypeField(fvPatchVectorField, expDirectionMixedFvPatchVectorField);
defineTypeNameAndDebug(Foam::sixDoFRigidBodyMotionRestraint, 0);
template<>
makePatchTypeField(fvPatchTensorField, expDirectionMixedFvPatchTensorField);
defineRunTimeSelectionTable(Foam::sixDoFRigidBodyMotionRestraint, dictionary);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionRestraint::sixDoFRigidBodyMotionRestraint
(
const dictionary& sDoFRBMRDict
)
:
name_(fileName(sDoFRBMRDict.name().name()).components(token::COLON).last()),
sDoFRBMRCoeffs_
(
sDoFRBMRDict.subDict
(
word(sDoFRBMRDict.lookup("sixDoFRigidBodyMotionRestraint"))
+ "Coeffs"
)
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionRestraint::~sixDoFRigidBodyMotionRestraint()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionRestraint::read
(
const dictionary& sDoFRBMRDict
)
{
sDoFRBMRCoeffs_ = sDoFRBMRDict.subDict(type() + "Coeffs");
return true;
}
} // End namespace Foam
// ************************************************************************* //

163
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H Normal file
View File

@ -0,0 +1,163 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Namespace
Foam::sixDoFRigidBodyMotionRestraints
Description
Namespace for six DoF motion restraints
Class
Foam::sixDoFRigidBodyMotionRestraint
Description
Base class for defining restraints for sixDoF motions
SourceFiles
sixDoFRigidBodyMotionRestraint.C
newDynamicFvMesh.C
\*---------------------------------------------------------------------------*/
#ifndef sixDoFRigidBodyMotionRestraint_H
#define sixDoFRigidBodyMotionRestraint_H
#include "Time.H"
#include "dictionary.H"
#include "autoPtr.H"
#include "vector.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class sixDoFRigidBodyMotion;
/*---------------------------------------------------------------------------*\
Class sixDoFRigidBodyMotionRestraint Declaration
\*---------------------------------------------------------------------------*/
class sixDoFRigidBodyMotionRestraint
{
protected:
// Protected data
//- Name of the constraint in dictionary
word name_;
//- Restraint model specific coefficient dictionary
dictionary sDoFRBMRCoeffs_;
public:
//- Runtime type information
TypeName("sixDoFRigidBodyMotionRestraint");
// Declare run-time constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
sixDoFRigidBodyMotionRestraint,
dictionary,
(const dictionary& sDoFRBMRDict),
(sDoFRBMRDict)
);
// Constructors
//- Construct from the sDoFRBMRDict dictionary and Time
sixDoFRigidBodyMotionRestraint
(
const dictionary& sDoFRBMRDict
);
//- Construct and return a clone
virtual autoPtr<sixDoFRigidBodyMotionRestraint> clone() const = 0;
// Selectors
//- Select constructed from the sDoFRBMRDict dictionary and Time
static autoPtr<sixDoFRigidBodyMotionRestraint> New
(
const dictionary& sDoFRBMRDict
);
// Destructor
virtual ~sixDoFRigidBodyMotionRestraint();
// Member Functions
//- Calculate the restraint position, force and moment.
// Global reference frame vectors.
virtual void restrain
(
const sixDoFRigidBodyMotion& motion,
vector& restraintPosition,
vector& restraintForce,
vector& restraintMoment
) const = 0;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMRDict) = 0;
// Access
//- Return access to the name of the restraint
inline const word& name() const
{
return name_;
}
// Return access to sDoFRBMRCoeffs
inline const dictionary& coeffDict() const
{
return sDoFRBMRCoeffs_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

156
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.C Normal file
View File

@ -0,0 +1,156 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "sphericalAngularSpring.H"
#include "addToRunTimeSelectionTable.H"
#include "sixDoFRigidBodyMotion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionRestraints
{
defineTypeNameAndDebug(sphericalAngularSpring, 0);
addToRunTimeSelectionTable
(
sixDoFRigidBodyMotionRestraint,
sphericalAngularSpring,
dictionary
);
};
};
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::
sphericalAngularSpring
(
const dictionary& sDoFRBMRDict
)
:
sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
refQ_(),
stiffness_(),
damping_()
{
read(sDoFRBMRDict);
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::
~sphericalAngularSpring()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void
Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::restrain
(
const sixDoFRigidBodyMotion& motion,
vector& restraintPosition,
vector& restraintForce,
vector& restraintMoment
) const
{
restraintMoment = vector::zero;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
vector axis = vector::zero;
axis[cmpt] = 1;
vector refDir = vector::zero;
refDir[(cmpt + 1) % 3] = 1;
vector newDir = motion.currentOrientation(refDir);
axis = (refQ_ & axis);
refDir = (refQ_ & refDir);
newDir -= (axis & newDir)*axis;
restraintMoment += -stiffness_*(refDir ^ newDir);
}
restraintMoment += -damping_*motion.omega();
restraintForce = vector::zero;
// Not needed to be altered as restraintForce is zero, but set to
// centreOfMass to be sure of no spurious moment
restraintPosition = motion.centreOfMass();
if (motion.report())
{
Info<< "Restraint " << this->name()
<< " force " << restraintForce
<< " moment " << restraintMoment
<< endl;
}
}
bool Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::read
(
const dictionary& sDoFRBMRDict
)
{
sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict);
refQ_ = sDoFRBMRCoeffs_.lookupOrDefault<tensor>("referenceOrientation", I);
if (mag(mag(refQ_) - sqrt(3.0)) > 1e-9)
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::"
"read"
"("
"const dictionary& sDoFRBMRDict"
")"
)
<< "referenceOrientation " << refQ_ << " is not a rotation tensor. "
<< "mag(referenceOrientation) - sqrt(3) = "
<< mag(refQ_) - sqrt(3.0) << nl
<< exit(FatalError);
}
sDoFRBMRCoeffs_.lookup("stiffness") >> stiffness_;
sDoFRBMRCoeffs_.lookup("damping") >> damping_;
return true;
}
// ************************************************************************* //

126
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.H Normal file
View File

@ -0,0 +1,126 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring
Description
sixDoFRigidBodyMotionRestraints model. Spherical angular spring.
SourceFiles
sphericalAngularSpring.C
\*---------------------------------------------------------------------------*/
#ifndef sphericalAngularSpring_H
#define sphericalAngularSpring_H
#include "sixDoFRigidBodyMotionRestraint.H"
#include "point.H"
#include "tensor.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionRestraints
{
/*---------------------------------------------------------------------------*\
Class sphericalAngularSpring Declaration
\*---------------------------------------------------------------------------*/
class sphericalAngularSpring
:
public sixDoFRigidBodyMotionRestraint
{
// Private data
//- Reference orientation where there is no moment
tensor refQ_;
//- Spring stiffness coefficient (Nm/rad)
scalar stiffness_;
//- Damping coefficient (Nms/rad)
scalar damping_;
public:
//- Runtime type information
TypeName("sphericalAngularSpring");
// Constructors
//- Construct from components
sphericalAngularSpring
(
const dictionary& sDoFRBMRDict
);
//- Construct and return a clone
virtual autoPtr<sixDoFRigidBodyMotionRestraint> clone() const
{
return autoPtr<sixDoFRigidBodyMotionRestraint>
(
new sphericalAngularSpring(*this)
);
}
// Destructor
virtual ~sphericalAngularSpring();
// Member Functions
//- Calculate the restraint position, force and moment.
// Global reference frame vectors.
virtual void restrain
(
const sixDoFRigidBodyMotion& motion,
vector& restraintPosition,
vector& restraintForce,
vector& restraintMoment
) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMRCoeff);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solidBodyMotionFunctions
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

226
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C Normal file
View File

@ -0,0 +1,226 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "tabulatedAxialAngularSpring.H"
#include "addToRunTimeSelectionTable.H"
#include "sixDoFRigidBodyMotion.H"
#include "transform.H"
#include "unitConversion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionRestraints
{
defineTypeNameAndDebug(tabulatedAxialAngularSpring, 0);
addToRunTimeSelectionTable
(
sixDoFRigidBodyMotionRestraint,
tabulatedAxialAngularSpring,
dictionary
);
};
};
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::
tabulatedAxialAngularSpring
(
const dictionary& sDoFRBMRDict
)
:
sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
refQ_(),
axis_(),
moment_(),
convertToDegrees_(),
damping_()
{
read(sDoFRBMRDict);
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::
~tabulatedAxialAngularSpring()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void
Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::restrain
(
const sixDoFRigidBodyMotion& motion,
vector& restraintPosition,
vector& restraintForce,
vector& restraintMoment
) const
{
vector refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 1, 0);
vector oldDir = refQ_ & refDir;
vector newDir = motion.currentOrientation(refDir);
if (mag(oldDir & axis_) > 0.95 || mag(newDir & axis_) > 0.95)
{
// Directions getting close to the axis, change reference
refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 0, 1);
vector oldDir = refQ_ & refDir;
vector newDir = motion.currentOrientation(refDir);
}
// Removing any axis component from oldDir and newDir and normalising
oldDir -= (axis_ & oldDir)*axis_;
oldDir /= mag(oldDir);
newDir -= (axis_ & newDir)*axis_;
newDir /= mag(newDir);
scalar theta = mag(acos(min(oldDir & newDir, 1.0)));
// Determining the sign of theta by comparing the cross product of
// the direction vectors with the axis
theta *= sign((oldDir ^ newDir) & axis_);
scalar moment;
if (convertToDegrees_)
{
moment = moment_(radToDeg(theta));
}
else
{
moment = moment_(theta);
}
// Damping of along axis angular velocity only
restraintMoment = moment*axis_ - damping_*(motion.omega() & axis_)*axis_;
restraintForce = vector::zero;
// Not needed to be altered as restraintForce is zero, but set to
// centreOfMass to be sure of no spurious moment
restraintPosition = motion.centreOfMass();
if (motion.report())
{
Info<< "Restraint " << this->name()
<< " angle " << theta
<< " force " << restraintForce
<< " moment " << restraintMoment
<< endl;
}
}
bool Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::read
(
const dictionary& sDoFRBMRDict
)
{
sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict);
refQ_ = sDoFRBMRCoeffs_.lookupOrDefault<tensor>("referenceOrientation", I);
if (mag(mag(refQ_) - sqrt(3.0)) > 1e-9)
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotionRestraints::"
"tabulatedAxialAngularSpring::read"
"("
"const dictionary& sDoFRBMRDict"
")"
)
<< "referenceOrientation " << refQ_ << " is not a rotation tensor. "
<< "mag(referenceOrientation) - sqrt(3) = "
<< mag(refQ_) - sqrt(3.0) << nl
<< exit(FatalError);
}
axis_ = sDoFRBMRCoeffs_.lookup("axis");
scalar magAxis(mag(axis_));
if (magAxis > VSMALL)
{
axis_ /= magAxis;
}
else
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotionRestraints::"
"tabulatedAxialAngularSpring::read"
"("
"const dictionary& sDoFRBMCDict"
")"
)
<< "axis has zero length"
<< abort(FatalError);
}
moment_ = interpolationTable<scalar>(sDoFRBMRCoeffs_);
word angleFormat = sDoFRBMRCoeffs_.lookup("angleFormat");
if (angleFormat == "degrees" || angleFormat == "degree")
{
convertToDegrees_ = true;
}
else if (angleFormat == "radians" || angleFormat == "radian")
{
convertToDegrees_ = false;
}
else
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotionRestraints::"
"tabulatedAxialAngularSpring::read"
"("
"const dictionary& sDoFRBMCDict"
")"
)
<< "angleFormat must be degree, degrees, radian or radians"
<< abort(FatalError);
}
sDoFRBMRCoeffs_.lookup("damping") >> damping_;
return true;
}
// ************************************************************************* //

136
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H Normal file
View File

@ -0,0 +1,136 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring
Description
sixDoFRigidBodyMotionRestraints model. Axial angular spring with moment
values drawn from an interpolation table. Linear damping.
SourceFiles
tabulatedAxialAngularSpring.C
\*---------------------------------------------------------------------------*/
#ifndef tabulatedAxialAngularSpring_H
#define tabulatedAxialAngularSpring_H
#include "sixDoFRigidBodyMotionRestraint.H"
#include "point.H"
#include "tensor.H"
#include "interpolationTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFRigidBodyMotionRestraints
{
/*---------------------------------------------------------------------------*\
Class tabulatedAxialAngularSpring Declaration
\*---------------------------------------------------------------------------*/
class tabulatedAxialAngularSpring
:
public sixDoFRigidBodyMotionRestraint
{
// Private data
//- Reference orientation where there is no moment
tensor refQ_;
//- Global unit axis around which the motion is sprung
vector axis_;
//- Spring moment interpolation table, depending on angleFormat
interpolationTable<scalar> moment_;
//- Boolean stating whether the angle around the axis needs to
// be converted to degrees before asking the
// interpolationTable for a moment value
bool convertToDegrees_;
//- Damping coefficient (Nms/rad)
scalar damping_;
public:
//- Runtime type information
TypeName("tabulatedAxialAngularSpring");
// Constructors
//- Construct from components
tabulatedAxialAngularSpring
(
const dictionary& sDoFRBMRDict
);
//- Construct and return a clone
virtual autoPtr<sixDoFRigidBodyMotionRestraint> clone() const
{
return autoPtr<sixDoFRigidBodyMotionRestraint>
(
new tabulatedAxialAngularSpring(*this)
);
}
// Destructor
virtual ~tabulatedAxialAngularSpring();
// Member Functions
//- Calculate the restraint position, force and moment.
// Global reference frame vectors.
virtual void restrain
(
const sixDoFRigidBodyMotion& motion,
vector& restraintPosition,
vector& restraintForce,
vector& restraintMoment
) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMRCoeff);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solidBodyMotionFunctions
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

0
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C → src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C
View File

0
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.H → src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.H
View File

0
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H → src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H
View File

0
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C → src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C
View File

173
src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C Normal file
View File

@ -0,0 +1,173 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.H"
#include "pointPatchFields.H"
#include "addToRunTimeSelectionTable.H"
#include "Time.H"
#include "fvMesh.H"
#include "volFields.H"
#include "uniformDimensionedFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
(
const pointPatch& p,
const DimensionedField<vector, pointMesh>& iF
)
:
fixedValuePointPatchField<vector>(p, iF),
motion_(),
p0_(p.localPoints()),
rhoInf_(1.0)
{}
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
(
const pointPatch& p,
const DimensionedField<vector, pointMesh>& iF,
const dictionary& dict
)
:
fixedValuePointPatchField<vector>(p, iF, dict),
motion_(dict),
rhoInf_(readScalar(dict.lookup("rhoInf")))
{
if (!dict.found("value"))
{
updateCoeffs();
}
if (dict.found("p0"))
{
p0_ = vectorField("p0", dict , p.size());
}
else
{
p0_ = p.localPoints();
}
}
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
(
const uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField& ptf,
const pointPatch& p,
const DimensionedField<vector, pointMesh>& iF,
const pointPatchFieldMapper& mapper
)
:
fixedValuePointPatchField<vector>(ptf, p, iF, mapper),
motion_(ptf.motion_),
p0_(ptf.p0_, mapper),
rhoInf_(ptf.rhoInf_)
{}
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
(
const uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField& ptf,
const DimensionedField<vector, pointMesh>& iF
)
:
fixedValuePointPatchField<vector>(ptf, iF),
motion_(ptf.motion_),
p0_(ptf.p0_),
rhoInf_(ptf.rhoInf_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
{
if (this->updated())
{
return;
}
const polyMesh& mesh = this->dimensionedInternalField().mesh()();
const Time& t = mesh.time();
motion_.updatePosition(t.deltaTValue());
vector gravity = vector::zero;
if (db().foundObject<uniformDimensionedVectorField>("g"))
{
uniformDimensionedVectorField g =
db().lookupObject<uniformDimensionedVectorField>("g");
gravity = g.value();
}
// Do not modify the accelerations, except with gravity, where available
motion_.updateForce(gravity*motion_.mass(), vector::zero, t.deltaTValue());
Field<vector>::operator=(motion_.currentPosition(p0_) - p0_);
fixedValuePointPatchField<vector>::updateCoeffs();
}
void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::write
(
Ostream& os
) const
{
pointPatchField<vector>::write(os);
motion_.write(os);
os.writeKeyword("rhoInf")
<< rhoInf_ << token::END_STATEMENT << nl;
p0_.writeEntry("p0", os);
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePointPatchTypeField
(
pointPatchVectorField,
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

158
src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.H Normal file
View File

@ -0,0 +1,158 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
Description
Foam::uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
SourceFiles
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
\*---------------------------------------------------------------------------*/
#ifndef uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField_H
#define uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField_H
#include "fixedValuePointPatchField.H"
#include "sixDoFRigidBodyMotion.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField Declaration
\*---------------------------------------------------------------------------*/
class uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
:
public fixedValuePointPatchField<vector>
{
// Private data
//- Six dof motion object
sixDoFRigidBodyMotion motion_;
//- Reference positions of points on the patch
pointField p0_;
//- Reference density required by the forces object for
// incompressible calculations. Retained here to give
// dictionary compatibility with other sixDoF patches.
scalar rhoInf_;
public:
//- Runtime type information
TypeName("uncoupledSixDoFRigidBodyDisplacement");
// Constructors
//- Construct from patch and internal field
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
(
const pointPatch&,
const DimensionedField<vector, pointMesh>&
);
//- Construct from patch, internal field and dictionary
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
(
const pointPatch&,
const DimensionedField<vector, pointMesh>&,
const dictionary&
);
//- Construct by mapping given patchField<vector> onto a new patch
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
(
const uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField&,
const pointPatch&,
const DimensionedField<vector, pointMesh>&,
const pointPatchFieldMapper&
);
//- Construct and return a clone
virtual autoPtr<pointPatchField<vector> > clone() const
{
return autoPtr<pointPatchField<vector> >
(
new uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
(
*this
)
);
}
//- Construct as copy setting internal field reference
uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
(
const uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField&,
const DimensionedField<vector, pointMesh>&
);
//- Construct and return a clone setting internal field reference
virtual autoPtr<pointPatchField<vector> > clone
(
const DimensionedField<vector, pointMesh>& iF
) const
{
return autoPtr<pointPatchField<vector> >
(
new uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField
(
*this,
iF
)
);
}
// Member functions
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

20
src/randomProcesses/fft/kShellIntegration.C
View File

@ -29,12 +29,7 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
graph kShellIntegration
Foam::graph Foam::kShellIntegration
(
const complexVectorField& Ek,
const Kmesh& K
@ -55,8 +50,8 @@ graph kShellIntegration
// now scale this to get the energy in a box of side l0
scalar l0(K.sizeOfBox()[0]*(scalar(K.nn()[0])/(scalar(K.nn()[0]) - 1.0)));
scalar factor = pow3(l0/constant::mathematical::twoPi);
scalar l0(K.sizeOfBox()[0]*(scalar(K.nn()[0])/(scalar(K.nn()[0])-1.0)));
scalar factor = pow((l0/(2.0*constant::mathematical::pi)),3.0);
y *= factor;
@ -72,7 +67,7 @@ graph kShellIntegration
// kShellMean : average over the points in a k-shell to evaluate the
// radial part of the energy spectrum.
graph kShellMean
Foam::graph Foam::kShellMean
(
const complexVectorField& Ek,
const Kmesh& K
@ -136,13 +131,8 @@ graph kShellMean
}
}
return graph("E(k)", "k", "E(k)", Ek1D, k1D);
return graph("E(k)", "k", "E(k)", k1D, Ek1D);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

1
src/sampling/Make/files
View File

@ -38,6 +38,7 @@ sampledSurface/sampledSurface/sampledSurface.C
sampledSurface/sampledSurfaces/sampledSurfaces.C
sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C
sampledSurface/sampledSurfacesFunctionObject/sampledSurfacesFunctionObject.C
sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
sampledSurface/thresholdCellFaces/thresholdCellFaces.C
sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.C

352
src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C Normal file
View File

@ -0,0 +1,352 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "sampledTriSurfaceMesh.H"
#include "dictionary.H"
#include "polyMesh.H"
#include "polyPatch.H"
#include "volFields.H"
#include "meshSearch.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(sampledTriSurfaceMesh, 0);
addToRunTimeSelectionTable
(
sampledSurface,
sampledTriSurfaceMesh,
word
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
(
const word& name,
const polyMesh& mesh,
const word& surfaceName
)
:
sampledSurface(name, mesh),
surface_
(
IOobject
(
name,
mesh.time().constant(), // instance
"triSurface", // local
mesh, // registry
IOobject::MUST_READ,
IOobject::NO_WRITE
)
),
needsUpdate_(true),
cellLabels_(0),
pointMap_(0),
faceMap_(0)
{}
Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
(
const word& name,
const polyMesh& mesh,
const dictionary& dict
)
:
sampledSurface(name, mesh, dict),
surface_
(
IOobject
(
dict.lookup("surface"),
mesh.time().constant(), // instance
"triSurface", // local
mesh, // registry
IOobject::MUST_READ,
IOobject::NO_WRITE
)
),
needsUpdate_(true),
cellLabels_(0),
pointMap_(0),
faceMap_(0)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sampledTriSurfaceMesh::~sampledTriSurfaceMesh()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::sampledTriSurfaceMesh::needsUpdate() const
{
return needsUpdate_;
}
bool Foam::sampledTriSurfaceMesh::expire()
{
// already marked as expired
if (needsUpdate_)
{
return false;
}
sampledSurface::clearGeom();
MeshStorage::clear();
cellLabels_.clear();
pointMap_.clear();
faceMap_.clear();
needsUpdate_ = true;
return true;
}
bool Foam::sampledTriSurfaceMesh::update()
{
if (!needsUpdate_)
{
return false;
}
// Find the cells the triangles of the surface are in.
// Does approximation by looking at the face centres only
const pointField fc = surface_.faceCentres();
cellLabels_.setSize(fc.size());
meshSearch meshSearcher(mesh(), false);
forAll(fc, triI)
{
cellLabels_[triI] = meshSearcher.findCell
(
fc[triI],
-1, // seed
true // use tree
);
}
boolList include(surface_.size(), true);
label nNotFound = 0;
forAll(cellLabels_, triI)
{
if (cellLabels_[triI] == -1)
{
include[triI] = false;
nNotFound++;
}
}
label nTotalNotFound = returnReduce(nNotFound, sumOp<label>());
if (debug)
{
Pout<< "surface:" << surface_.size()
<< " included:" << surface_.size()-nNotFound
<< " total not found" << nTotalNotFound << endl;
}
// Add to master all triangles are outside all meshes.
{
boolList onAnyProc(surface_.size(), false);
Pstream::listCombineGather(onAnyProc, orEqOp<bool>());
if (Pstream::master())
{
label nAdded = 0;
forAll(onAnyProc, triI)
{
if (!onAnyProc[triI])
{
include[triI] = true;
nAdded++;
}
}
if (debug)
{
Pout<< "nAdded to master:" << nAdded << endl;
}
}
}
// Now subset the surface
triSurface localSurface
(
surface_.subsetMesh
(
include,
pointMap_,
faceMap_
)
);
// And convert into faces.
faceList& faces = this->storedFaces();
faces.setSize(localSurface.size());
forAll(localSurface, i)
{
faces[i] = localSurface[i].triFaceFace();
}
this->storedPoints() = localSurface.points();
if (debug)
{
print(Pout);
Pout << endl;
}
needsUpdate_ = false;
return true;
}
Foam::tmp<Foam::scalarField>
Foam::sampledTriSurfaceMesh::sample
(
const volScalarField& vField
) const
{
return sampleField(vField);
}
Foam::tmp<Foam::vectorField>
Foam::sampledTriSurfaceMesh::sample
(
const volVectorField& vField
) const
{
return sampleField(vField);
}
Foam::tmp<Foam::sphericalTensorField>
Foam::sampledTriSurfaceMesh::sample
(
const volSphericalTensorField& vField
) const
{
return sampleField(vField);
}
Foam::tmp<Foam::symmTensorField>
Foam::sampledTriSurfaceMesh::sample
(
const volSymmTensorField& vField
) const
{
return sampleField(vField);
}
Foam::tmp<Foam::tensorField>
Foam::sampledTriSurfaceMesh::sample
(
const volTensorField& vField
) const
{
return sampleField(vField);
}
Foam::tmp<Foam::scalarField>
Foam::sampledTriSurfaceMesh::interpolate
(
const interpolation<scalar>& interpolator
) const
{
return interpolateField(interpolator);
}
Foam::tmp<Foam::vectorField>
Foam::sampledTriSurfaceMesh::interpolate
(
const interpolation<vector>& interpolator
) const
{
return interpolateField(interpolator);
}
Foam::tmp<Foam::sphericalTensorField>
Foam::sampledTriSurfaceMesh::interpolate
(
const interpolation<sphericalTensor>& interpolator
) const
{
return interpolateField(interpolator);
}
Foam::tmp<Foam::symmTensorField>
Foam::sampledTriSurfaceMesh::interpolate
(
const interpolation<symmTensor>& interpolator
) const
{
return interpolateField(interpolator);
}
Foam::tmp<Foam::tensorField>
Foam::sampledTriSurfaceMesh::interpolate
(
const interpolation<tensor>& interpolator
) const
{
return interpolateField(interpolator);
}
void Foam::sampledTriSurfaceMesh::print(Ostream& os) const
{
os << "sampledTriSurfaceMesh: " << name() << " :"
<< " surface:" << surface_.objectRegistry::name()
<< " faces:" << faces().size()
<< " points:" << points().size();
}
// ************************************************************************* //

241
src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H Normal file
View File

@ -0,0 +1,241 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::sampledTriSurfaceMesh
Description
A sampledSurface from a triSurfaceMesh. It samples on the points/triangles
of the triSurface.
In parallel every processor just operates on the part of the surface
where the face centres are inside the mesh. It is then up to the
caller to stitch the partial surfaces together.
No check is done for triangle centres being on multiple processors
(should never happen but truncation errors ...)
Any value on a triangle outside the mesh will be set to 0.
SourceFiles
sampledTriSurfaceMesh.C
\*---------------------------------------------------------------------------*/
#ifndef sampledTriSurfaceMesh_H
#define sampledTriSurfaceMesh_H
#include "sampledSurface.H"
#include "triSurfaceMesh.H"
#include "MeshedSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class sampledTriSurfaceMesh Declaration
\*---------------------------------------------------------------------------*/
class sampledTriSurfaceMesh
:
public sampledSurface,
public MeshedSurface<face>
{
//- Private typedefs for convenience
typedef MeshedSurface<face> MeshStorage;
// Private data
//- Surface to sample on
const triSurfaceMesh surface_;
//- Track if the surface needs an update
mutable bool needsUpdate_;
//- Local cell labels
labelList cellLabels_;
//- From local surface back to surface_
labelList pointMap_;
//- From local surface back to surface_
labelList faceMap_;
// Private Member Functions
//- sample field on faces
template <class Type>
tmp<Field<Type> > sampleField
(
const GeometricField<Type, fvPatchField, volMesh>& vField
) const;
template <class Type>
tmp<Field<Type> >
interpolateField(const interpolation<Type>&) const;
public:
//- Runtime type information
TypeName("sampledTriSurfaceMesh");
// Constructors
//- Construct from components
sampledTriSurfaceMesh
(
const word& name,
const polyMesh& mesh,
const word& surfaceName
);
//- Construct from dictionary
sampledTriSurfaceMesh
(
const word& name,
const polyMesh& mesh,
const dictionary& dict
);
// Destructor
virtual ~sampledTriSurfaceMesh();
// Member Functions
//- Does the surface need an update?
virtual bool needsUpdate() const;
//- Mark the surface as needing an update.
// May also free up unneeded data.
// Return false if surface was already marked as expired.
virtual bool expire();
//- Update the surface as required.
// Do nothing (and return false) if no update was needed
virtual bool update();
//- Points of surface
virtual const pointField& points() const
{
return MeshStorage::points();
}
//- Faces of surface
virtual const faceList& faces() const
{
return MeshStorage::faces();
}
//- sample field on surface
virtual tmp<scalarField> sample
(
const volScalarField&
) const;
//- sample field on surface
virtual tmp<vectorField> sample
(
const volVectorField&
) const;
//- sample field on surface
virtual tmp<sphericalTensorField> sample
(
const volSphericalTensorField&
) const;
//- sample field on surface
virtual tmp<symmTensorField> sample
(
const volSymmTensorField&
) const;
//- sample field on surface
virtual tmp<tensorField> sample
(
const volTensorField&
) const;
//- interpolate field on surface
virtual tmp<scalarField> interpolate
(
const interpolation<scalar>&
) const;
//- interpolate field on surface
virtual tmp<vectorField> interpolate
(
const interpolation<vector>&
) const;
//- interpolate field on surface
virtual tmp<sphericalTensorField> interpolate
(
const interpolation<sphericalTensor>&
) const;
//- interpolate field on surface
virtual tmp<symmTensorField> interpolate
(
const interpolation<symmTensor>&
) const;
//- interpolate field on surface
virtual tmp<tensorField> interpolate
(
const interpolation<tensor>&
) const;
//- Write
virtual void print(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "sampledTriSurfaceMeshTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

91
src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C Normal file
View File

@ -0,0 +1,91 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "sampledTriSurfaceMesh.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template <class Type>
Foam::tmp<Foam::Field<Type> >
Foam::sampledTriSurfaceMesh::sampleField
(
const GeometricField<Type, fvPatchField, volMesh>& vField
) const
{
// One value per face
tmp<Field<Type> > tvalues(new Field<Type>(faceMap_.size()));
Field<Type>& values = tvalues();
forAll(faceMap_, i)
{
label cellI = cellLabels_[faceMap_[i]];
if (cellI != -1)
{
values[i] = vField[cellI];
}
else
{
values[i] = pTraits<Type>::zero;
}
}
return tvalues;
}
template <class Type>
Foam::tmp<Foam::Field<Type> >
Foam::sampledTriSurfaceMesh::interpolateField
(
const interpolation<Type>& interpolator
) const
{
// One value per vertex
tmp<Field<Type> > tvalues(new Field<Type>(pointMap_.size()));
Field<Type>& values = tvalues();
forAll(pointMap_, i)
{
label origPointI = pointMap_[i];
label origTriI = surface_.pointFaces()[origPointI][0];
label cellI = cellLabels_[origTriI];
if (cellI != -1)
{
values[i] = interpolator.interpolate(points()[i], cellI);
}
else
{
values[i] = pTraits<Type>::zero;
}
}
return tvalues;
}
// ************************************************************************* //

2
tutorials/discreteMethods/dsmcFoam/freeSpacePeriodic/constant/dsmcProperties
View File

@ -59,7 +59,6 @@ moleculeProperties
mass 46.5e-27;
diameter 4.17e-10;
internalDegreesOfFreedom 2;
viscosityCoefficient 1.656e-5;
omega 0.74;
}
@ -68,7 +67,6 @@ moleculeProperties
mass 53.12e-27;
diameter 4.07e-10;
internalDegreesOfFreedom 2;
viscosityCoefficient 1.919e-5;
omega 0.77;
}
}

2
tutorials/discreteMethods/dsmcFoam/freeSpaceStream/constant/dsmcProperties
View File

@ -69,7 +69,6 @@ moleculeProperties
mass 46.5e-27;
diameter 4.17e-10;
internalDegreesOfFreedom 2;
viscosityCoefficient 1.656e-5;
omega 0.74;
}
@ -78,7 +77,6 @@ moleculeProperties
mass 53.12e-27;
diameter 4.07e-10;
internalDegreesOfFreedom 2;
viscosityCoefficient 1.919e-5;
omega 0.77;
}
}

1
tutorials/discreteMethods/dsmcFoam/supersonicCorner/constant/dsmcProperties
View File

@ -80,7 +80,6 @@ moleculeProperties
mass 66.3e-27;
diameter 4.17e-10;
internalDegreesOfFreedom 0;
viscosityCoefficient 2.117e-5;
omega 0.81;
}
}

2
tutorials/discreteMethods/dsmcFoam/wedge15Ma5/constant/dsmcProperties
View File

@ -71,7 +71,6 @@ moleculeProperties
mass 46.5e-27;
diameter 4.17e-10;
internalDegreesOfFreedom 2;
viscosityCoefficient 1.656e-5;
omega 0.74;
}
@ -80,7 +79,6 @@ moleculeProperties
mass 53.12e-27;
diameter 4.07e-10;
internalDegreesOfFreedom 2;
viscosityCoefficient 1.919e-5;
omega 0.77;
}
}

17
tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
View File

@ -116,9 +116,9 @@ castellatedMeshControls
// Refinement parameters
// ~~~~~~~~~~~~~~~~~~~~~
// While refining maximum number of cells per processor. This is basically
// the number of cells that fit on a processor. If you choose this too small
// it will do just more refinement iterations to obtain a similar mesh.
// If local number of cells is >= maxLocalCells on any processor
// switches from from refinement followed by balancing
// (current method) to (weighted) balancing before refinement.
maxLocalCells 1000000;
// Overall cell limit (approximately). Refinement will stop immediately
@ -255,6 +255,8 @@ snapControls
// Settings for the layer addition.
addLayersControls
{
// Are the thickness parameters below relative to the undistorted
// size of the refined cell outside layer (true) or absolute sizes (false).
relativeSizes true;
// Per final patch (so not geometry!) the layer information
@ -277,11 +279,14 @@ addLayersControls
// is the
// thickness of the layer furthest away from the wall.
// Relative to undistorted size of cell outside layer.
// is the thickness of the layer furthest away from the wall.
// See relativeSizes parameter.
finalLayerThickness 0.5;
//- Minimum thickness of cell layer. If for any reason layer
// cannot be above minThickness do not add layer.
// Relative to undistorted size of cell outside layer.
// See relativeSizes parameter.
minThickness 0.25;
//- If points get not extruded do nGrow layers of connected faces that are
@ -323,7 +328,10 @@ addLayersControls
// Create buffer region for new layer terminations
nBufferCellsNoExtrude 0;
// Overall max number of layer addition iterations
// Overall max number of layer addition iterations. The mesher will exit
// if it reaches this number of iterations; possibly with an illegal
// mesh.
nLayerIter 50;
}
@ -349,6 +357,7 @@ meshQualityControls
minFlatness 0.5;
//- Minimum pyramid volume. Is absolute volume of cell pyramid.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
minVol 1e-13;

6
tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/snappyHexMeshDict
View File

@ -210,6 +210,12 @@ castellatedMeshControls
// NOTE: This point should never be on a face, always inside a cell, even
// after refinement.
locationInMesh (0.01 0.01 0.01);
// Whether any faceZones (as specified in the refinementSurfaces)
// are only on the boundary of corresponding cellZones or also allow
// free-standing zone faces. Not used if there are no faceZones.
allowFreeStandingZoneFaces false;
}

18
tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
View File

@ -51,9 +51,9 @@ castellatedMeshControls
// Refinement parameters
// ~~~~~~~~~~~~~~~~~~~~~
// While refining maximum number of cells per processor. This is basically
// the number of cells that fit on a processor. If you choose this too small
// it will do just more refinement iterations to obtain a similar mesh.
// If local number of cells is >= maxLocalCells on any processor
// switches from from refinement followed by balancing
// (current method) to (weighted) balancing before refinement.
maxLocalCells 1000000;
// Overall cell limit (approximately). Refinement will stop immediately
@ -69,6 +69,13 @@ castellatedMeshControls
// (unless the number of cells to refine is 0)
minRefinementCells 10;
// Allow a certain level of imbalance during refining
// (since balancing is quite expensive)
// Expressed as fraction of perfect balance (= overall number of cells /
// nProcs). 0=balance always.
maxLoadUnbalance 0.10;
// Number of buffer layers between different levels.
// 1 means normal 2:1 refinement restriction, larger means slower
// refinement.
@ -180,6 +187,8 @@ snapControls
// Settings for the layer addition.
addLayersControls
{
// Are the thickness parameters below relative to the undistorted
// size of the refined cell outside layer (true) or absolute sizes (false).
relativeSizes true;
// Per final patch (so not geometry!) the layer information
@ -466,6 +475,8 @@ addLayersControls
// is the
// thickness of the layer furthest away from the wall.
// Relative to undistorted size of cell outside layer.
// is the thickness of the layer furthest away from the wall.
// See relativeSizes parameter.
finalLayerThickness 0.3;
//- Minimum thickness of cell layer. If for any reason layer
@ -539,6 +550,7 @@ meshQualityControls
minFlatness 0.5;
//- Minimum pyramid volume. Is absolute volume of cell pyramid.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
minVol 1e-13;

3
tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement
View File

@ -37,7 +37,8 @@ boundaryField
centreOfMass (0.5 0.5 0.5);
momentOfInertia (0.08622222 0.8622222 0.144);
mass 9.6;
rhoInf 1;
rhoInf 1; // for forces calculation
// See sixDoFRigidBodyMotionState
Q (1 0 0 0 1 0 0 0 1);
v (0 0 0);
a (0 0 0);