Merge branch 'master' of github.com:OpenFOAM/OpenFOAM-2.3.x

This commit is contained in:
william 2014-06-05 16:09:01 +01:00
commit 646df63c12
17 changed files with 639 additions and 349 deletions

View file

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

View file

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

View file

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

View file

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -28,9 +28,6 @@ License
#include "polyMesh.H" #include "polyMesh.H"
#include "mapPolyMesh.H" #include "mapPolyMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::faceMapper::calcAddressing() const 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()); insertedFaceLabelsPtr_ = new labelList(mesh_.nFaces());
labelList& insertedFaces = *insertedFaceLabelsPtr_; labelList& insertedFaces = *insertedFaceLabelsPtr_;
@ -413,13 +410,4 @@ const Foam::labelList& Foam::faceMapper::oldPatchSizes() const
} }
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// ************************************************************************* // // ************************************************************************* //

View file

@ -921,7 +921,6 @@ Foam::polyMesh::cellTree() const
} }
// Add boundary patches. Constructor helper
void Foam::polyMesh::addPatches void Foam::polyMesh::addPatches
( (
const List<polyPatch*>& p, const List<polyPatch*>& p,
@ -968,7 +967,6 @@ void Foam::polyMesh::addPatches
} }
// Add mesh zones. Constructor helper
void Foam::polyMesh::addZones void Foam::polyMesh::addZones
( (
const List<pointZone*>& pz, 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 const Foam::pointField& Foam::polyMesh::oldPoints() const
{ {
if (oldPointsPtr_.empty()) if (oldPointsPtr_.empty())
@ -1129,11 +1126,14 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
points_ = newPoints; points_ = newPoints;
bool moveError = false;
if (debug) if (debug)
{ {
// Check mesh motion // Check mesh motion
if (checkMeshMotion(points_, true)) if (checkMeshMotion(points_, true))
{ {
moveError = true;
Info<< "tmp<scalarField> polyMesh::movePoints" Info<< "tmp<scalarField> polyMesh::movePoints"
<< "(const pointField&) : " << "(const pointField&) : "
<< "Moving the mesh with given points will " << "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); const_cast<Time&>(time()).functionObjects().movePoints(*this);
if (debug && moveError)
{
write();
}
return sweptVols; 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 void Foam::polyMesh::removeFiles(const fileName& instanceDir) const
{ {
fileName meshFilesPath = thisDb().time().path()/instanceDir/meshDir(); fileName meshFilesPath = thisDb().time().path()/instanceDir/meshDir();

View file

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -40,56 +40,90 @@ namespace Foam
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // // * * * * * * * * * * * * 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); if (returnReduce(size(), sumOp<label>()))
patchNormal_ /= mag(patchNormal_); {
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_ = centreNormal_ =
vector vector
( (
sign(patchNormal_.x())*(max(mag(patchNormal_.x()), 0.5) - 0.5), sign(n_.x())*(max(mag(n_.x()), 0.5) - 0.5),
sign(patchNormal_.y())*(max(mag(patchNormal_.y()), 0.5) - 0.5), sign(n_.y())*(max(mag(n_.y()), 0.5) - 0.5),
sign(patchNormal_.z())*(max(mag(patchNormal_.z()), 0.5) - 0.5) sign(n_.z())*(max(mag(n_.z()), 0.5) - 0.5)
); );
centreNormal_ /= mag(centreNormal_); centreNormal_ /= mag(centreNormal_);
cosAngle_ = centreNormal_ & n_;
const scalar cnCmptSum = const scalar cnCmptSum =
centreNormal_.x() + centreNormal_.y() + centreNormal_.z(); centreNormal_.x() + centreNormal_.y() + centreNormal_.z();
if (mag(cnCmptSum) < (1 - SMALL)) if (mag(cnCmptSum) < (1 - SMALL))
{ {
FatalErrorIn("wedgePolyPatch::initTransforms()") FatalErrorIn("wedgePolyPatch::calcGeometry(PstreamBuffers&)")
<< "wedge " << name() << "wedge " << name()
<< " centre plane does not align with a coordinate plane by " << " centre plane does not align with a coordinate plane by "
<< 1 - mag(cnCmptSum) << 1 - mag(cnCmptSum)
<< exit(FatalError); << exit(FatalError);
} }
axis_ = centreNormal_ ^ patchNormal_; axis_ = centreNormal_ ^ n_;
scalar magAxis = mag(axis_); scalar magAxis = mag(axis_);
if (magAxis < SMALL) if (magAxis < SMALL)
{ {
FatalErrorIn("wedgePolyPatch::initTransforms()") FatalErrorIn("wedgePolyPatch::calcGeometry(PstreamBuffers&)")
<< "wedge " << name() << "wedge " << name()
<< " plane aligns with a coordinate plane." << nl << " plane aligns with a coordinate plane." << nl
<< " The wedge plane should make a small angle (~2.5deg)" << " The wedge plane should make a small angle (~2.5deg)"
" with the coordinate plane" << nl " with the coordinate plane" << nl
<< " and the the pair of wedge planes should be symmetric" << " and the the pair of wedge planes should be symmetric"
<< " about the coordinate plane." << nl << " about the coordinate plane." << nl
<< " Normal of face " << 0 << " is " << patchNormal_ << " Normal of wedge plane is " << n_
<< " , implied coordinate plane direction is " << centreNormal_ << " , implied coordinate plane direction is " << centreNormal_
<< exit(FatalError); << exit(FatalError);
} }
axis_ /= magAxis; axis_ /= magAxis;
faceT_ = rotationTensor(centreNormal_, patchNormal_); faceT_ = rotationTensor(centreNormal_, n_);
cellT_ = faceT_ & faceT_; cellT_ = faceT_ & faceT_;
} }
} }
@ -107,10 +141,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
const word& patchType const word& patchType
) )
: :
polyPatch(name, size, start, index, bm, patchType) polyPatch(name, size, start, index, bm, patchType),
{ axis_(vector::rootMax),
initTransforms(); centreNormal_(vector::rootMax),
} n_(vector::rootMax),
cosAngle_(0.0),
faceT_(tensor::zero),
cellT_(tensor::zero)
{}
Foam::wedgePolyPatch::wedgePolyPatch Foam::wedgePolyPatch::wedgePolyPatch
@ -122,10 +160,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
const word& patchType const word& patchType
) )
: :
polyPatch(name, dict, index, bm, patchType) polyPatch(name, dict, index, bm, patchType),
{ axis_(vector::rootMax),
initTransforms(); centreNormal_(vector::rootMax),
} n_(vector::rootMax),
cosAngle_(0.0),
faceT_(tensor::zero),
cellT_(tensor::zero)
{}
Foam::wedgePolyPatch::wedgePolyPatch Foam::wedgePolyPatch::wedgePolyPatch
@ -134,10 +176,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
const polyBoundaryMesh& bm const polyBoundaryMesh& bm
) )
: :
polyPatch(pp, bm) polyPatch(pp, bm),
{ axis_(pp.axis_),
initTransforms(); centreNormal_(pp.centreNormal_),
} n_(pp.n_),
cosAngle_(pp.cosAngle_),
faceT_(pp.faceT_),
cellT_(pp.cellT_)
{}
Foam::wedgePolyPatch::wedgePolyPatch Foam::wedgePolyPatch::wedgePolyPatch
@ -149,10 +195,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
const label newStart const label newStart
) )
: :
polyPatch(pp, bm, index, newSize, newStart) polyPatch(pp, bm, index, newSize, newStart),
{ axis_(pp.axis_),
initTransforms(); centreNormal_(pp.centreNormal_),
} n_(pp.n_),
cosAngle_(pp.cosAngle_),
faceT_(pp.faceT_),
cellT_(pp.cellT_)
{}
Foam::wedgePolyPatch::wedgePolyPatch Foam::wedgePolyPatch::wedgePolyPatch
@ -164,10 +214,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
const label newStart const label newStart
) )
: :
polyPatch(pp, bm, index, mapAddressing, newStart) polyPatch(pp, bm, index, mapAddressing, newStart),
{ axis_(pp.axis_),
initTransforms(); centreNormal_(pp.centreNormal_),
} n_(pp.n_),
cosAngle_(pp.cosAngle_),
faceT_(pp.faceT_),
cellT_(pp.cellT_)
{}
// ************************************************************************* // // ************************************************************************* //

View file

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

View file

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

View file

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

View file

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

View file

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -34,6 +34,76 @@ namespace Foam
} }
// * * * * * * * * * * * * * Protected Data Members * * * * * * * * * * * * * //
Foam::IOobject Foam::displacementMotionSolver::points0IO
(
const polyMesh& mesh
) const
{
const word instance =
time().findInstance
(
mesh.meshDir(),
"points0",
IOobject::READ_IF_PRESENT
);
if (instance != time().constant())
{
// points0 written to a time folder
return
IOobject
(
"points0",
instance,
polyMesh::meshSubDir,
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
);
}
else
{
// check that points0 are actually in constant directory
IOobject io
(
"points0",
instance,
polyMesh::meshSubDir,
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
);
if (io.headerOk())
{
return io;
}
else
{
// copy of original mesh points
return
IOobject
(
"points",
instance,
polyMesh::meshSubDir,
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::displacementMotionSolver::displacementMotionSolver Foam::displacementMotionSolver::displacementMotionSolver
@ -49,29 +119,14 @@ Foam::displacementMotionSolver::displacementMotionSolver
IOobject IOobject
( (
"pointDisplacement", "pointDisplacement",
mesh.time().timeName(), time().timeName(),
mesh, mesh,
IOobject::MUST_READ, IOobject::MUST_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
pointMesh::New(mesh) pointMesh::New(mesh)
), ),
points0_ points0_(pointIOField(points0IO(mesh)))
(
pointIOField
(
IOobject
(
"points",
mesh.time().constant(),
polyMesh::meshSubDir,
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
)
)
{ {
if (points0_.size() != mesh.nPoints()) if (points0_.size() != mesh.nPoints())
{ {
@ -81,7 +136,8 @@ Foam::displacementMotionSolver::displacementMotionSolver
"displacementMotionSolver\n" "displacementMotionSolver\n"
"(\n" "(\n"
" const polyMesh&,\n" " const polyMesh&,\n"
" const IOdictionary&\n" " const IOdictionary&,\n"
" const word&\n"
")" ")"
) << "Number of points in mesh " << mesh.nPoints() ) << "Number of points in mesh " << mesh.nPoints()
<< " differs from number of points " << points0_.size() << " differs from number of points " << points0_.size()
@ -90,7 +146,7 @@ Foam::displacementMotionSolver::displacementMotionSolver
IOobject IOobject
( (
"points", "points",
mesh.time().constant(), time().constant(),
polyMesh::meshSubDir, polyMesh::meshSubDir,
mesh, mesh,
IOobject::MUST_READ, IOobject::MUST_READ,
@ -118,7 +174,7 @@ void Foam::displacementMotionSolver::movePoints(const pointField&)
void Foam::displacementMotionSolver::updateMesh(const mapPolyMesh& mpm) void Foam::displacementMotionSolver::updateMesh(const mapPolyMesh& mpm)
{ {
// pointMesh already updates pointFields. // pointMesh already updates pointFields
motionSolver::updateMesh(mpm); motionSolver::updateMesh(mpm);
@ -136,7 +192,7 @@ void Foam::displacementMotionSolver::updateMesh(const mapPolyMesh& mpm)
// Note: boundBox does reduce // Note: boundBox does reduce
const vector span0 = boundBox(points0_).span(); const vector span0 = boundBox(points0_).span();
const vector span = boundBox(points).span(); const vector span = boundBox(points).span();
vector scaleFactors(cmptDivide(span0, span)); vector scaleFactors(cmptDivide(span0, span));
@ -156,11 +212,11 @@ void Foam::displacementMotionSolver::updateMesh(const mapPolyMesh& mpm)
} }
else else
{ {
// New point. Assume motion is scaling. // New point - assume motion is scaling
newPoints0[pointI] = points0_[oldPointI] + cmptMultiply newPoints0[pointI] = points0_[oldPointI] + cmptMultiply
( (
scaleFactors, scaleFactors,
points[pointI]-points[masterPointI] points[pointI] - points[masterPointI]
); );
} }
} }
@ -170,12 +226,21 @@ void Foam::displacementMotionSolver::updateMesh(const mapPolyMesh& mpm)
( (
"displacementMotionSolver::updateMesh" "displacementMotionSolver::updateMesh"
"(const mapPolyMesh&)" "(const mapPolyMesh&)"
) << "Cannot work out coordinates of introduced vertices." ) << "Cannot determine co-ordinates of introduced vertices."
<< " New vertex " << pointI << " at coordinate " << " New vertex " << pointI << " at co-ordinate "
<< points[pointI] << exit(FatalError); << points[pointI] << exit(FatalError);
} }
} }
twoDCorrectPoints(newPoints0);
points0_.transfer(newPoints0); points0_.transfer(newPoints0);
// points0 changed - set to write and check-in to database
points0_.rename("points0");
points0_.writeOpt() = IOobject::AUTO_WRITE;
points0_.instance() = time().timeName();
points0_.checkIn();
} }

View file

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -40,6 +40,7 @@ SourceFiles
#include "motionSolver.H" #include "motionSolver.H"
#include "pointFields.H" #include "pointFields.H"
#include "pointIOField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -63,21 +64,22 @@ protected:
//- Point motion field //- Point motion field
mutable pointVectorField pointDisplacement_; mutable pointVectorField pointDisplacement_;
private:
// Private data
//- Starting points //- Starting points
pointField points0_; pointIOField points0_;
// Protected Member Functions
//- Return IO object for points0
IOobject points0IO(const polyMesh& mesh) const;
private:
// Private Member Functions // Private Member Functions
//- Disallow default bitwise copy construct //- Disallow default bitwise copy construct
displacementMotionSolver displacementMotionSolver(const displacementMotionSolver&);
(
const displacementMotionSolver&
);
//- Disallow default bitwise assignment //- Disallow default bitwise assignment
void operator=(const displacementMotionSolver&); void operator=(const displacementMotionSolver&);

View file

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

View file

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -149,7 +149,7 @@ public:
(-rAB + rA - rB) (-rAB + rA - rB)
*(-rAB - rA + rB) *(-rAB - rA + rB)
*(-rAB + rA + rB) *(-rAB + rA + rB)
*( rAB + rA + rA) *( rAB + rA + rB)
); );
} }

View file

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

View file

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -56,7 +56,7 @@ namespace Foam
class polyMesh; class polyMesh;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class twoDPointCorrector Declaration Class twoDPointCorrector Declaration
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
class twoDPointCorrector class twoDPointCorrector
@ -74,6 +74,15 @@ class twoDPointCorrector
//- Indices of edges normal to plane //- Indices of edges normal to plane
mutable labelList* normalEdgeIndicesPtr_; 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 // Private Member Functions
@ -90,6 +99,9 @@ class twoDPointCorrector
//- Clear addressing //- Clear addressing
void clearAddressing() const; 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 // Static data members

View file

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