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

ENH: mesh motion updates

Browse Source
This commit is contained in:
andy 2014-06-03 14:42:39 +01:00 committed by Andrew Heather
parent f46e99668a
commit 709b92d907
14 changed files with 534 additions and 311 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
View File

@ -19,7 +19,7 @@ Foam::label Foam::findOppositeWedge
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
scalar wppCosAngle = wpp.centreNormal()&wpp.patchNormal();
scalar wppCosAngle = wpp.cosAngle();
forAll(patches, patchI)
{
@ -30,13 +30,11 @@ Foam::label Foam::findOppositeWedge
&& isA<wedgePolyPatch>(patches[patchI])
)
{
const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
(
patches[patchI]
);
const wedgePolyPatch& pp =
refCast<const wedgePolyPatch>(patches[patchI]);
// Calculate (cos of) angle to wpp (not pp!) centre normal
scalar ppCosAngle = wpp.centreNormal()&pp.patchNormal();
scalar ppCosAngle = wpp.centreNormal() & pp.n();
if
(
@ -73,12 +71,10 @@ bool Foam::checkWedges
{
if (patches[patchI].size() && isA<wedgePolyPatch>(patches[patchI]))
{
const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
(
patches[patchI]
);
const wedgePolyPatch& pp =
refCast<const wedgePolyPatch>(patches[patchI]);
scalar wedgeAngle = acos(pp.centreNormal()&pp.patchNormal());
scalar wedgeAngle = acos(pp.cosAngle());
if (report)
{
@ -100,10 +96,8 @@ bool Foam::checkWedges
return true;
}
const wedgePolyPatch& opp = refCast<const wedgePolyPatch>
(
patches[oppositePatchI]
);
const wedgePolyPatch& opp =
refCast<const wedgePolyPatch>(patches[oppositePatchI]);
if (mag(opp.axis() & pp.axis()) < (1-1e-3))
@ -140,7 +134,7 @@ bool Foam::checkWedges
forAll(pp.meshPoints(), i)
{
const point& pt = p[pp.meshPoints()[i]];
scalar d = mag((pt-p0) & pp.patchNormal());
scalar d = mag((pt - p0) & pp.n());
if (d > sqrt(SMALL))
{
@ -385,10 +379,8 @@ bool Foam::checkCoupledPoints
{
if (patches[patchI].coupled())
{
const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
(
patches[patchI]
);
const coupledPolyPatch& cpp =
refCast<const coupledPolyPatch>(patches[patchI]);
if (cpp.owner())
{

20
src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -21,9 +21,6 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
Wedge front and back plane patch
\*---------------------------------------------------------------------------*/
#include "wedgePointPatch.H"
@ -46,6 +43,19 @@ namespace Foam
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::wedgePointPatch::wedgePointPatch
(
const polyPatch& patch,
const pointBoundaryMesh& bm
)
:
facePointPatch(patch, bm),
wedgePolyPatch_(refCast<const wedgePolyPatch>(patch))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::wedgePointPatch::applyConstraint
@ -54,7 +64,7 @@ void Foam::wedgePointPatch::applyConstraint
pointConstraint& pc
) const
{
pc.applyConstraint(pointNormals()[pointi]);
pc.applyConstraint(wedgePolyPatch_.n());
}

18
src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -51,6 +51,11 @@ class wedgePointPatch
:
public facePointPatch
{
// Private data
//- Local reference cast into the symmetryPlane patch
const wedgePolyPatch& wedgePolyPatch_;
public:
@ -65,10 +70,7 @@ public:
(
const polyPatch& patch,
const pointBoundaryMesh& bm
)
:
facePointPatch(patch, bm)
{}
);
// Member Functions
@ -85,6 +87,12 @@ public:
const label pointi,
pointConstraint&
) const;
//- Return symmetry plane normal
const vector& n() const
{
return wedgePolyPatch_.n();
}
};

16
src/OpenFOAM/meshes/polyMesh/mapPolyMesh/faceMapper/faceMapper.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -28,9 +28,6 @@ License
#include "polyMesh.H"
#include "mapPolyMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::faceMapper::calcAddressing() const
@ -168,7 +165,7 @@ void Foam::faceMapper::calcAddressing() const
}
// Grab inserted points (for them the size of addressing is still zero)
// Grab inserted faces (for them the size of addressing is still zero)
insertedFaceLabelsPtr_ = new labelList(mesh_.nFaces());
labelList& insertedFaces = *insertedFaceLabelsPtr_;
@ -413,13 +410,4 @@ const Foam::labelList& Foam::faceMapper::oldPatchSizes() const
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// ************************************************************************* //

13
src/OpenFOAM/meshes/polyMesh/polyMesh.C
View File

@ -921,7 +921,6 @@ Foam::polyMesh::cellTree() const
}
// Add boundary patches. Constructor helper
void Foam::polyMesh::addPatches
(
const List<polyPatch*>& p,
@ -968,7 +967,6 @@ void Foam::polyMesh::addPatches
}
// Add mesh zones. Constructor helper
void Foam::polyMesh::addZones
(
const List<pointZone*>& pz,
@ -1084,7 +1082,6 @@ const Foam::labelList& Foam::polyMesh::faceNeighbour() const
}
// Return old mesh motion points
const Foam::pointField& Foam::polyMesh::oldPoints() const
{
if (oldPointsPtr_.empty())
@ -1129,11 +1126,14 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
points_ = newPoints;
bool moveError = false;
if (debug)
{
// Check mesh motion
if (checkMeshMotion(points_, true))
{
moveError = true;
Info<< "tmp<scalarField> polyMesh::movePoints"
<< "(const pointField&) : "
<< "Moving the mesh with given points will "
@ -1176,6 +1176,12 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
const_cast<Time&>(time()).functionObjects().movePoints(*this);
if (debug && moveError)
{
write();
}
return sweptVols;
}
@ -1219,7 +1225,6 @@ Foam::label& Foam::polyMesh::comm()
}
// Remove all files and some subdirs (eg, sets)
void Foam::polyMesh::removeFiles(const fileName& instanceDir) const
{
fileName meshFilesPath = thisDb().time().path()/instanceDir/meshDir();

122
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -40,56 +40,90 @@ namespace Foam
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
void Foam::wedgePolyPatch::initTransforms()
void Foam::wedgePolyPatch::calcGeometry(PstreamBuffers&)
{
if (size() > 0)
if (axis_ != vector::rootMax)
{
const pointField& points = this->points();
return;
}
patchNormal_ = operator[](0).normal(points);
patchNormal_ /= mag(patchNormal_);
if (returnReduce(size(), sumOp<label>()))
{
const vectorField& nf(faceNormals());
n_ = gAverage(nf);
if (debug)
{
Info<< "Patch " << name() << " calculated average normal "
<< n_ << endl;
}
// Check the wedge is planar
forAll(nf, faceI)
{
if (magSqr(n_ - nf[faceI]) > SMALL)
{
// only issue warning instead of error so that the case can
// still be read for post-processing
WarningIn
(
"wedgePolyPatch::calcGeometry(PstreamBuffers&)"
)
<< "Wedge patch '" << name() << "' is not planar." << nl
<< "At local face at "
<< primitivePatch::faceCentres()[faceI]
<< " the normal " << nf[faceI]
<< " differs from the average normal " << n_
<< " by " << magSqr(n_ - nf[faceI]) << nl
<< "Either correct the patch or split it into planar parts"
<< endl;
}
}
centreNormal_ =
vector
(
sign(patchNormal_.x())*(max(mag(patchNormal_.x()), 0.5) - 0.5),
sign(patchNormal_.y())*(max(mag(patchNormal_.y()), 0.5) - 0.5),
sign(patchNormal_.z())*(max(mag(patchNormal_.z()), 0.5) - 0.5)
sign(n_.x())*(max(mag(n_.x()), 0.5) - 0.5),
sign(n_.y())*(max(mag(n_.y()), 0.5) - 0.5),
sign(n_.z())*(max(mag(n_.z()), 0.5) - 0.5)
);
centreNormal_ /= mag(centreNormal_);
cosAngle_ = centreNormal_ & n_;
const scalar cnCmptSum =
centreNormal_.x() + centreNormal_.y() + centreNormal_.z();
if (mag(cnCmptSum) < (1 - SMALL))
{
FatalErrorIn("wedgePolyPatch::initTransforms()")
FatalErrorIn("wedgePolyPatch::calcGeometry(PstreamBuffers&)")
<< "wedge " << name()
<< " centre plane does not align with a coordinate plane by "
<< 1 - mag(cnCmptSum)
<< exit(FatalError);
}
axis_ = centreNormal_ ^ patchNormal_;
axis_ = centreNormal_ ^ n_;
scalar magAxis = mag(axis_);
if (magAxis < SMALL)
{
FatalErrorIn("wedgePolyPatch::initTransforms()")
FatalErrorIn("wedgePolyPatch::calcGeometry(PstreamBuffers&)")
<< "wedge " << name()
<< " plane aligns with a coordinate plane." << nl
<< " The wedge plane should make a small angle (~2.5deg)"
" with the coordinate plane" << nl
<< " and the the pair of wedge planes should be symmetric"
<< " about the coordinate plane." << nl
<< " Normal of face " << 0 << " is " << patchNormal_
<< " Normal of wedge plane is " << n_
<< " , implied coordinate plane direction is " << centreNormal_
<< exit(FatalError);
}
axis_ /= magAxis;
faceT_ = rotationTensor(centreNormal_, patchNormal_);
faceT_ = rotationTensor(centreNormal_, n_);
cellT_ = faceT_ & faceT_;
}
}
@ -107,10 +141,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
const word& patchType
)
:
polyPatch(name, size, start, index, bm, patchType)
{
initTransforms();
}
polyPatch(name, size, start, index, bm, patchType),
axis_(vector::rootMax),
centreNormal_(vector::rootMax),
n_(vector::rootMax),
cosAngle_(0.0),
faceT_(tensor::zero),
cellT_(tensor::zero)
{}
Foam::wedgePolyPatch::wedgePolyPatch
@ -122,10 +160,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
const word& patchType
)
:
polyPatch(name, dict, index, bm, patchType)
{
initTransforms();
}
polyPatch(name, dict, index, bm, patchType),
axis_(vector::rootMax),
centreNormal_(vector::rootMax),
n_(vector::rootMax),
cosAngle_(0.0),
faceT_(tensor::zero),
cellT_(tensor::zero)
{}
Foam::wedgePolyPatch::wedgePolyPatch
@ -134,10 +176,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
const polyBoundaryMesh& bm
)
:
polyPatch(pp, bm)
{
initTransforms();
}
polyPatch(pp, bm),
axis_(pp.axis_),
centreNormal_(pp.centreNormal_),
n_(pp.n_),
cosAngle_(pp.cosAngle_),
faceT_(pp.faceT_),
cellT_(pp.cellT_)
{}
Foam::wedgePolyPatch::wedgePolyPatch
@ -149,10 +195,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
const label newStart
)
:
polyPatch(pp, bm, index, newSize, newStart)
{
initTransforms();
}
polyPatch(pp, bm, index, newSize, newStart),
axis_(pp.axis_),
centreNormal_(pp.centreNormal_),
n_(pp.n_),
cosAngle_(pp.cosAngle_),
faceT_(pp.faceT_),
cellT_(pp.cellT_)
{}
Foam::wedgePolyPatch::wedgePolyPatch
@ -164,10 +214,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
const label newStart
)
:
polyPatch(pp, bm, index, mapAddressing, newStart)
{
initTransforms();
}
polyPatch(pp, bm, index, mapAddressing, newStart),
axis_(pp.axis_),
centreNormal_(pp.centreNormal_),
n_(pp.n_),
cosAngle_(pp.cosAngle_),
faceT_(pp.faceT_),
cellT_(pp.cellT_)
{}
// ************************************************************************* //

26
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -59,7 +59,10 @@ class wedgePolyPatch
vector centreNormal_;
//- Normal to the patch
vector patchNormal_;
vector n_;
//- Cosine of the wedge angle
scalar cosAngle_;
//- Face transformation tensor
tensor faceT_;
@ -67,8 +70,13 @@ class wedgePolyPatch
//- Neighbour-cell transformation tensor
tensor cellT_;
//- Calculate the above tensors
void initTransforms();
protected:
// Protected Member Functions
//- Calculate the patch geometry
virtual void calcGeometry(PstreamBuffers&);
public:
@ -180,9 +188,15 @@ public:
}
//- Return the normal to the patch
const vector& patchNormal() const
const vector& n() const
{
return patchNormal_;
return n_;
}
//- Return the cosine of the wedge angle
scalar cosAngle() const
{
return cosAngle_;
}
//- Return face transformation tensor

118
src/dynamicMesh/layerAdditionRemoval/addCellLayer.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -35,7 +35,6 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
//- Calculate the offset to the next layer
Foam::tmp<Foam::vectorField> Foam::layerAdditionRemoval::extrusionDir() const
{
const polyMesh& mesh = topoChanger().mesh();
@ -79,6 +78,7 @@ Foam::tmp<Foam::vectorField> Foam::layerAdditionRemoval::extrusionDir() const
extrusionDir = minLayerThickness_*masterFaceLayer.pointNormals();
}
return textrusionDir;
}
@ -120,7 +120,7 @@ void Foam::layerAdditionRemoval::addCellLayer
// Get the extrusion direction for the added points
tmp<vectorField> tpointOffsets = extrusionDir();
const vectorField pointOffsets(extrusionDir());
// Add the new points
labelList addedPoints(mp.size());
@ -135,7 +135,7 @@ void Foam::layerAdditionRemoval::addCellLayer
polyAddPoint
(
points[mp[pointI]] // point
+ addDelta_*tpointOffsets()[pointI],
+ addDelta_*pointOffsets[pointI],
mp[pointI], // master point
-1, // zone for point
true // supports a cell
@ -143,14 +143,18 @@ void Foam::layerAdditionRemoval::addCellLayer
);
}
// Pout<< "mp: " << mp << " addedPoints: " << addedPoints << endl;
if (debug > 1)
{
Pout<< "mp: " << mp << " addedPoints: " << addedPoints << endl;
}
// Create the cells
const labelList& mc =
mesh.faceZones()[faceZoneID_.index()].masterCells();
const labelList& sc =
mesh.faceZones()[faceZoneID_.index()].slaveCells();
// Pout<< "mc: " << mc << " sc: " << sc << endl;
const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
@ -231,10 +235,13 @@ void Foam::layerAdditionRemoval::addCellLayer
)
);
// Pout<< "adding face: " << newFace
// << " own: " << mc[faceI]
// << " nei: " << addedCells[faceI]
// << endl;
if (debug > 1)
{
Pout<< "adding face: " << newFace
<< " own: " << mc[faceI]
<< " nei: " << addedCells[faceI]
<< endl;
}
}
// Modify the faces from the master zone for the new neighbour
@ -267,10 +274,14 @@ void Foam::layerAdditionRemoval::addCellLayer
)
);
// Pout<< "Modifying a boundary face. Face: " << curfaceID
// << " flip: " << mfFlip[faceI]
// << endl;
if (debug > 1)
{
Pout<< "Modifying a boundary face. Face: " << curfaceID
<< " flip: " << mfFlip[faceI]
<< endl;
}
}
// If slave cell is owner, the face remains the same (but with
// a new neighbour - the newly created cell). Otherwise, the
// face is flipped.
@ -293,10 +304,13 @@ void Foam::layerAdditionRemoval::addCellLayer
)
);
// Pout<< "modify face, no flip " << curfaceID
// << " own: " << own[curfaceID]
// << " nei: " << addedCells[faceI]
// << endl;
if (debug > 1)
{
Pout<< "modify face, no flip " << curfaceID
<< " own: " << own[curfaceID]
<< " nei: " << addedCells[faceI]
<< endl;
}
}
else
{
@ -317,10 +331,13 @@ void Foam::layerAdditionRemoval::addCellLayer
)
);
// Pout<< "modify face, with flip " << curfaceID
// << " own: " << own[curfaceID]
// << " nei: " << addedCells[faceI]
// << endl;
if (debug > 1)
{
Pout<< "modify face, with flip " << curfaceID
<< " own: " << own[curfaceID]
<< " nei: " << addedCells[faceI]
<< endl;
}
}
}
@ -362,10 +379,13 @@ void Foam::layerAdditionRemoval::addCellLayer
)
);
// Pout<< "Add internal face off edge: " << newFace
// << " own: " << addedCells[edgeFaces[curEdgeID][0]]
// << " nei: " << addedCells[edgeFaces[curEdgeID][1]]
// << endl;
if (debug > 1)
{
Pout<< "Add internal face off edge: " << newFace
<< " own: " << addedCells[edgeFaces[curEdgeID][0]]
<< " nei: " << addedCells[edgeFaces[curEdgeID][1]]
<< endl;
}
}
// Prepare creation of faces from boundary edges.
@ -448,10 +468,13 @@ void Foam::layerAdditionRemoval::addCellLayer
)
);
// Pout<< "add boundary face: " << newFace
// << " into patch " << patchID
// << " own: " << addedCells[edgeFaces[curEdgeID][0]]
// << endl;
if (debug > 1)
{
Pout<< "add boundary face: " << newFace
<< " into patch " << patchID
<< " own: " << addedCells[edgeFaces[curEdgeID][0]]
<< endl;
}
}
// Modify the remaining faces of the master cells to reconnect to the new
@ -562,12 +585,15 @@ void Foam::layerAdditionRemoval::addCellLayer
)
);
// Pout<< "modifying stick-out face. Internal Old face: "
// << oldFace
// << " new face: " << newFace
// << " own: " << own[curFaceID]
// << " nei: " << nei[curFaceID]
// << endl;
if (debug > 1)
{
Pout<< "modifying stick-out face. Internal Old face: "
<< oldFace
<< " new face: " << newFace
<< " own: " << own[curFaceID]
<< " nei: " << nei[curFaceID]
<< endl;
}
}
else
{
@ -588,22 +614,24 @@ void Foam::layerAdditionRemoval::addCellLayer
)
);
// Pout<< "modifying stick-out face. Boundary Old face: "
// << oldFace
// << " new face: " << newFace
// << " own: " << own[curFaceID]
// << " patch: "
// << mesh.boundaryMesh().whichPatch(curFaceID)
// << endl;
if (debug > 1)
{
Pout<< "modifying stick-out face. Boundary Old face: "
<< oldFace
<< " new face: " << newFace
<< " own: " << own[curFaceID]
<< " patch: "
<< mesh.boundaryMesh().whichPatch(curFaceID)
<< endl;
}
}
}
}
if (debug)
{
Pout<< "void layerAdditionRemoval::addCellLayer("
<< "polyTopoChange& ref) const "
<< " for object " << name() << " : "
Pout<< "void layerAdditionRemoval::addCellLayer(polyTopoChange&) const "
<< " for object " << name() << ": "
<< "Finished adding cell layer" << endl;
}
}

90
src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -21,9 +21,6 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
Remove a layer of cells and prepare addressing data
\*---------------------------------------------------------------------------*/
#include "layerAdditionRemoval.H"
@ -87,6 +84,7 @@ bool Foam::layerAdditionRemoval::validCollapse() const
}
}
void Foam::layerAdditionRemoval::removeCellLayer
(
polyTopoChange& ref
@ -115,7 +113,7 @@ void Foam::layerAdditionRemoval::removeCellLayer
// Remove all the cells from the master layer
const labelList& mc =
topoChanger().mesh().faceZones()[faceZoneID_.index()].masterCells();
mesh.faceZones()[faceZoneID_.index()].masterCells();
forAll(mc, faceI)
{
@ -205,7 +203,10 @@ void Foam::layerAdditionRemoval::removeCellLayer
labelList ftm = facesToModify.toc();
//Pout<< "faces to modify: " << ftm << endl;
if (debug > 1)
{
Pout<< "faces to modify: " << ftm << endl;
}
forAll(ftm, faceI)
{
@ -228,9 +229,12 @@ void Foam::layerAdditionRemoval::removeCellLayer
}
}
//Pout<< "face label: " << curFaceID
// << " old face: " << faces[curFaceID]
// << " new face: " << newFace << endl;
if (debug > 1)
{
Pout<< "face label: " << curFaceID
<< " old face: " << faces[curFaceID]
<< " new face: " << newFace << endl;
}
// Get face zone and its flip
label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
@ -238,22 +242,15 @@ void Foam::layerAdditionRemoval::removeCellLayer
if (modifiedFaceZone >= 0)
{
modifiedFaceZoneFlip =
mesh.faceZones()[modifiedFaceZone].flipMap()
[
mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
];
const faceZone& fz = mesh.faceZones()[modifiedFaceZone];
modifiedFaceZoneFlip = fz.flipMap()[fz.whichFace(curFaceID)];
}
label newNei;
label newNeighbour = -1;
if (curFaceID < mesh.nInternalFaces())
{
newNei = nei[curFaceID];
}
else
{
newNei = -1;
newNeighbour = nei[curFaceID];
}
// Modify the face
@ -264,7 +261,7 @@ void Foam::layerAdditionRemoval::removeCellLayer
newFace, // modified face
curFaceID, // label of face being modified
own[curFaceID], // owner
newNei, // neighbour
newNeighbour, // neighbour
false, // face flip
mesh.boundaryMesh().whichPatch(curFaceID),// patch for face
false, // remove from zone
@ -285,20 +282,34 @@ void Foam::layerAdditionRemoval::removeCellLayer
// of the cell to be removed
label masterSideCell = own[mf[faceI]];
if (mesh.isInternalFace(mf[faceI]) && masterSideCell == mc[faceI])
if (masterSideCell == mc[faceI])
{
// Owner cell of the face is being removed.
// Grab the neighbour instead
masterSideCell = nei[mf[faceI]];
if (mesh.isInternalFace(mf[faceI]))
{
// Owner cell of the face is being removed.
// Grab the neighbour instead
masterSideCell = nei[mf[faceI]];
}
else
{
masterSideCell = -1;
}
}
label slaveSideCell = own[ftc[faceI]];
if (mesh.isInternalFace(ftc[faceI]) && slaveSideCell == mc[faceI])
if (slaveSideCell == mc[faceI])
{
// Owner cell of the face is being removed.
// Grab the neighbour instead
slaveSideCell = nei[ftc[faceI]];
if (mesh.isInternalFace(ftc[faceI]))
{
// Owner cell of the face is being removed.
// Grab the neighbour instead
slaveSideCell = nei[ftc[faceI]];
}
else
{
slaveSideCell = -1;
}
}
// Find out if the face needs to be flipped
@ -374,15 +385,18 @@ void Foam::layerAdditionRemoval::removeCellLayer
zoneFlip = !zoneFlip;
}
// Pout<< "Modifying face " << mf[faceI]
// << " newFace: " << newFace << nl
// << " newOwner: " << newOwner
// << " newNeighbour: " << newNeighbour
// << " flipFace: " << flipFace
// << " newPatchID: " << newPatchID
// << " newZoneID: " << newZoneID << nl
// << " oldOwn: " << own[mf[faceI]]
// << " oldNei: " << nei[mf[faceI]] << endl;
if (debug > 1)
{
Pout<< "Modifying face " << mf[faceI]
<< " newFace: " << newFace << nl
<< " newOwner: " << newOwner
<< " newNeighbour: " << newNeighbour
<< " flipFace: " << flipFace
<< " newPatchID: " << newPatchID
<< " newZoneID: " << newZoneID << nl
<< " oldOwn: " << own[mf[faceI]]
<< " oldNei: " << nei[mf[faceI]] << endl;
}
ref.setAction
(

48
src/dynamicMesh/layerAdditionRemoval/setLayerPairing.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -88,10 +88,14 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
facesPairingPtr_ = new labelList(mf.size(), -1);
labelList& ftc = *facesPairingPtr_;
// Pout<< "meshPoints: " << meshPoints << nl
// << "localPoints: "
// << mesh.faceZones()[faceZoneID_.index()]().localPoints()
// << endl;
if (debug > 1)
{
Pout<< "meshPoints: " << meshPoints << nl
<< "localPoints: "
<< mesh.faceZones()[faceZoneID_.index()]().localPoints()
<< endl;
}
// For all faces, create the mapping
label nPointErrors = 0;
@ -119,12 +123,15 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
continue;
}
// Pout<< "curMasterFace: " << faces[mf[faceI]] << nl
// << "cell shape: " << mesh.cellShapes()[mc[faceI]] << nl
// << "curLocalFace: " << curLocalFace << nl
// << "lidFace: " << lidFace
// << " master index: " << lidFace.masterIndex()
// << " oppositeIndex: " << lidFace.oppositeIndex() << endl;
if (debug > 1)
{
Pout<< "curMasterFace: " << faces[mf[faceI]] << nl
<< "cell shape: " << mesh.cellShapes()[mc[faceI]] << nl
<< "curLocalFace: " << curLocalFace << nl
<< "lidFace: " << lidFace
<< " master index: " << lidFace.masterIndex()
<< " oppositeIndex: " << lidFace.oppositeIndex() << endl;
}
// Grab the opposite face for face collapse addressing
ftc[faceI] = lidFace.oppositeIndex();
@ -145,17 +152,20 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
if (ptc[clp] != lidFace[pointI])
{
nPointErrors++;
// Pout<< "Topological error in cell layer pairing. "
// << "This mesh is either topologically incorrect "
// << "or the master afce layer is not defined "
// << "consistently. Please check the "
// << "face zone flip map." << nl
// << "First index: " << ptc[clp]
// << " new index: " << lidFace[pointI] << endl;
if (debug > 1)
{
Pout<< "Topological error in cell layer pairing. "
<< "This mesh is either topologically incorrect "
<< "or the master face layer is not defined "
<< "consistently. Please check the "
<< "face zone flip map." << nl
<< "First index: " << ptc[clp]
<< " new index: " << lidFace[pointI] << endl;
}
}
}
}
// Pout<< "ptc: " << ptc << endl;
}
reduce(nPointErrors, sumOp<label>());

178
src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -30,6 +30,8 @@ License
#include "PointEdgeWave.H"
#include "syncTools.H"
#include "interpolationTable.H"
#include "mapPolyMesh.H"
#include "pointConstraints.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -111,11 +113,13 @@ void Foam::displacementLayeredMotionMotionSolver::calcZoneMask
0
);
Info<< "On cellZone " << cz.name()
<< " marked " << returnReduce(nPoints, sumOp<label>())
<< " points and " << returnReduce(nEdges, sumOp<label>())
<< " edges." << endl;
if (debug)
{
Info<< "On cellZone " << cz.name()
<< " marked " << returnReduce(nPoints, sumOp<label>())
<< " points and " << returnReduce(nEdges, sumOp<label>())
<< " edges." << endl;
}
}
}
@ -239,14 +243,21 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
{
FatalIOErrorIn
(
"displacementLayeredMotionMotionSolver::faceZoneEvaluate(..)",
"displacementLayeredMotionMotionSolver::faceZoneEvaluate"
"("
"const faceZone&, "
"const labelList&, "
"const dictionary&, "
"const PtrList<pointVectorField>&, "
"const label"
") const",
*this
) << "slip can only be used on second faceZonePatch of pair."
) << "slip can only be used on second faceZone patch of pair. "
<< "FaceZone:" << fz.name()
<< exit(FatalIOError);
}
// Use field set by previous bc
fld = vectorField(patchDisp[patchI-1], meshPoints);
fld = vectorField(patchDisp[patchI - 1], meshPoints);
}
else if (type == "follow")
{
@ -269,7 +280,14 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
{
FatalIOErrorIn
(
"displacementLayeredMotionMotionSolver::faceZoneEvaluate(..)",
"displacementLayeredMotionMotionSolver::faceZoneEvaluate"
"("
"const faceZone&, "
"const labelList&, "
"const dictionary&, "
"const PtrList<pointVectorField>&, "
"const label"
") const",
*this
) << "Unknown faceZonePatch type " << type << " for faceZone "
<< fz.name() << exit(FatalIOError);
@ -295,9 +313,9 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
FatalIOErrorIn
(
"displacementLayeredMotionMotionSolver::"
"correctBoundaryConditions(..)",
"cellZoneSolve(const label, const dictionary&)",
*this
) << "Can only handle 2 faceZones (= patches) per cellZone. "
) << "Two faceZones (patches) must be specifed per cellZone. "
<< " cellZone:" << cellZoneI
<< " patches:" << patchesDict.toc()
<< exit(FatalIOError);
@ -317,7 +335,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
FatalIOErrorIn
(
"displacementLayeredMotionMotionSolver::"
"correctBoundaryConditions(..)",
"cellZoneSolve(const label, const dictionary&)",
*this
) << "Cannot find faceZone " << faceZoneName
<< endl << "Valid zones are " << mesh().faceZones().names()
@ -386,14 +404,17 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
patchI
);
Info<< "For cellZone:" << cellZoneI
<< " for faceZone:" << fz.name() << " nPoints:" << tseed().size()
<< " have patchField:"
<< " max:" << gMax(tseed())
<< " min:" << gMin(tseed())
<< " avg:" << gAverage(tseed())
<< endl;
if (debug)
{
Info<< "For cellZone:" << cellZoneI
<< " for faceZone:" << fz.name()
<< " nPoints:" << tseed().size()
<< " have patchField:"
<< " max:" << gMax(tseed())
<< " min:" << gMin(tseed())
<< " avg:" << gAverage(tseed())
<< endl;
}
// Set distance and transported value
walkStructured
@ -417,16 +438,15 @@ Info<< "For cellZone:" << cellZoneI
// Solve
// ~~~~~
// solving the interior is just interpolating
if (debug)
{
// Get normalised distance
// Normalised distance
pointScalarField distance
(
IOobject
(
"distance",
mesh().cellZones()[cellZoneI].name() + ":distance",
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
@ -434,40 +454,66 @@ Info<< "For cellZone:" << cellZoneI
false
),
pointMesh::New(mesh()),
dimensionedScalar("distance", dimLength, 0.0)
dimensionedScalar("zero", dimLength, 0.0)
);
forAll(distance, pointI)
{
scalar d1 = patchDist[0][pointI];
scalar d2 = patchDist[1][pointI];
if (d1 + d2 > SMALL)
{
scalar s = d1/(d1 + d2);
distance[pointI] = s;
}
}
Info<< "Writing " << pointScalarField::typeName
<< distance.name() << " to "
<< mesh().time().timeName() << endl;
distance.write();
}
const word interpolationScheme = zoneDict.lookup("interpolationScheme");
if (interpolationScheme == "oneSided")
{
forAll(pointDisplacement_, pointI)
{
if (isZonePoint[pointI])
{
pointDisplacement_[pointI] = patchDisp[0][pointI];
}
}
}
else if (interpolationScheme == "linear")
{
forAll(pointDisplacement_, pointI)
{
if (isZonePoint[pointI])
{
scalar d1 = patchDist[0][pointI];
scalar d2 = patchDist[1][pointI];
if (d1+d2 > SMALL)
{
scalar s = d1/(d1+d2);
distance[pointI] = s;
}
scalar s = d1/(d1 + d2 + VSMALL);
const vector& pd1 = patchDisp[0][pointI];
const vector& pd2 = patchDisp[1][pointI];
pointDisplacement_[pointI] = (1 - s)*pd1 + s*pd2;
}
}
Info<< "Writing distance pointScalarField to "
<< mesh().time().timeName() << endl;
distance.write();
}
// Average
forAll(pointDisplacement_, pointI)
else
{
if (isZonePoint[pointI])
{
scalar d1 = patchDist[0][pointI];
scalar d2 = patchDist[1][pointI];
scalar s = d1/(d1+d2+VSMALL);
pointDisplacement_[pointI] =
(1-s)*patchDisp[0][pointI]
+ s*patchDisp[1][pointI];
}
FatalErrorIn
(
"displacementLayeredMotionMotionSolver::"
"cellZoneSolve(const label, const dictionary&)"
)
<< "Invalid interpolationScheme: " << interpolationScheme
<< ". Valid schemes are 'oneSided' and 'linear'"
<< exit(FatalError);
}
}
@ -502,23 +548,19 @@ Foam::displacementLayeredMotionMotionSolver::curPoints() const
points0() + pointDisplacement_.internalField()
);
twoDCorrectPoints(tcurPoints());
return tcurPoints;
}
void Foam::displacementLayeredMotionMotionSolver::solve()
{
// The points have moved so before interpolation update
// the motionSolver accordingly
// The points have moved so before interpolation update the motionSolver
movePoints(mesh().points());
// Apply boundary conditions
pointDisplacement_.boundaryField().updateCoeffs();
// Apply all regions (=cellZones)
// Solve motion on all regions (=cellZones)
const dictionary& regionDicts = coeffDict().subDict("regions");
forAllConstIter(dictionary, regionDicts, regionIter)
{
@ -527,14 +569,13 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
label zoneI = mesh().cellZones().findZoneID(cellZoneName);
Info<< "solve : zone:" << cellZoneName << " index:" << zoneI
<< endl;
Info<< "solving for zone: " << cellZoneName << endl;
if (zoneI == -1)
{
FatalIOErrorIn
(
"displacementLayeredMotionMotionSolver::solve(..)",
"displacementLayeredMotionMotionSolver::solve()",
*this
) << "Cannot find cellZone " << cellZoneName
<< endl << "Valid zones are " << mesh().cellZones().names()
@ -545,7 +586,9 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
}
// Update pointDisplacement for solved values
pointDisplacement_.correctBoundaryConditions();
const pointConstraints& pcs =
pointConstraints::New(pointDisplacement_.mesh());
pcs.constrainDisplacement(pointDisplacement_, false);
}
@ -555,6 +598,27 @@ void Foam::displacementLayeredMotionMotionSolver::updateMesh
)
{
displacementMotionSolver::updateMesh(mpm);
const vectorField displacement(this->newPoints() - points0_);
forAll(points0_, pointI)
{
label oldPointI = mpm.pointMap()[pointI];
if (oldPointI >= 0)
{
label masterPointI = mpm.reversePointMap()[oldPointI];
if ((masterPointI != pointI))
{
// newly inserted point in this cellZone
// need to set point0 so that it represents the position that
// it would have had if it had existed for all time
points0_[pointI] -= displacement[pointI];
}
}
}
}

134
src/meshTools/twoDPointCorrector/twoDPointCorrector.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -28,26 +28,20 @@ License
#include "wedgePolyPatch.H"
#include "emptyPolyPatch.H"
#include "SubField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
#include "meshTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const scalar twoDPointCorrector::edgeOrthogonalityTol = 1.0 - 1e-4;
const Foam::scalar Foam::twoDPointCorrector::edgeOrthogonalityTol = 1.0 - 1e-4;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void twoDPointCorrector::calcAddressing() const
void Foam::twoDPointCorrector::calcAddressing() const
{
// Find geometry normal
planeNormalPtr_ = new vector(0, 0, 0);
vector& pn = *planeNormalPtr_;
bool isWedge = false;
// Algorithm:
// Attempt to find wedge patch and work out the normal from it.
// If not found, find an empty patch with faces in it and use the
@ -61,9 +55,15 @@ void twoDPointCorrector::calcAddressing() const
{
if (isA<wedgePolyPatch>(patches[patchI]))
{
isWedge = true;
isWedge_ = true;
pn = refCast<const wedgePolyPatch>(patches[patchI]).centreNormal();
const wedgePolyPatch& wp =
refCast<const wedgePolyPatch>(patches[patchI]);
pn = wp.centreNormal();
wedgeAxis_ = wp.axis();
wedgeAngle_ = mag(acos(wp.cosAngle()));
if (polyMesh::debug)
{
@ -75,7 +75,7 @@ void twoDPointCorrector::calcAddressing() const
}
// Try to find an empty patch with faces
if (!isWedge)
if (!isWedge_)
{
forAll(patches, patchI)
{
@ -96,11 +96,8 @@ void twoDPointCorrector::calcAddressing() const
if (mag(pn) < VSMALL)
{
FatalErrorIn
(
"twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh, "
"const vector& n)"
) << "Cannot determine normal vector from patches."
FatalErrorIn("twoDPointCorrector::calcAddressing()")
<< "Cannot determine normal vector from patches."
<< abort(FatalError);
}
else
@ -125,9 +122,9 @@ void twoDPointCorrector::calcAddressing() const
forAll(meshEdges, edgeI)
{
vector edgeVector =
meshEdges[edgeI].vec(meshPoints)/
(meshEdges[edgeI].mag(meshPoints) + VSMALL);
const edge& e = meshEdges[edgeI];
vector edgeVector = e.vec(meshPoints)/(e.mag(meshPoints) + VSMALL);
if (mag(edgeVector & pn) > edgeOrthogonalityTol)
{
@ -141,15 +138,12 @@ void twoDPointCorrector::calcAddressing() const
// Construction check: number of points in a read 2-D or wedge geometry
// should be odd and the number of edges normal to the plane should be
// exactly half the number of points
if (!isWedge)
if (!isWedge_)
{
if (meshPoints.size() % 2 != 0)
{
WarningIn
(
"twoDPointCorrector::twoDPointCorrector("
"const polyMesh& mesh, const vector& n)"
) << "the number of vertices in the geometry "
WarningIn("twoDPointCorrector::calcAddressing()")
<< "the number of vertices in the geometry "
<< "is odd - this should not be the case for a 2-D case. "
<< "Please check the geometry."
<< endl;
@ -157,11 +151,8 @@ void twoDPointCorrector::calcAddressing() const
if (2*nNormalEdges != meshPoints.size())
{
WarningIn
(
"twoDPointCorrector::twoDPointCorrector("
"const polyMesh& mesh, const vector& n)"
) << "The number of points in the mesh is "
WarningIn("twoDPointCorrector::calcAddressing()")
<< "The number of points in the mesh is "
<< "not equal to twice the number of edges normal to the plane "
<< "- this may be OK only for wedge geometries.\n"
<< " Please check the geometry or adjust "
@ -174,21 +165,38 @@ void twoDPointCorrector::calcAddressing() const
}
void twoDPointCorrector::clearAddressing() const
void Foam::twoDPointCorrector::clearAddressing() const
{
deleteDemandDrivenData(planeNormalPtr_);
deleteDemandDrivenData(normalEdgeIndicesPtr_);
}
void Foam::twoDPointCorrector::snapToWedge
(
const vector& n,
const point& A,
point& p
) const
{
scalar ADash = mag(A - wedgeAxis_*(wedgeAxis_ & A));
vector pDash = ADash*tan(wedgeAngle_)*planeNormal();
p = A + sign(n & p)*pDash;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh)
Foam::twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh)
:
MeshObject<polyMesh, Foam::UpdateableMeshObject, twoDPointCorrector>(mesh),
required_(mesh_.nGeometricD() == 2),
planeNormalPtr_(NULL),
normalEdgeIndicesPtr_(NULL)
normalEdgeIndicesPtr_(NULL),
isWedge_(false),
wedgeAxis_(vector::zero),
wedgeAngle_(0.0)
{}
@ -203,7 +211,7 @@ Foam::twoDPointCorrector::~twoDPointCorrector()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
direction twoDPointCorrector::normalDir() const
Foam::direction Foam::twoDPointCorrector::normalDir() const
{
const vector& pn = planeNormal();
@ -231,8 +239,7 @@ direction twoDPointCorrector::normalDir() const
}
// Return plane normal
const vector& twoDPointCorrector::planeNormal() const
const Foam::vector& Foam::twoDPointCorrector::planeNormal() const
{
if (!planeNormalPtr_)
{
@ -243,8 +250,7 @@ const vector& twoDPointCorrector::planeNormal() const
}
// Return indices of normal edges.
const labelList& twoDPointCorrector::normalEdgeIndices() const
const Foam::labelList& Foam::twoDPointCorrector::normalEdgeIndices() const
{
if (!normalEdgeIndicesPtr_)
{
@ -255,7 +261,7 @@ const labelList& twoDPointCorrector::normalEdgeIndices() const
}
void twoDPointCorrector::correctPoints(pointField& p) const
void Foam::twoDPointCorrector::correctPoints(pointField& p) const
{
if (!required_) return;
@ -277,16 +283,25 @@ void twoDPointCorrector::correctPoints(pointField& p) const
point& pEnd = p[meshEdges[neIndices[edgeI]].end()];
// calculate average point position
const point A = 0.5*(pStart + pEnd);
point A = 0.5*(pStart + pEnd);
meshTools::constrainToMeshCentre(mesh_, A);
// correct point locations
pStart = A + pn*(pn & (pStart - A));
pEnd = A + pn*(pn & (pEnd - A));
if (isWedge_)
{
snapToWedge(pn, A, pStart);
snapToWedge(pn, A, pEnd);
}
else
{
// correct point locations
pStart = A + pn*(pn & (pStart - A));
pEnd = A + pn*(pn & (pEnd - A));
}
}
}
void twoDPointCorrector::correctDisplacement
void Foam::twoDPointCorrector::correctDisplacement
(
const pointField& p,
vectorField& disp
@ -316,29 +331,36 @@ void twoDPointCorrector::correctDisplacement
point pEnd = p[endPointI] + disp[endPointI];
// calculate average point position
const point A = 0.5*(pStart + pEnd);
point A = 0.5*(pStart + pEnd);
meshTools::constrainToMeshCentre(mesh_, A);
// correct point locations
disp[startPointI] = (A + pn*(pn & (pStart - A))) - p[startPointI];
disp[endPointI] = (A + pn*(pn & (pEnd - A))) - p[endPointI];
if (isWedge_)
{
snapToWedge(pn, A, pStart);
snapToWedge(pn, A, pEnd);
disp[startPointI] = pStart - p[startPointI];
disp[endPointI] = pEnd - p[endPointI];
}
else
{
// correct point locations
disp[startPointI] = (A + pn*(pn & (pStart - A))) - p[startPointI];
disp[endPointI] = (A + pn*(pn & (pEnd - A))) - p[endPointI];
}
}
}
void twoDPointCorrector::updateMesh(const mapPolyMesh&)
void Foam::twoDPointCorrector::updateMesh(const mapPolyMesh&)
{
clearAddressing();
}
bool twoDPointCorrector::movePoints()
bool Foam::twoDPointCorrector::movePoints()
{
return true;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

16
src/meshTools/twoDPointCorrector/twoDPointCorrector.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -56,7 +56,7 @@ namespace Foam
class polyMesh;
/*---------------------------------------------------------------------------*\
Class twoDPointCorrector Declaration
Class twoDPointCorrector Declaration
\*---------------------------------------------------------------------------*/
class twoDPointCorrector
@ -74,6 +74,15 @@ class twoDPointCorrector
//- Indices of edges normal to plane
mutable labelList* normalEdgeIndicesPtr_;
//- Flag to indicate a wedge geometry
mutable bool isWedge_;
//- Wedge axis (if wedge geometry)
mutable vector wedgeAxis_;
//- Wedge angle (if wedge geometry)
mutable scalar wedgeAngle_;
// Private Member Functions
@ -90,6 +99,9 @@ class twoDPointCorrector
//- Clear addressing
void clearAddressing() const;
//- Snap a point to the wedge patch(es)
void snapToWedge(const vector& n, const point& A, point& p) const;
// Static data members

14
src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
View File

@ -30,8 +30,8 @@ License
bool Foam::sixDoFRigidBodyMotion::read(const dictionary& dict)
{
dict.lookup("momentOfInertia") >> momentOfInertia_;
dict.lookup("mass") >> mass_;
dict.lookup("momentOfInertia") >> momentOfInertia_;
aRelax_ = dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0);
aDamp_ = dict.lookupOrDefault<scalar>("accelerationDamping", 1.0);
report_ = dict.lookupOrDefault<Switch>("report", false);
@ -50,16 +50,18 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
{
motionState_.write(os);
os.writeKeyword("initialCentreOfRotation")
<< initialCentreOfRotation_ << token::END_STATEMENT << nl;
os.writeKeyword("initialOrientation")
os.writeKeyword("centreOfMass")
<< initialCentreOfMass_ << token::END_STATEMENT << nl;
os.writeKeyword("orientation")
<< initialQ_ << token::END_STATEMENT << nl;
os.writeKeyword("momentOfInertia")
<< momentOfInertia_ << token::END_STATEMENT << nl;
os.writeKeyword("mass")
<< mass_ << token::END_STATEMENT << nl;
os.writeKeyword("momentOfInertia")
<< momentOfInertia_ << token::END_STATEMENT << nl;
os.writeKeyword("accelerationRelaxation")
<< aRelax_ << token::END_STATEMENT << nl;
os.writeKeyword("accelerationDamping")
<< aDamp_ << token::END_STATEMENT << nl;
os.writeKeyword("report")
<< report_ << token::END_STATEMENT << nl;