diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index 32e260f5fa99435c2c0dabb44fe173fa36dc3d96..8677f08c65919bdaebb8fb16de3cc1ab4da3c20b 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -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())
{
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.C b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.C
index c335911c55a42c247970d1e071c33bef2abb2d6d..56dba9971a2af0729ca22671a86cb699fee818c0 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.C
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.C
@@ -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());
}
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.H
index a62acfaa90d29bc9258772b4e5042ccae32563a3..0a437c321882b2f75a4a3646edd987f4bb582829 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.H
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.H
@@ -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();
+ }
};
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/faceMapper/faceMapper.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/faceMapper/faceMapper.C
index e95fbe95f9149a8d439c82691291a497e1ecd68d..803e523d41db27ee9506eb61dae23bcb99a08426 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/faceMapper/faceMapper.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/faceMapper/faceMapper.C
@@ -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 * * * * * * * * * * * * * //
-
-
// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index cd70eb39991b7dededca438995b21f74bd1db6bb..668a595321efd913be23f707f7dfa193426003fe 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -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();
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C
index b7c9783cb984dad0c2cf0691c9afd93fac55b36b..e12d22e5a5c43bf388bd44d5737af2f534f3e01e 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C
@@ -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;
+ }
+
+ if (returnReduce(size(), sumOp<label>()))
+ {
+ const vectorField& nf(faceNormals());
+ n_ = gAverage(nf);
+
+ if (debug)
+ {
+ Info<< "Patch " << name() << " calculated average normal "
+ << n_ << endl;
+ }
- patchNormal_ = operator[](0).normal(points);
- patchNormal_ /= mag(patchNormal_);
+
+ // 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_)
+{}
// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H
index 2caa492d9f143c60b6a656c75fdf859bede26b2f..9e434cd636db58ef7fe32d429ccada434b83449f 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H
@@ -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 n_;
+ }
+
+ //- Return the cosine of the wedge angle
+ scalar cosAngle() const
{
- return patchNormal_;
+ return cosAngle_;
}
//- Return face transformation tensor
diff --git a/src/dynamicMesh/layerAdditionRemoval/addCellLayer.C b/src/dynamicMesh/layerAdditionRemoval/addCellLayer.C
index c05731781fd8d80194118c61d9f6cf5c05237448..0fffb4aa359377305fc11fcaf114568e7c888931 100644
--- a/src/dynamicMesh/layerAdditionRemoval/addCellLayer.C
+++ b/src/dynamicMesh/layerAdditionRemoval/addCellLayer.C
@@ -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;
}
}
diff --git a/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C b/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
index 955e07a1dab7ad62d164462b5874c4411c60c1e0..5569cc876ef92778caffc921e3b2400484376c7c 100644
--- a/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
+++ b/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
@@ -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
(
diff --git a/src/dynamicMesh/layerAdditionRemoval/setLayerPairing.C b/src/dynamicMesh/layerAdditionRemoval/setLayerPairing.C
index d41241617757dac034629e5819506247dc7f5913..0e2996e45ecda8fd82041d87f4237f4870dea5d1 100644
--- a/src/dynamicMesh/layerAdditionRemoval/setLayerPairing.C
+++ b/src/dynamicMesh/layerAdditionRemoval/setLayerPairing.C
@@ -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>());
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
index 3ea1d3ed39e0fd91908d8537eb659311c12e9642..137283c1ba89785cd517f785e944a2b7eea1ccc4 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
@@ -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,41 +454,67 @@ Info<< "For cellZone:" << cellZoneI
false
),
pointMesh::New(mesh()),
- dimensionedScalar("distance", dimLength, 0.0)
+ dimensionedScalar("zero", dimLength, 0.0)
);
+
forAll(distance, pointI)
{
- if (isZonePoint[pointI])
+ scalar d1 = patchDist[0][pointI];
+ scalar d2 = patchDist[1][pointI];
+ if (d1 + d2 > SMALL)
{
- 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);
+ distance[pointI] = s;
}
}
- Info<< "Writing distance pointScalarField to "
+
+ Info<< "Writing " << pointScalarField::typeName
+ << distance.name() << " to "
<< mesh().time().timeName() << endl;
distance.write();
}
- // Average
- forAll(pointDisplacement_, pointI)
+
+ const word interpolationScheme = zoneDict.lookup("interpolationScheme");
+
+ if (interpolationScheme == "oneSided")
{
- if (isZonePoint[pointI])
+ forAll(pointDisplacement_, pointI)
{
- scalar d1 = patchDist[0][pointI];
- scalar d2 = patchDist[1][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];
+ scalar s = d1/(d1 + d2 + VSMALL);
- scalar s = d1/(d1+d2+VSMALL);
+ const vector& pd1 = patchDisp[0][pointI];
+ const vector& pd2 = patchDisp[1][pointI];
- pointDisplacement_[pointI] =
- (1-s)*patchDisp[0][pointI]
- + s*patchDisp[1][pointI];
+ pointDisplacement_[pointI] = (1 - s)*pd1 + s*pd2;
+ }
}
}
+ else
+ {
+ 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];
+ }
+ }
+ }
}
diff --git a/src/meshTools/twoDPointCorrector/twoDPointCorrector.C b/src/meshTools/twoDPointCorrector/twoDPointCorrector.C
index ad0413e61b22fff214e03a3fc9a4f95b499255f8..798d0f8899b5a107061d0df3d96b89904b7e8030 100644
--- a/src/meshTools/twoDPointCorrector/twoDPointCorrector.C
+++ b/src/meshTools/twoDPointCorrector/twoDPointCorrector.C
@@ -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;
+
+ const wedgePolyPatch& wp =
+ refCast<const wedgePolyPatch>(patches[patchI]);
+
+ pn = wp.centreNormal();
- pn = refCast<const wedgePolyPatch>(patches[patchI]).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
-
// ************************************************************************* //
diff --git a/src/meshTools/twoDPointCorrector/twoDPointCorrector.H b/src/meshTools/twoDPointCorrector/twoDPointCorrector.H
index 3cdfc53e8505865635f7e313a410af6292368346..c30a0282a0fd297110efd264a34857142a86c3a3 100644
--- a/src/meshTools/twoDPointCorrector/twoDPointCorrector.H
+++ b/src/meshTools/twoDPointCorrector/twoDPointCorrector.H
@@ -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
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
index 23e2c35ac6a46869eb47e6599655f5f5252de6e7..4da01ef285f9e4e85093383896c0c75966f3fd5f 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
@@ -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;