diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C index 7fcf53e935..fef49a3fe6 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C @@ -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() + ); + } + + if (!overwrite) { const_cast(mesh.time())++; @@ -444,6 +457,7 @@ int main(int argc, char *argv[]) layerDict, motionDict, layerParams, + preBalance, decomposer, distributor ); diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index 783469ee7d..dab2d1d064 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -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. diff --git a/applications/utilities/mesh/manipulation/setSet/setSet.C b/applications/utilities/mesh/manipulation/setSet/setSet.C index cf24a1c447..bb187dc361 100644 --- a/applications/utilities/mesh/manipulation/setSet/setSet.C +++ b/applications/utilities/mesh/manipulation/setSet/setSet.C @@ -45,6 +45,7 @@ Description #include "cellZoneSet.H" #include "faceZoneSet.H" #include "pointZoneSet.H" +#include "timeSelector.H" #include @@ -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 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,135 +856,171 @@ int main(int argc, char *argv[]) // Print current sets printAllSets(mesh, Info); - - - std::ifstream* fileStreamPtr(NULL); - - if (args.optionFound("batch")) - { - fileName batchFile(args.option("batch")); - - Info<< "Reading commands from file " << batchFile << endl; - - // we cannot handle .gz files - if (!isFile(batchFile, false)) - { - FatalErrorIn(args.executable()) - << "Cannot open file " << batchFile << exit(FatalError); - } - - fileStreamPtr = new std::ifstream(batchFile.c_str()); - } -#if READLINE != 0 - else if (!read_history(historyFile)) + // Read history if interactive +# if READLINE != 0 + if (!batch && !read_history(historyFile)) { Info<< "Successfully read history from " << historyFile << endl; } -#endif +# endif - Info<< "Please type 'help', 'quit' or a set command after prompt." << endl; - bool ok = true; + // Exit status + int status = 0; - FatalError.throwExceptions(); - FatalIOError.throwExceptions(); - do + forAll(timeDirs, timeI) { - string rawLine; + runTime.setTime(timeDirs[timeI], timeI); + Info<< "Time = " << runTime.timeName() << endl; - // Type: cellSet, faceSet, pointSet - word setType; - // Name of destination set. - word setName; - // Action (new, invert etc.) - word actionName; + // Handle geometry/topology changes + meshReadUpdate(mesh); - commandStatus stat = INVALID; - if (fileStreamPtr) + // Main command read & execute loop + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + autoPtr fileStreamPtr(NULL); + + if (batch) { - if (!fileStreamPtr->good()) + fileName batchFile(args.option("batch")); + + Info<< "Reading commands from file " << batchFile << endl; + + // we cannot handle .gz files + if (!isFile(batchFile, false)) { - Info<< "End of batch file" << endl; - break; + FatalErrorIn(args.executable()) + << "Cannot open file " << batchFile << exit(FatalError); } - std::getline(*fileStreamPtr, rawLine); - - if (rawLine.size()) - { - Info<< "Doing:" << rawLine << endl; - } + fileStreamPtr.reset(new IFstream(batchFile)); } - else + + Info<< "Please type 'help', 'quit' or a set command after prompt." + << endl; + + // Whether to quit + bool quit = false; + + FatalError.throwExceptions(); + FatalIOError.throwExceptions(); + + do { -# if READLINE != 0 + string rawLine; + + // Type: cellSet, faceSet, pointSet + word setType; + // Name of destination set. + word setName; + // Action (new, invert etc.) + word actionName; + + commandStatus stat = INVALID; + + if (fileStreamPtr.valid()) { - char* linePtr = readline("readline>"); - - rawLine = string(linePtr); - - if (*linePtr) + if (!fileStreamPtr().good()) { - add_history(linePtr); - write_history(historyFile); + Info<< "End of batch file" << endl; + // No error. + break; } - free(linePtr); // readline uses malloc, not new. - } -# else - { - Info<< "Command>" << flush; - std::getline(std::cin, rawLine); - } -# endif - } + fileStreamPtr().getLine(rawLine); - if (rawLine.empty() || rawLine[0] == '#') - { - continue; - } - - IStringStream is(rawLine + ' '); - - // Type: cellSet, faceSet, pointSet, faceZoneSet - is >> setType; - - stat = parseType(runTime, mesh, setType, is); - - if (stat == VALIDSETCMD || stat == VALIDZONECMD) - { - if (is >> setName) - { - if (is >> actionName) + if (rawLine.size()) { - stat = parseAction(actionName); + Info<< "Doing:" << rawLine << endl; } } - } - ok = true; + else + { +# if READLINE != 0 + { + char* linePtr = readline("readline>"); - if (stat == QUIT) + rawLine = string(linePtr); + + if (*linePtr) + { + add_history(linePtr); + write_history(historyFile); + } + + free(linePtr); // readline uses malloc, not new. + } +# else + { + Info<< "Command>" << flush; + std::getline(std::cin, rawLine); + } +# endif + } + + if (rawLine.empty() || rawLine[0] == '#') + { + continue; + } + + IStringStream is(rawLine + ' '); + + // Type: cellSet, faceSet, pointSet, faceZoneSet + is >> setType; + + stat = parseType(runTime, mesh, setType, is); + + if (stat == VALIDSETCMD || stat == VALIDZONECMD) + { + if (is >> setName) + { + if (is >> actionName) + { + stat = parseAction(actionName); + } + } + } + + if (stat == QUIT) + { + // Make sure to quit + quit = true; + } + else if (stat == VALIDSETCMD || stat == VALIDZONECMD) + { + 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 (!quit); + + if (quit) { break; } - else if (stat == VALIDSETCMD || stat == VALIDZONECMD) - { - ok = doCommand(mesh, setType, setName, actionName, writeVTK, is); - } - - } while (ok); - - - if (fileStreamPtr) - { - delete fileStreamPtr; } Info<< "\nEnd\n" << endl; - return 0; + return status; } diff --git a/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C b/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C index c384eba03d..a534c5e546 100644 --- a/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C +++ b/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C @@ -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; diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C index 001db69e20..36086551cc 100644 --- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C +++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C @@ -400,25 +400,30 @@ void Foam::vtkPV3Foam::updateInfoSets Info<< " Foam::vtkPV3Foam::updateInfoSets" << endl; } - // Add names of sets - IOobjectList objects + // Add names of sets. Search for last time directory with a sets + // subdirectory. Take care not to search beyond the last mesh. + + word facesInstance = dbPtr_().findInstance ( - dbPtr_(), - dbPtr_().findInstance(meshDir_, "faces", IOobject::READ_IF_PRESENT), - meshDir_/"sets" + meshDir_, + "faces", + IOobject::READ_IF_PRESENT ); + 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 " - << dbPtr_().findInstance - ( - meshDir_, - "faces", - IOobject::READ_IF_PRESENT - ) - << endl; + << objects.names() << " from " << setsInstance << endl; } diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict index cb8d07babc..c69f609c39 100644 --- a/applications/utilities/postProcessing/sampling/sample/sampleDict +++ b/applications/utilities/postProcessing/sampling/sample/sampleDict @@ -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; + } ); diff --git a/applications/utilities/surface/surfaceInertia/surfaceInertia.C b/applications/utilities/surface/surfaceInertia/surfaceInertia.C index f3bcb4599b..95edabed0c 100644 --- a/applications/utilities/surface/surfaceInertia/surfaceInertia.C +++ b/applications/utilities/surface/surfaceInertia/surfaceInertia.C @@ -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" ); diff --git a/src/OpenFOAM/db/Time/Time.H b/src/OpenFOAM/db/Time/Time.H index d90cab32ec..9f935091d8 100644 --- a/src/OpenFOAM/db/Time/Time.H +++ b/src/OpenFOAM/db/Time/Time.H @@ -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 diff --git a/src/OpenFOAM/db/Time/findInstance.C b/src/OpenFOAM/db/Time/findInstance.C index a19d3bda63..da7c491fe0 100644 --- a/src/OpenFOAM/db/Time/findInstance.C +++ b/src/OpenFOAM/db/Time/findInstance.C @@ -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); } diff --git a/src/OpenFOAM/dimensionSet/dimensionSets.C b/src/OpenFOAM/dimensionSet/dimensionSets.C index e89f3c944f..693bd27c91 100644 --- a/src/OpenFOAM/dimensionSet/dimensionSets.C +++ b/src/OpenFOAM/dimensionSet/dimensionSets.C @@ -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); // ************************************************************************* // diff --git a/src/OpenFOAM/dimensionSet/dimensionSets.H b/src/OpenFOAM/dimensionSet/dimensionSets.H index 93d9c1aadc..1369cd5e87 100644 --- a/src/OpenFOAM/dimensionSet/dimensionSets.H +++ b/src/OpenFOAM/dimensionSet/dimensionSets.H @@ -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; diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C index 2e90e4d62c..6f711cc603 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C +++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C @@ -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,8 +159,8 @@ Foam::plane::plane(const point& basePoint, const vector& normalVector) else { FatalErrorIn("plane::plane(const point&, const vector&)") - << "plane normal has got zero length" - << abort(FatalError); + << "plane normal has zero length" + << abort(FatalError); } } @@ -217,8 +230,8 @@ Foam::plane::plane(const dictionary& dict) "plane::plane(const dictionary&)", dict ) - << "Invalid plane type: " << planeType - << abort(FatalIOError); + << "Invalid plane type: " << planeType + << abort(FatalIOError); } } @@ -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); } } diff --git a/src/dynamicMesh/setUpdater/setUpdaterTemplates.C b/src/dynamicMesh/setUpdater/setUpdaterTemplates.C index 0239728de2..9c7fcc1a8f 100644 --- a/src/dynamicMesh/setUpdater/setUpdaterTemplates.C +++ b/src/dynamicMesh/setUpdater/setUpdaterTemplates.C @@ -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" ); diff --git a/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchField.C deleted file mode 100644 index 88fee229bf..0000000000 --- a/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchField.C +++ /dev/null @@ -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 -expDirectionMixedFvPatchField::expDirectionMixedFvPatchField -( - const fvPatch& p, - const DimensionedField& iF -) -: - fvPatchField(p, iF), - refValue_(p.size()), - refGrad_(p.size()), - valueFraction_(p.size()) -{} - - -template -expDirectionMixedFvPatchField::expDirectionMixedFvPatchField -( - const expDirectionMixedFvPatchField& ptf, - const fvPatch& p, - const DimensionedField& iF, - const fvPatchFieldMapper& mapper -) -: - fvPatchField(ptf, p, iF, mapper), - refValue_(ptf.refValue_, mapper), - refGrad_(ptf.refGrad_, mapper), - valueFraction_(ptf.valueFraction_, mapper) -{} - - -template -expDirectionMixedFvPatchField::expDirectionMixedFvPatchField -( - const fvPatch& p, - const DimensionedField& iF, - const dictionary& dict -) -: - fvPatchField(p, iF, dict), - refValue_("refValue", dict, p.size()), - refGrad_("refGradient", dict, p.size()), - valueFraction_("valueFraction", dict, p.size()) -{ - evaluate(); -} - - -template -expDirectionMixedFvPatchField::expDirectionMixedFvPatchField -( - const expDirectionMixedFvPatchField& ptf -) -: - fvPatchField(ptf), - refValue_(ptf.refValue_), - refGrad_(ptf.refGrad_), - valueFraction_(ptf.valueFraction_) -{} - - -template -expDirectionMixedFvPatchField::expDirectionMixedFvPatchField -( - const expDirectionMixedFvPatchField& ptf, - const DimensionedField& iF -) -: - fvPatchField(ptf, iF), - refValue_(ptf.refValue_), - refGrad_(ptf.refGrad_), - valueFraction_(ptf.valueFraction_) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template -void expDirectionMixedFvPatchField::autoMap -( - const fvPatchFieldMapper& m -) -{ - fvPatchField::autoMap(m); - refValue_.autoMap(m); - refGrad_.autoMap(m); - valueFraction_.autoMap(m); -} - - -// Reverse-map the given fvPatchField onto this fvPatchField -template -void expDirectionMixedFvPatchField::rmap -( - const fvPatchField& ptf, - const labelList& addr -) -{ - fvPatchField::rmap(ptf, addr); - - const expDirectionMixedFvPatchField& edmptf = - refCast >(ptf); - - refValue_.rmap(edmptf.refValue_, addr); - refGrad_.rmap(edmptf.refGrad_, addr); - valueFraction_.rmap(edmptf.valueFraction_, addr); -} - - -template -tmp > expDirectionMixedFvPatchField::snGrad() const -{ - const vectorField& nHat = patch().faceNormals(); - - Field gradValue = - patchInternalField() + refGrad_/patch().deltaCoeffs(); - - Field mixedValue = - nHat*(nHat & refValue_) - + gradValue - nHat*(nHat & gradValue); - - return - valueFraction_* - (mixedValue - patchInternalField())*patch().deltaCoeffs() - + (1.0 - valueFraction_)*refGrad_; -} - - -template -void expDirectionMixedFvPatchField::evaluate(const Pstream::commsTypes) -{ - if (!updated()) - { - updateCoeffs(); - } - - const vectorField& nHat = patch().faceNormals(); - - Field gradValue = - patchInternalField() + refGrad_/patch().deltaCoeffs(); - - Field mixedValue = - nHat*(nHat & refValue_) - + gradValue - nHat*(nHat & gradValue); - - Field::operator= - ( - valueFraction_*mixedValue + (1.0 - valueFraction_)*gradValue - ); - - fvPatchField::evaluate(); -} - - -template -tmp > expDirectionMixedFvPatchField::valueInternalCoeffs -( - const tmp& -) const -{ - return Type(pTraits::one)*(1.0 - valueFraction_); -} - - -template -tmp > expDirectionMixedFvPatchField::valueBoundaryCoeffs -( - const tmp& -) const -{ - const vectorField& nHat = patch().faceNormals(); - - Field gradValue = - patchInternalField() + refGrad_/patch().deltaCoeffs(); - - Field mixedValue = - nHat*(nHat & refValue_) - + gradValue - nHat*(nHat & gradValue); - - return - valueFraction_*mixedValue - + (1.0 - valueFraction_)*refGrad_/patch().deltaCoeffs(); -} - - -template -tmp > expDirectionMixedFvPatchField:: -gradientInternalCoeffs() const -{ - return -Type(pTraits::one)*valueFraction_*patch().deltaCoeffs(); -} - - -template -tmp > expDirectionMixedFvPatchField:: -gradientBoundaryCoeffs() const -{ - const vectorField& nHat = patch().faceNormals(); - - Field gradValue = - patchInternalField() + refGrad_/patch().deltaCoeffs(); - - Field mixedValue = - nHat*(nHat & refValue_) - + gradValue - nHat*(nHat & gradValue); - - return - valueFraction_*patch().deltaCoeffs()*mixedValue - + (1.0 - valueFraction_)*refGrad_; -} - - -template -void expDirectionMixedFvPatchField::write(Ostream& os) const -{ - fvPatchField::write(os); - refValue_.writeEntry("refValue", os); - refGrad_.writeEntry("refGradient", os); - valueFraction_.writeEntry("valueFraction", os); - writeEntry("value", os); -} - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchField.H deleted file mode 100644 index 009e7a9cb6..0000000000 --- a/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchField.H +++ /dev/null @@ -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 expDirectionMixedFvPatchField -: - public fvPatchField -{ - // Private data - - //- Value field - Field refValue_; - - //- Normal gradient field - Field 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& - ); - - //- Construct from patch, internal field and dictionary - expDirectionMixedFvPatchField - ( - const fvPatch&, - const DimensionedField&, - const dictionary& - ); - - //- Construct by mapping given expDirectionMixedFvPatchField - // onto a new patch - expDirectionMixedFvPatchField - ( - const expDirectionMixedFvPatchField&, - const fvPatch&, - const DimensionedField&, - const fvPatchFieldMapper& - ); - - //- Construct as copy - expDirectionMixedFvPatchField - ( - const expDirectionMixedFvPatchField& - ); - - //- Construct and return a clone - virtual tmp > clone() const - { - return tmp > - ( - new expDirectionMixedFvPatchField(*this) - ); - } - - //- Construct as copy setting internal field reference - expDirectionMixedFvPatchField - ( - const expDirectionMixedFvPatchField&, - const DimensionedField& - ); - - //- Construct and return a clone setting internal field reference - virtual tmp > clone - ( - const DimensionedField& iF - ) const - { - return tmp > - ( - new expDirectionMixedFvPatchField(*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& refValue() - { - return refValue_; - } - - virtual const Field& refValue() const - { - return refValue_; - } - - virtual Field& refGrad() - { - return refGrad_; - } - - virtual const Field& 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&, - const labelList& - ); - - - // Evaluation functions - - //- Return gradient at boundary - virtual tmp > 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 > valueInternalCoeffs - ( - const tmp& - ) const; - - //- Return the matrix source coefficients corresponding to the - // evaluation of the value of this patchField with given weights - virtual tmp > valueBoundaryCoeffs - ( - const tmp& - ) const; - - //- Return the matrix diagonal coefficients corresponding to the - // evaluation of the gradient of this patchField - virtual tmp > gradientInternalCoeffs() const; - - //- Return the matrix source coefficients corresponding to the - // evaluation of the gradient of this patchField - virtual tmp > gradientBoundaryCoeffs() const; - - - //- Write - virtual void write(Ostream&) const; - - - // Member operators - - virtual void operator=(const fvPatchField&) {} - virtual void operator+=(const fvPatchField&) {} - virtual void operator-=(const fvPatchField&) {} - virtual void operator*=(const fvPatchField&) {} - virtual void operator/=(const fvPatchField&) {} - - virtual void operator=(const Field&) {} - virtual void operator+=(const Field&) {} - virtual void operator-=(const Field&) {} - virtual void operator*=(const Field&) {} - virtual void operator/=(const Field&) {} - - 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 - -// ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/AverageIOField.C b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/AverageIOField.C index 2f079100a5..d0c312ff00 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/AverageIOField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/AverageIOField.C @@ -51,7 +51,7 @@ Foam::AverageIOField::AverageIOField : regIOobject(io), Field(size), - average_(0) + average_(pTraits::zero) {} diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C index cc110c1d40..76e589c039 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C @@ -539,6 +539,7 @@ void Foam::autoHexMeshDriver::doMesh() shrinkDict, motionDict, layerParams, + true, // pre-balance decomposerPtr_(), distributorPtr_() ); diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C index acd1673392..d416d11033 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C @@ -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 diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H index fa02d3d0ee..748f26d85e 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H @@ -535,6 +535,7 @@ public: const dictionary& shrinkDict, const dictionary& motionDict, const layerParameters& layerParams, + const bool preBalance, // balance before adding? decompositionMethod& decomposer, fvMeshDistribute& distributor ); diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C index cff924d864..586dd926ce 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C @@ -205,14 +205,36 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine const_cast(mesh.time())++; } - meshRefiner_.refineAndBalance + + if ( - "feature refinement iteration " + name(iter), - decomposer_, - distributor_, - cellsToRefine, - refineParams.maxLoadUnbalance() - ); + returnReduce + ( + (mesh.nCells() >= refineParams.maxLocalCells()), + orOp() + ) + ) + { + meshRefiner_.balanceAndRefine + ( + "feature refinement iteration " + name(iter), + decomposer_, + distributor_, + cellsToRefine, + refineParams.maxLoadUnbalance() + ); + } + else + { + meshRefiner_.refineAndBalance + ( + "feature refinement iteration " + name(iter), + decomposer_, + distributor_, + cellsToRefine, + refineParams.maxLoadUnbalance() + ); + } } } return iter; @@ -306,14 +328,36 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine const_cast(mesh.time())++; } - meshRefiner_.refineAndBalance + + if ( - "surface refinement iteration " + name(iter), - decomposer_, - distributor_, - cellsToRefine, - refineParams.maxLoadUnbalance() - ); + returnReduce + ( + (mesh.nCells() >= refineParams.maxLocalCells()), + orOp() + ) + ) + { + meshRefiner_.balanceAndRefine + ( + "surface refinement iteration " + name(iter), + decomposer_, + distributor_, + cellsToRefine, + refineParams.maxLoadUnbalance() + ); + } + else + { + meshRefiner_.refineAndBalance + ( + "surface refinement iteration " + name(iter), + decomposer_, + distributor_, + cellsToRefine, + 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,14 +535,35 @@ Foam::label Foam::autoRefineDriver::shellRefine const_cast(mesh.time())++; } - meshRefiner_.refineAndBalance + if ( - "shell refinement iteration " + name(iter), - decomposer_, - distributor_, - cellsToRefine, - refineParams.maxLoadUnbalance() - ); + returnReduce + ( + (mesh.nCells() >= refineParams.maxLocalCells()), + orOp() + ) + ) + { + meshRefiner_.balanceAndRefine + ( + "shell refinement iteration " + name(iter), + decomposer_, + distributor_, + cellsToRefine, + refineParams.maxLoadUnbalance() + ); + } + else + { + meshRefiner_.refineAndBalance + ( + "shell refinement iteration " + name(iter), + decomposer_, + distributor_, + cellsToRefine, + refineParams.maxLoadUnbalance() + ); + } } meshRefiner_.userFaceData().clear(); diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H index 6347570cac..7e857c0ba2 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H @@ -656,6 +656,16 @@ public: const scalar maxLoadUnbalance ); + //- Balance before refining some cells + autoPtr balanceAndRefine + ( + const string& msg, + decompositionMethod& decomposer, + fvMeshDistribute& distributor, + const labelList& cellsToRefine, + const scalar maxLoadUnbalance + ); + // Baffle handling diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C index f7adea091d..8d6627882c 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C @@ -598,6 +598,7 @@ Foam::List 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(); diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C index 7f2d0cd617..9fbd5923a0 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C @@ -203,6 +203,7 @@ Foam::Map 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::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; diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C index c77a634d30..fde478a91a 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C @@ -1245,90 +1245,110 @@ Foam::autoPtr Foam::meshRefinement::refine } -//// Do refinement of consistent set of cells followed by truncation and -//// load balancing. -//Foam::autoPtr -//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 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::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 distMap; + + if (Pstream::nProcs() > 1) + { + scalar nIdealCells = + mesh_.globalData().nTotalCells() + / Pstream::nProcs(); + + scalar unbalance = returnReduce + ( + mag(1.0-mesh_.nCells()/nIdealCells), + maxOp() + ); + + 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::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 { diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C index d634603649..7ced800c09 100644 --- a/src/meshTools/meshSearch/meshSearch.C +++ b/src/meshTools/meshSearch/meshSearch.C @@ -813,7 +813,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFace { if (useTreeSearch) { - const indexedOctree& tree = boundaryTree(); + const indexedOctree& tree = boundaryTree(); scalar span = tree.bb().mag(); diff --git a/src/meshTools/momentOfInertia/momentOfInertia.C b/src/meshTools/momentOfInertia/momentOfInertia.C new file mode 100644 index 0000000000..5207850172 --- /dev/null +++ b/src/meshTools/momentOfInertia/momentOfInertia.C @@ -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]; +} + + + +// ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFieldsFwd.H b/src/meshTools/momentOfInertia/momentOfInertia.H similarity index 71% rename from src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFieldsFwd.H rename to src/meshTools/momentOfInertia/momentOfInertia.H index 82899878ff..66a8d357ad 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFieldsFwd.H +++ b/src/meshTools/momentOfInertia/momentOfInertia.H @@ -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 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 +); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/meshTools/momentOfInertia/volInt.ps.gz b/src/meshTools/momentOfInertia/volInt.ps.gz new file mode 100644 index 0000000000..19fb295421 Binary files /dev/null and b/src/meshTools/momentOfInertia/volInt.ps.gz differ diff --git a/src/meshTools/momentOfInertia/volumeIntegration/README b/src/meshTools/momentOfInertia/volumeIntegration/README new file mode 100644 index 0000000000..dadca677b7 --- /dev/null +++ b/src/meshTools/momentOfInertia/volumeIntegration/README @@ -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 + + 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. diff --git a/src/meshTools/momentOfInertia/volumeIntegration/cube b/src/meshTools/momentOfInertia/volumeIntegration/cube new file mode 100644 index 0000000000..07ab26ad16 --- /dev/null +++ b/src/meshTools/momentOfInertia/volumeIntegration/cube @@ -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 + + diff --git a/src/meshTools/momentOfInertia/volumeIntegration/icosa b/src/meshTools/momentOfInertia/volumeIntegration/icosa new file mode 100644 index 0000000000..16339ece89 --- /dev/null +++ b/src/meshTools/momentOfInertia/volumeIntegration/icosa @@ -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 + diff --git a/src/meshTools/momentOfInertia/volumeIntegration/tetra b/src/meshTools/momentOfInertia/volumeIntegration/tetra new file mode 100644 index 0000000000..159e5a2bf8 --- /dev/null +++ b/src/meshTools/momentOfInertia/volumeIntegration/tetra @@ -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 + + + diff --git a/src/meshTools/momentOfInertia/volumeIntegration/volInt.c b/src/meshTools/momentOfInertia/volumeIntegration/volInt.c new file mode 100644 index 0000000000..7601655488 --- /dev/null +++ b/src/meshTools/momentOfInertia/volumeIntegration/volInt.c @@ -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 +#include +#include + +/* + ============================================================================ + 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 \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]); + +} diff --git a/src/meshTools/searchableSurface/searchableSurfaceCollection.C b/src/meshTools/searchableSurface/searchableSurfaceCollection.C index 9279087c19..bef7f79efa 100644 --- a/src/meshTools/searchableSurface/searchableSurfaceCollection.C +++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.C @@ -681,19 +681,22 @@ void Foam::searchableSurfaceCollection::getField List > 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); - const labelList& map = infoMap[surfI]; - forAll(map, i) + if (surfValues.size()) { - values[map[i]] = surfValues[i]; + // Size values only when we have a surface that supports it. + values.setSize(info.size()); + + const labelList& map = infoMap[surfI]; + forAll(map, i) + { + values[map[i]] = surfValues[i]; + } } } } diff --git a/src/meshTools/sets/topoSets/cellSet.C b/src/meshTools/sets/topoSets/cellSet.C index dabf987a12..3fdb4b55dd 100644 --- a/src/meshTools/sets/topoSets/cellSet.C +++ b/src/meshTools/sets/topoSets/cellSet.C @@ -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 diff --git a/src/meshTools/sets/topoSets/topoSet.C b/src/meshTools/sets/topoSets/topoSet.C index 67d9bc4202..0ff02273ea 100644 --- a/src/meshTools/sets/topoSets/topoSet.C +++ b/src/meshTools/sets/topoSets/topoSet.C @@ -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 ) ), diff --git a/src/meshTools/sets/topoSets/topoSet.H b/src/meshTools/sets/topoSets/topoSet.H index 124681305a..bb6d02e64a 100644 --- a/src/meshTools/sets/topoSets/topoSet.H +++ b/src/meshTools/sets/topoSets/topoSet.H @@ -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&); diff --git a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C index 53931689cb..b1b666198c 100644 --- a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C +++ b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C @@ -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); } } diff --git a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.H b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.H index 93bcc7b9a4..7429eefacf 100644 --- a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.H +++ b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.H @@ -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, diff --git a/src/parallel/metisDecomp/metisDecomp.C b/src/parallel/metisDecomp/metisDecomp.C index 9066e92248..fb2fddb868 100644 --- a/src/parallel/metisDecomp/metisDecomp.C +++ b/src/parallel/metisDecomp/metisDecomp.C @@ -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); } } diff --git a/src/parallel/metisDecomp/metisDecomp.H b/src/parallel/metisDecomp/metisDecomp.H index c585686f60..12610f7d82 100644 --- a/src/parallel/metisDecomp/metisDecomp.H +++ b/src/parallel/metisDecomp/metisDecomp.H @@ -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, diff --git a/src/parallel/parMetisDecomp/parMetisDecomp.C b/src/parallel/parMetisDecomp/parMetisDecomp.C index cf7434b5c9..f8e590c957 100644 --- a/src/parallel/parMetisDecomp/parMetisDecomp.C +++ b/src/parallel/parMetisDecomp/parMetisDecomp.C @@ -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); } } diff --git a/src/parallel/parMetisDecomp/parMetisDecomp.H b/src/parallel/parMetisDecomp/parMetisDecomp.H index dda8173557..69c9500a47 100644 --- a/src/parallel/parMetisDecomp/parMetisDecomp.H +++ b/src/parallel/parMetisDecomp/parMetisDecomp.H @@ -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, diff --git a/src/postProcessing/functionObjects/forces/Make/files b/src/postProcessing/functionObjects/forces/Make/files index 511568cab1..0cde9d248e 100644 --- a/src/postProcessing/functionObjects/forces/Make/files +++ b/src/postProcessing/functionObjects/forces/Make/files @@ -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 diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C index 8fb0408e67..dce0389466 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C @@ -160,7 +160,27 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs() t.deltaTValue() ); - Field::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::operator=(motion_.currentPosition(p0_) - p0_); fixedValuePointPatchField::updateCoeffs(); } diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C deleted file mode 100644 index dae1ba9f3d..0000000000 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C +++ /dev/null @@ -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::sixDoFRigidBodyMotion::generatePositions(const pointField& pts) const -{ - return (centreOfMass() + (Q() & (pts - refCentreOfMass_))); -} - - -// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C new file mode 100644 index 0000000000..b89971b17f --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C @@ -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("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, 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, 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; +} + + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H similarity index 60% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H index 16ea1d42b2..ecde4ac6d7 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H @@ -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 restraints_; + + //- Constaints on the motion + PtrList 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& + restraints() const; + + //- Return access to the constraints + inline const PtrList& + 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 generatePositions(const pointField& pts) const; + //- Transform the given reference state pointField by the + // current motion state + inline tmp 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; diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.C new file mode 100644 index 0000000000..9993dfcece --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.C @@ -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; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.H new file mode 100644 index 0000000000..5ba3dcb60d --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.H @@ -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 clone() const + { + return autoPtr + ( + 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 + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.C new file mode 100644 index 0000000000..635267a259 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.C @@ -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; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.H new file mode 100644 index 0000000000..01cc8ff52f --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.H @@ -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 clone() const + { + return autoPtr + ( + 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 + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.C new file mode 100644 index 0000000000..74112ab3b4 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.C @@ -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 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("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; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.H new file mode 100644 index 0000000000..0d46cddea9 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.H @@ -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 clone() const + { + return autoPtr + ( + 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 + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.C new file mode 100644 index 0000000000..6d5ba4962b --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.C @@ -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; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.H new file mode 100644 index 0000000000..45f35dbbb4 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.H @@ -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 clone() const + { + return autoPtr + ( + 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 + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.C new file mode 100644 index 0000000000..5d03ab411f --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.C @@ -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; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.H new file mode 100644 index 0000000000..933ba2db09 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.H @@ -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 clone() const + { + return autoPtr + ( + 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 + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/newSixDoFRigidBodyMotionConstraint.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/newSixDoFRigidBodyMotionConstraint.C new file mode 100644 index 0000000000..890a83b2c0 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/newSixDoFRigidBodyMotionConstraint.C @@ -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::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(cstrIter()(sDoFRBMCDict)); +} + + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C new file mode 100644 index 0000000000..082bbb34af --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C @@ -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("relaxationFactor", 1) + ) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraint::~sixDoFRigidBodyMotionConstraint() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +bool Foam::sixDoFRigidBodyMotionConstraint::read +( + const dictionary& sDoFRBMCDict +) +{ + tolerance_ = (readScalar(sDoFRBMCDict.lookup("tolerance"))); + + relaxationFactor_ = sDoFRBMCDict.lookupOrDefault + ( + "relaxationFactor", + 1 + ); + + sDoFRBMCCoeffs_ = sDoFRBMCDict.subDict(type() + "Coeffs"); + + return true; +} + + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H new file mode 100644 index 0000000000..22e5770202 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H @@ -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 clone() const = 0; + + + // Selectors + + //- Select constructed from the sDoFRBMCDict dictionary and Time + static autoPtr 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 + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H similarity index 66% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H index aeb75bdea4..f112ef42ef 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H @@ -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::sixDoFRigidBodyMotion::restraints() const { - return motionState_.centreOfMass(); + return restraints_; +} + + +inline const Foam::PtrList& +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::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_; +} + // ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C similarity index 61% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C index 005cbbe6aa..55245c3bfe 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C @@ -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 diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C new file mode 100644 index 0000000000..5aa198f2b7 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C @@ -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("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; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.H new file mode 100644 index 0000000000..780598afb2 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.H @@ -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 clone() const + { + return autoPtr + ( + 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 + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.C new file mode 100644 index 0000000000..8067233e1c --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.C @@ -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; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.H new file mode 100644 index 0000000000..c5dc9bc074 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.H @@ -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 clone() const + { + return autoPtr + ( + 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 + +// ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFields.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/newSixDoFRigidBodyMotionRestraint.C similarity index 51% rename from src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFields.H rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/newSixDoFRigidBodyMotionRestraint.C index 8040ff664f..7fee5eb84e 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFields.H +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/newSixDoFRigidBodyMotionRestraint.C @@ -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::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(cstrIter()(sDoFRBMRDict)); +} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif // ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFields.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C similarity index 58% rename from src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFields.C rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C index d0bd6d5fe6..7cfe723c01 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/expDirectionMixed/expDirectionMixedFvPatchFields.C +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C @@ -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 // ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H new file mode 100644 index 0000000000..1bdabb57e3 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H @@ -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 clone() const = 0; + + + // Selectors + + //- Select constructed from the sDoFRBMRDict dictionary and Time + static autoPtr 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 + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.C new file mode 100644 index 0000000000..a961e5591a --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.C @@ -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; cmptname() + << " force " << restraintForce + << " moment " << restraintMoment + << endl; + } +} + + +bool Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::read +( + const dictionary& sDoFRBMRDict +) +{ + sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict); + + refQ_ = sDoFRBMRCoeffs_.lookupOrDefault("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; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.H new file mode 100644 index 0000000000..6ee6636b29 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.H @@ -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 clone() const + { + return autoPtr + ( + 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 + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C new file mode 100644 index 0000000000..159419d44a --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C @@ -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("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(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; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H new file mode 100644 index 0000000000..09cc336e0e --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H @@ -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 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 clone() const + { + return autoPtr + ( + 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 + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C similarity index 100% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.H similarity index 100% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.H rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.H diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H similarity index 100% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C similarity index 100% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C new file mode 100644 index 0000000000..d35eefd671 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C @@ -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& iF +) +: + fixedValuePointPatchField(p, iF), + motion_(), + p0_(p.localPoints()), + rhoInf_(1.0) +{} + + +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField:: +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField +( + const pointPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValuePointPatchField(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& iF, + const pointPatchFieldMapper& mapper +) +: + fixedValuePointPatchField(ptf, p, iF, mapper), + motion_(ptf.motion_), + p0_(ptf.p0_, mapper), + rhoInf_(ptf.rhoInf_) +{} + + +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField:: +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField +( + const uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField& ptf, + const DimensionedField& iF +) +: + fixedValuePointPatchField(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("g")) + { + uniformDimensionedVectorField g = + db().lookupObject("g"); + + gravity = g.value(); + } + + // Do not modify the accelerations, except with gravity, where available + motion_.updateForce(gravity*motion_.mass(), vector::zero, t.deltaTValue()); + + Field::operator=(motion_.currentPosition(p0_) - p0_); + + fixedValuePointPatchField::updateCoeffs(); +} + + +void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::write +( + Ostream& os +) const +{ + pointPatchField::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 + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.H new file mode 100644 index 0000000000..fb4177d104 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.H @@ -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 +{ + // 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& + ); + + //- Construct from patch, internal field and dictionary + uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + ( + const pointPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given patchField onto a new patch + uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + ( + const uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField&, + const pointPatch&, + const DimensionedField&, + const pointPatchFieldMapper& + ); + + //- Construct and return a clone + virtual autoPtr > clone() const + { + return autoPtr > + ( + new uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + ( + *this + ) + ); + } + + //- Construct as copy setting internal field reference + uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + ( + const uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual autoPtr > clone + ( + const DimensionedField& iF + ) const + { + return autoPtr > + ( + 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 + +// ************************************************************************* // diff --git a/src/randomProcesses/fft/kShellIntegration.C b/src/randomProcesses/fft/kShellIntegration.C index 846de03802..9cfdcb8b95 100644 --- a/src/randomProcesses/fft/kShellIntegration.C +++ b/src/randomProcesses/fft/kShellIntegration.C @@ -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 - // ************************************************************************* // - diff --git a/src/sampling/Make/files b/src/sampling/Make/files index e9d78df7a0..2c7920d29e 100644 --- a/src/sampling/Make/files +++ b/src/sampling/Make/files @@ -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 diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C new file mode 100644 index 0000000000..58b200832b --- /dev/null +++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C @@ -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