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

ENH: cyclics: transformPosition for lagrangian positions. Rotational transform calculation

Browse Source 
For rotational:
- 180 degree transformation now correct
- use furthest away face to determine rotation
- do not use calculated faceNormals at all
This commit is contained in:
mattijs
2011-03-17 15:56:46 +00:00
parent feb65cd3b9
commit fbfac5e83c
11 changed files with 252 additions and 194 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H
View File

@ -254,7 +254,10 @@ public:
}
//- Transform a patch-based position from other side to this side
virtual void transformPosition(pointField& l) const = 0;
virtual void transformPosition(pointField&) const = 0;
//- Transform a patch-based position from other side to this side
virtual void transformPosition(point&, const label facei) const = 0;
//- Are the planes separated.
virtual bool separated() const
@ -297,12 +300,12 @@ public:
virtual void calcGeometry
(
const primitivePatch& referPatch,
const UList<point>& thisCtrs,
const UList<point>& thisAreas,
const UList<point>& thisCc,
const UList<point>& nbrCtrs,
const UList<point>& nbrAreas,
const UList<point>& nbrCc
const pointField& thisCtrs,
const vectorField& thisAreas,
const pointField& thisCc,
const pointField& nbrCtrs,
const vectorField& nbrAreas,
const pointField& nbrCc
) = 0;
//- Initialize ordering for primitivePatch. Does not

294
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
View File

@ -77,7 +77,6 @@ void Foam::cyclicPolyPatch::calcTransforms()
if (size())
{
// Half0
const cyclicPolyPatch& half0 = *this;
vectorField half0Areas(half0.size());
forAll(half0, facei)
@ -108,10 +107,10 @@ void Foam::cyclicPolyPatch::calcTransforms()
void Foam::cyclicPolyPatch::calcTransforms
(
const primitivePatch& half0,
const UList<point>& half0Ctrs,
const UList<point>& half0Areas,
const UList<point>& half1Ctrs,
const UList<point>& half1Areas
const pointField& half0Ctrs,
const vectorField& half0Areas,
const pointField& half1Ctrs,
const vectorField& half1Areas
)
{
if (debug && owner())
@ -165,23 +164,9 @@ void Foam::cyclicPolyPatch::calcTransforms
if (half0Ctrs.size() > 0)
{
scalarField half0Tols
(
calcFaceTol
(
half0,
half0.points(),
static_cast<const pointField&>(half0Ctrs)
)
);
vectorField half0Normals(half0Areas.size());
vectorField half1Normals(half1Areas.size());
//- Additional warning about faces non-aligned with rotation axis
//scalar maxCos = -GREAT;
//label maxFacei = -1;
forAll(half0, facei)
{
scalar magSf = mag(half0Areas[facei]);
@ -221,75 +206,74 @@ void Foam::cyclicPolyPatch::calcTransforms
{
half0Normals[facei] = half0Areas[facei] / magSf;
half1Normals[facei] = half1Areas[facei] / nbrMagSf;
//if (transform_ == ROTATIONAL)
//{
// scalar cos = mag(half0Normals[facei] & rotationAxis_);
// if (cos > maxCos)
// {
// maxCos = cos;
// maxFacei = facei;
// }
//}
}
}
//if (maxCos > sqrt(SMALL))
//{
// WarningIn
// (
// "cyclicPolyPatch::calcTransforms()"
// ) << "on patch " << name()
// << " face:" << maxFacei << " fc:" << half0Ctrs[maxFacei]
// << " is not perpendicular to the rotationAxis." << endl
// << "This will cause problems with topology changes." << endl
// << "rotation axis : " << rotationAxis_ << endl
// << "face normal : " << half0Normals[maxFacei] << endl
// << "cosine of angle : " << maxCos << endl;
//}
// Calculate transformation tensors
calcTransformTensors
(
static_cast<const pointField&>(half0Ctrs),
static_cast<const pointField&>(half1Ctrs),
half0Normals,
half1Normals,
half0Tols,
matchTol,
transform_
);
if (transform_ == ROTATIONAL && !parallel() && forwardT().size() > 1)
if (transform_ == ROTATIONAL)
{
// Get index of maximum area face to minimise truncation errors.
label max0I = findMaxArea(half0.points(), half0);
// Calculate using the given rotation axis and centre. Do not
// use calculated normals.
label face0 = getConsistentRotationFace(half0Ctrs);
label face1 = face0;
const tensor fwdT = forwardT()[max0I];
const_cast<tensorField&>(forwardT()) = tensorField(1, fwdT);
const tensor revT = reverseT()[max0I];
const_cast<tensorField&>(reverseT()) = tensorField(1, revT);
const bool coll = collocated()[max0I];
const_cast<boolList&>(collocated()).setSize(1);
const_cast<boolList&>(collocated())[0] = coll;
vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
n0 /= mag(n0) + VSMALL;
n1 /= mag(n1) + VSMALL;
WarningIn
if (debug)
{
Pout<< "cyclicPolyPatch::calcTransforms :"
<< " Specified rotation :"
<< " n0:" << n0 << " n1:" << n1 << endl;
}
// Extended tensor from two local coordinate systems calculated
// using normal and rotation axis
const tensor E0
(
"cyclicPolyPatch::calcTransforms\n"
" (\n"
" const primitivePatch&,\n"
" const UList<point>&,\n"
" const UList<point>&,\n"
" const UList<point>&,\n"
" const UList<point>&\n"
" )"
) << "For patch " << name()
<< " calculated non-uniform transform tensor even though"
<< " the transform type is " << transformTypeNames[transform_]
<< "." << nl
<< " Setting the transformation tensor to be a uniform"
<< " rotation calculated from face " << max0I
<< endl;
rotationAxis_,
(n0 ^ rotationAxis_),
n0
);
const tensor E1
(
rotationAxis_,
(-n1 ^ rotationAxis_),
-n1
);
const tensor revT(E1.T() & E0);
const_cast<tensorField&>(forwardT()) = tensorField(1, revT.T());
const_cast<tensorField&>(reverseT()) = tensorField(1, revT);
const_cast<vectorField&>(separation()).setSize(0);
const_cast<boolList&>(collocated()) = boolList(1, false);
}
else
{
scalarField half0Tols
(
calcFaceTol
(
half0,
half0.points(),
static_cast<const pointField&>(half0Ctrs)
)
);
calcTransformTensors
(
static_cast<const pointField&>(half0Ctrs),
static_cast<const pointField&>(half1Ctrs),
half0Normals,
half1Normals,
half0Tols,
matchTol,
transform_
);
}
}
}
@ -333,14 +317,39 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
<< " n0:" << n0 << " n1:" << n1 << endl;
}
// Rotation (around origin)
const tensor reverseT(rotationTensor(n0, -n1));
// Extended tensor from two local coordinate systems calculated
// using normal and rotation axis
const tensor E0
(
rotationAxis_,
(n0 ^ rotationAxis_),
n0
);
const tensor E1
(
rotationAxis_,
(-n1 ^ rotationAxis_),
-n1
);
const tensor revT(E1.T() & E0);
// Rotation
forAll(half0Ctrs, faceI)
{
half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
half0Ctrs[faceI] =
Foam::transform
(
revT,
half0Ctrs[faceI] - rotationCentre_
)
+ rotationCentre_;
anchors0[faceI] =
Foam::transform
(
revT,
anchors0[faceI] - rotationCentre_
)
+ rotationCentre_;
}
break;
@ -387,19 +396,19 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
}
// Rotation (around origin)
const tensor reverseT(rotationTensor(n0, -n1));
const tensor revT(rotationTensor(n0, -n1));
// Rotation
forAll(half0Ctrs, faceI)
{
half0Ctrs[faceI] = Foam::transform
(
reverseT,
revT,
half0Ctrs[faceI]
);
anchors0[faceI] = Foam::transform
(
reverseT,
revT,
anchors0[faceI]
);
}
@ -437,41 +446,20 @@ Foam::label Foam::cyclicPolyPatch::getConsistentRotationFace
const pointField& faceCentres
) const
{
const scalarField magRadSqr
(
magSqr((faceCentres - rotationCentre_) ^ rotationAxis_)
);
scalarField axisLen
(
(faceCentres - rotationCentre_) & rotationAxis_
);
axisLen -= min(axisLen);
const scalarField magLenSqr
(
magRadSqr + axisLen*axisLen
);
// Determine a face furthest away from the axis
label rotFace = -1;
scalar maxMagLenSqr = -GREAT;
scalar maxMagRadSqr = -GREAT;
forAll(faceCentres, i)
{
if (magLenSqr[i] >= maxMagLenSqr)
{
if (magRadSqr[i] > maxMagRadSqr)
{
rotFace = i;
maxMagLenSqr = magLenSqr[i];
maxMagRadSqr = magRadSqr[i];
}
}
}
const scalarField magRadSqr =
magSqr((faceCentres - rotationCentre_) ^ rotationAxis_);
label rotFace = findMax(magRadSqr);
if (debug)
{
Pout<< "getConsistentRotationFace(const pointField&)" << nl
<< " rotFace = " << rotFace << nl
<< " point = " << faceCentres[rotFace] << endl;
Info<< "getConsistentRotationFace(const pointField&)" << nl
<< " rotFace = " << rotFace << nl
<< " point = " << faceCentres[rotFace] << nl
<< " distance = " << Foam::sqrt(magRadSqr[rotFace])
<< endl;
}
return rotFace;
@ -586,6 +574,17 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
{
dict.lookup("rotationAxis") >> rotationAxis_;
dict.lookup("rotationCentre") >> rotationCentre_;
scalar magRot = mag(rotationAxis_);
if (magRot < SMALL)
{
FatalIOErrorIn("cyclicPolyPatch::cyclicPolyPatch(..)", dict)
<< "Illegal rotationAxis " << rotationAxis_ << endl
<< "Please supply a non-zero vector."
<< exit(FatalIOError);
}
rotationAxis_ /= magRot;
break;
}
case TRANSLATIONAL:
@ -730,7 +729,16 @@ void Foam::cyclicPolyPatch::transformPosition(pointField& l) const
{
if (!parallel())
{
l = Foam::transform(forwardT(), l);
if (transform_ == ROTATIONAL)
{
l =
Foam::transform(forwardT(), l-rotationCentre_)
+ rotationCentre_;
}
else
{
l = Foam::transform(forwardT(), l);
}
}
else if (separated())
{
@ -750,6 +758,40 @@ void Foam::cyclicPolyPatch::transformPosition(pointField& l) const
}
void Foam::cyclicPolyPatch::transformPosition(point& l, const label facei) const
{
if (!parallel())
{
const tensor& T =
(
forwardT().size() == 1
? forwardT()[0]
: forwardT()[facei]
);
if (transform_ == ROTATIONAL)
{
l = Foam::transform(T, l-rotationCentre_) + rotationCentre_;
}
else
{
l = Foam::transform(T, l);
}
}
else if (separated())
{
const vector& s =
(
separation().size() == 1
? separation()[0]
: separation()[facei]
);
l -= s;
}
}
void Foam::cyclicPolyPatch::initGeometry(PstreamBuffers& pBufs)
{
polyPatch::initGeometry(pBufs);
@ -759,9 +801,9 @@ void Foam::cyclicPolyPatch::initGeometry(PstreamBuffers& pBufs)
void Foam::cyclicPolyPatch::initGeometry
(
const primitivePatch& referPatch,
UList<point>& nbrCtrs,
UList<point>& nbrAreas,
UList<point>& nbrCc
pointField& nbrCtrs,
vectorField& nbrAreas,
pointField& nbrCc
)
{}
@ -769,12 +811,12 @@ void Foam::cyclicPolyPatch::initGeometry
void Foam::cyclicPolyPatch::calcGeometry
(
const primitivePatch& referPatch,
const UList<point>& thisCtrs,
const UList<point>& thisAreas,
const UList<point>& thisCc,
const UList<point>& nbrCtrs,
const UList<point>& nbrAreas,
const UList<point>& nbrCc
const pointField& thisCtrs,
const vectorField& thisAreas,
const pointField& thisCc,
const pointField& nbrCtrs,
const vectorField& nbrAreas,
const pointField& nbrCc
)
{
calcTransforms

29
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
View File

@ -110,10 +110,10 @@ class cyclicPolyPatch
void calcTransforms
(
const primitivePatch& half0,
const UList<point>& half0Ctrs,
const UList<point>& half0Areas,
const UList<point>& half1Ctrs,
const UList<point>& half1Areas
const pointField& half0Ctrs,
const vectorField& half0Areas,
const pointField& half1Ctrs,
const vectorField& half1Areas
);
// Face ordering
@ -149,9 +149,9 @@ protected:
virtual void initGeometry
(
const primitivePatch& referPatch,
UList<point>& nbrCtrs,
UList<point>& nbrAreas,
UList<point>& nbrCc
pointField& nbrCtrs,
vectorField& nbrAreas,
pointField& nbrCc
);
//- Calculate the patch geometry
@ -161,12 +161,12 @@ protected:
virtual void calcGeometry
(
const primitivePatch& referPatch,
const UList<point>& thisCtrs,
const UList<point>& thisAreas,
const UList<point>& thisCc,
const UList<point>& nbrCtrs,
const UList<point>& nbrAreas,
const UList<point>& nbrCc
const pointField& thisCtrs,
const vectorField& thisAreas,
const pointField& thisCc,
const pointField& nbrCtrs,
const vectorField& nbrAreas,
const pointField& nbrCc
);
//- Initialise the patches for moving points
@ -342,6 +342,9 @@ public:
//- Transform a patch-based position from other side to this side
virtual void transformPosition(pointField& l) const;
//- Transform a patch-based position from other side to this side
virtual void transformPosition(point&, const label facei) const;
// Transformation

14
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -711,12 +711,12 @@ void Foam::oldCyclicPolyPatch::initGeometry(PstreamBuffers& pBufs)
void Foam::oldCyclicPolyPatch::calcGeometry
(
const primitivePatch& referPatch,
const UList<point>& thisCtrs,
const UList<point>& thisAreas,
const UList<point>& thisCc,
const UList<point>& nbrCtrs,
const UList<point>& nbrAreas,
const UList<point>& nbrCc
const pointField& thisCtrs,
const vectorField& thisAreas,
const pointField& thisCc,
const pointField& nbrCtrs,
const vectorField& nbrAreas,
const pointField& nbrCc
)
{}

20
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -269,16 +269,22 @@ public:
notImplemented("transformPosition(pointField&)");
}
//- Transform a patch-based position from other side to this side
virtual void transformPosition(point&, const label facei) const
{
notImplemented("transformPosition(point&, const label)");
}
//- Calculate the patch geometry
virtual void calcGeometry
(
const primitivePatch& referPatch,
const UList<point>& thisCtrs,
const UList<point>& thisAreas,
const UList<point>& thisCc,
const UList<point>& nbrCtrs,
const UList<point>& nbrAreas,
const UList<point>& nbrCc
const pointField& thisCtrs,
const vectorField& thisAreas,
const pointField& thisCc,
const pointField& nbrCtrs,
const vectorField& nbrAreas,
const pointField& nbrCc
);
//- Initialize ordering for primitivePatch. Does not

18
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -93,12 +93,12 @@ protected:
virtual void calcGeometry
(
const primitivePatch& referPatch,
const UList<point>& thisCtrs,
const UList<point>& thisAreas,
const UList<point>& thisCc,
const UList<point>& nbrCtrs,
const UList<point>& nbrAreas,
const UList<point>& nbrCc
const pointField& thisCtrs,
const vectorField& thisAreas,
const pointField& thisCc,
const pointField& nbrCtrs,
const vectorField& nbrAreas,
const pointField& nbrCc
)
{
notImplemented("processorPolyPatch::calcGeometry(..)");
@ -300,6 +300,10 @@ public:
virtual void transformPosition(pointField& l) const
{}
//- Transform a patch-based position from other side to this side
virtual void transformPosition(point&, const label facei) const
{}
//- Initialize ordering for primitivePatch. Does not
// refer to *this (except for name() and type() etc.)
virtual void initOrder(PstreamBuffers&, const primitivePatch&) const;

18
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.H
View File

@ -81,12 +81,12 @@ protected:
virtual void calcGeometry
(
const primitivePatch& referPatch,
const UList<point>& thisCtrs,
const UList<point>& thisAreas,
const UList<point>& thisCc,
const UList<point>& nbrCtrs,
const UList<point>& nbrAreas,
const UList<point>& nbrCc
const pointField& thisCtrs,
const vectorField& thisAreas,
const pointField& thisCc,
const pointField& nbrCtrs,
const vectorField& nbrAreas,
const pointField& nbrCc
)
{
notImplemented("processorCyclicPolyPatch::calcGeometry(..)");
@ -279,6 +279,12 @@ public:
referPatch().transformPosition(l);
}
//- Transform a patch-based position from other side to this side
virtual void transformPosition(point& l, const label facei) const
{
referPatch().transformPosition(l, facei);
}
//- Are the planes separated.
virtual bool separated() const
{

3
src/finiteVolume/Make/files
View File

@ -114,7 +114,8 @@ derivedFvPatchFields = $(fvPatchFields)/derived
$(derivedFvPatchFields)/activeBaffleVelocity/activeBaffleVelocityFvPatchVectorField.C
$(derivedFvPatchFields)/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
$(derivedFvPatchFields)/advective/advectiveFvPatchFields.C
$(derivedFvPatchFields)/codedFixedValue/codedFixedValueFvPatchScalarField.C
/* $(derivedFvPatchFields)/codedFixedValue/codedFixedValueFvPatchScalarField.C */
$(derivedFvPatchFields)/codedFixedValue/codedFixedValueFvPatchFields.C
$(derivedFvPatchFields)/directMappedField/directMappedFieldFvPatchFields.C
$(derivedFvPatchFields)/directMappedFixedInternalValue/directMappedFixedInternalValueFvPatchFields.C
$(derivedFvPatchFields)/directMappedFixedPushedInternalValue/directMappedFixedPushedInternalValueFvPatchFields.C

6
src/lagrangian/basic/particle/particle.C
View File

@ -116,12 +116,6 @@ Foam::particle::particle(const particle& p, const polyMesh& mesh)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::particle::transformPosition(const tensor& T)
{
position_ = transform(T, position_);
}
void Foam::particle::transformProperties(const tensor&)
{}

4
src/lagrangian/basic/particle/particle.H
View File

@ -491,10 +491,6 @@ public:
// Transformations
//- Transform the position the particle
// according to the given transformation tensor
virtual void transformPosition(const tensor& T);
//- Transform the physical properties of the particle
// according to the given transformation tensor
virtual void transformProperties(const tensor& T);

23
src/lagrangian/basic/particle/particleTemplates.C
View File

@ -52,11 +52,15 @@ void Foam::particle::correctAfterParallelTransfer
TrackData& td
)
{
const processorPolyPatch& ppp =
refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
const coupledPolyPatch& ppp =
refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchI]);
cellI_ = ppp.faceCells()[faceI_];
// Have patch transform the position
ppp.transformPosition(position_, faceI_);
// Transform the properties
if (!ppp.parallel())
{
const tensor& T =
@ -65,8 +69,6 @@ void Foam::particle::correctAfterParallelTransfer
? ppp.forwardT()[0]
: ppp.forwardT()[faceI_]
);
transformPosition(T);
transformProperties(T);
}
else if (ppp.separated())
@ -77,7 +79,6 @@ void Foam::particle::correctAfterParallelTransfer
? ppp.separation()[0]
: ppp.separation()[faceI_]
);
position_ -= s;
transformProperties(-s);
}
@ -958,19 +959,22 @@ void Foam::particle::hitCyclicPatch
tetPtI_ = mesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
const cyclicPolyPatch& receiveCpp = cpp.neighbPatch();
label patchFacei = receiveCpp.whichFace(faceI_);
// Now the particle is on the receiving side
// Have patch transform the position
receiveCpp.transformPosition(position_, patchFacei);
// Transform the properties
if (!receiveCpp.parallel())
{
const tensor& T =
(
receiveCpp.forwardT().size() == 1
? receiveCpp.forwardT()[0]
: receiveCpp.forwardT()[receiveCpp.whichFace(faceI_)]
: receiveCpp.forwardT()[patchFacei]
);
transformPosition(T);
transformProperties(T);
}
else if (receiveCpp.separated())
@ -979,9 +983,8 @@ void Foam::particle::hitCyclicPatch
(
(receiveCpp.separation().size() == 1)
? receiveCpp.separation()[0]
: receiveCpp.separation()[receiveCpp.whichFace(faceI_)]
: receiveCpp.separation()[patchFacei]
);
position_ -= s;
transformProperties(-s);
}
}