ENH: cyclicAMI: consistent parallel transformations

This commit is contained in:
mattijs
2011-11-11 13:47:26 +00:00
parent ef4a575788
commit fac508c56d
4 changed files with 123 additions and 84 deletions

View File

@ -251,6 +251,7 @@ void Foam::cyclicPolyPatch::calcTransforms
if (debug) if (debug)
{ {
Pout<< "cyclicPolyPatch::calcTransforms :" Pout<< "cyclicPolyPatch::calcTransforms :"
<< " patch:" << name()
<< " Max area error:" << 100*maxAreaDiff << "% at face:" << " Max area error:" << 100*maxAreaDiff << "% at face:"
<< maxAreaFacei << " at:" << half0Ctrs[maxAreaFacei] << maxAreaFacei << " at:" << half0Ctrs[maxAreaFacei]
<< " coupled face at:" << half1Ctrs[maxAreaFacei] << " coupled face at:" << half1Ctrs[maxAreaFacei]
@ -264,17 +265,15 @@ void Foam::cyclicPolyPatch::calcTransforms
{ {
// Calculate using the given rotation axis and centre. Do not // Calculate using the given rotation axis and centre. Do not
// use calculated normals. // use calculated normals.
label face0 = getConsistentRotationFace(half0Ctrs); vector n0 = findFaceMaxRadius(half0Ctrs);
label face1 = face0; vector n1 = findFaceMaxRadius(half1Ctrs);
vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
n0 /= mag(n0) + VSMALL; n0 /= mag(n0) + VSMALL;
n1 /= mag(n1) + VSMALL; n1 /= mag(n1) + VSMALL;
if (debug) if (debug)
{ {
Pout<< "cyclicPolyPatch::calcTransforms :" Pout<< "cyclicPolyPatch::calcTransforms :"
<< " patch:" << name()
<< " Specified rotation :" << " Specified rotation :"
<< " n0:" << n0 << " n1:" << n1 << endl; << " n0:" << n0 << " n1:" << n1 << endl;
} }
@ -330,6 +329,7 @@ void Foam::cyclicPolyPatch::calcTransforms
if (debug) if (debug)
{ {
Pout<< "cyclicPolyPatch::calcTransforms :" Pout<< "cyclicPolyPatch::calcTransforms :"
<< " patch:" << name()
<< " Specified separation vector : " << " Specified separation vector : "
<< separationVector_ << endl; << separationVector_ << endl;
} }
@ -426,17 +426,15 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
{ {
case ROTATIONAL: case ROTATIONAL:
{ {
label face0 = getConsistentRotationFace(half0Ctrs); vector n0 = findFaceMaxRadius(half0Ctrs);
label face1 = getConsistentRotationFace(half1Ctrs); vector n1 = findFaceMaxRadius(half1Ctrs);
vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
n0 /= mag(n0) + VSMALL; n0 /= mag(n0) + VSMALL;
n1 /= mag(n1) + VSMALL; n1 /= mag(n1) + VSMALL;
if (debug) if (debug)
{ {
Pout<< "cyclicPolyPatch::getCentresAndAnchors :" Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
<< " patch:" << name()
<< " Specified rotation :" << " Specified rotation :"
<< " n0:" << n0 << " n1:" << n1 << endl; << " n0:" << n0 << " n1:" << n1 << endl;
} }
@ -485,6 +483,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
if (debug) if (debug)
{ {
Pout<< "cyclicPolyPatch::getCentresAndAnchors :" Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
<< " patch:" << name()
<< "Specified translation : " << separationVector_ << "Specified translation : " << separationVector_
<< endl; << endl;
} }
@ -516,6 +515,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
if (debug) if (debug)
{ {
Pout<< "cyclicPolyPatch::getCentresAndAnchors :" Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
<< " patch:" << name()
<< " Detected rotation :" << " Detected rotation :"
<< " n0:" << n0 << " n1:" << n1 << endl; << " n0:" << n0 << " n1:" << n1 << endl;
} }
@ -548,6 +548,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
if (debug) if (debug)
{ {
Pout<< "cyclicPolyPatch::getCentresAndAnchors :" Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
<< " patch:" << name()
<< " Detected translation :" << " Detected translation :"
<< " n0:" << n0 << " n1:" << n1 << " n0:" << n0 << " n1:" << n1
<< " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl; << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
@ -566,30 +567,29 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
} }
Foam::label Foam::cyclicPolyPatch::getConsistentRotationFace Foam::vector Foam::cyclicPolyPatch::findFaceMaxRadius
( (
const pointField& faceCentres const pointField& faceCentres
) const ) const
{ {
// Determine a face furthest away from the axis // Determine a face furthest away from the axis
const scalarField magRadSqr const vectorField n((faceCentres - rotationCentre_) ^ rotationAxis_);
(
magSqr((faceCentres - rotationCentre_) ^ rotationAxis_)
);
label rotFace = findMax(magRadSqr); const scalarField magRadSqr(magSqr(n));
label faceI = findMax(magRadSqr);
if (debug) if (debug)
{ {
Info<< "getConsistentRotationFace(const pointField&)" << nl Info<< "findFaceMaxRadius(const pointField&) : patch: " << name() << nl
<< " rotFace = " << rotFace << nl << " rotFace = " << faceI << nl
<< " point = " << faceCentres[rotFace] << nl << " point = " << faceCentres[faceI] << nl
<< " distance = " << Foam::sqrt(magRadSqr[rotFace]) << " distance = " << Foam::sqrt(magRadSqr[faceI])
<< endl; << endl;
} }
return rotFace; return n[faceI];
} }

View File

@ -130,9 +130,8 @@ class cyclicPolyPatch
scalarField& tols scalarField& tols
) const; ) const;
//- For rotational cases, try to find a unique face on each side //- Return normal of face at max distance from rotation axis
// of the cyclic. vector findFaceMaxRadius(const pointField& faceCentres) const;
label getConsistentRotationFace(const pointField&) const;
protected: protected:

View File

@ -30,6 +30,7 @@ License
#include "Time.H" #include "Time.H"
#include "addToRunTimeSelectionTable.H" #include "addToRunTimeSelectionTable.H"
#include "faceAreaIntersect.H" #include "faceAreaIntersect.H"
#include "ops.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -44,63 +45,68 @@ namespace Foam
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
Foam::label Foam::cyclicAMIPolyPatch::findFaceMaxRadius Foam::vector Foam::cyclicAMIPolyPatch::findFaceMaxRadius
( (
const pointField& faceCentres const pointField& faceCentres
) const ) const
{ {
// Determine a face furthest away from the axis // Determine a face furthest away from the axis
const scalarField magRadSqr const vectorField n((faceCentres - rotationCentre_) ^ rotationAxis_);
(
magSqr((faceCentres - rotationCentre_) ^ rotationAxis_) const scalarField magRadSqr(magSqr(n));
);
label faceI = findMax(magRadSqr); label faceI = findMax(magRadSqr);
if (debug) if (debug)
{ {
Info<< "findFaceMaxRadius(const pointField&)" << nl Info<< "findFaceMaxRadius(const pointField&) : patch: " << name() << nl
<< " rotFace = " << faceI << nl << " rotFace = " << faceI << nl
<< " point = " << faceCentres[faceI] << nl << " point = " << faceCentres[faceI] << nl
<< " distance = " << Foam::sqrt(magRadSqr[faceI]) << " distance = " << Foam::sqrt(magRadSqr[faceI])
<< endl; << endl;
} }
return faceI; return n[faceI];
} }
// * * * * * * * * * * * Protecetd Member Functions * * * * * * * * * * * * // // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
void Foam::cyclicAMIPolyPatch::calcTransforms() void Foam::cyclicAMIPolyPatch::calcTransforms()
{ {
if (size()) // Half0
const cyclicAMIPolyPatch& half0 = *this;
vectorField half0Areas(half0.size());
forAll(half0, facei)
{ {
// Half0 half0Areas[facei] = half0[facei].normal(half0.points());
const cyclicAMIPolyPatch& half0 = *this; }
vectorField half0Areas(half0.size());
forAll(half0, facei)
{
half0Areas[facei] = half0[facei].normal(half0.points());
}
// Half1 // Half1
const cyclicAMIPolyPatch& half1 = neighbPatch(); const cyclicAMIPolyPatch& half1 = neighbPatch();
vectorField half1Areas(half1.size()); vectorField half1Areas(half1.size());
forAll(half1, facei) forAll(half1, facei)
{ {
half1Areas[facei] = half1[facei].normal(half1.points()); half1Areas[facei] = half1[facei].normal(half1.points());
} }
calcTransforms calcTransforms
( (
half0, half0,
half0.faceCentres(), half0.faceCentres(),
half0Areas, half0Areas,
half1.faceCentres(), half1.faceCentres(),
half1Areas half1Areas
); );
if (debug)
{
Pout<< "calcTransforms() : patch: " << name() << nl
<< " forwardT = " << forwardT() << nl
<< " reverseT = " << reverseT() << nl
<< " separation = " << separation() << nl
<< " collocated = " << collocated() << nl << endl;
} }
} }
@ -128,26 +134,30 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
// Calculate transformation tensors // Calculate transformation tensors
if ((half0Ctrs.size() <= 0) || (half1Ctrs.size() <= 0))
{
return;
}
switch (transform_) switch (transform_)
{ {
case ROTATIONAL: case ROTATIONAL:
{ {
label face0 = findFaceMaxRadius(half0Ctrs); point n0 = vector::zero;
label face1 = findFaceMaxRadius(half1Ctrs); if (half0Ctrs.size())
{
n0 = findFaceMaxRadius(half0Ctrs);
}
point n1 = vector::zero;
if (half1Ctrs.size())
{
n1 = -findFaceMaxRadius(half1Ctrs);
}
reduce(n0, maxMagSqrOp<point>());
reduce(n1, maxMagSqrOp<point>());
vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
n0 /= mag(n0) + VSMALL; n0 /= mag(n0) + VSMALL;
n1 /= mag(n1) + VSMALL; n1 /= mag(n1) + VSMALL;
if (debug) if (debug)
{ {
Pout<< "cyclicAMIPolyPatch::calcTransforms :" Pout<< "cyclicAMIPolyPatch::calcTransforms : patch:" << name()
<< " Specified rotation :" << " Specified rotation :"
<< " n0:" << n0 << " n1:" << n1 << endl; << " n0:" << n0 << " n1:" << n1 << endl;
} }
@ -178,8 +188,8 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
{ {
if (debug) if (debug)
{ {
Pout<< "cyclicAMIPolyPatch::calcTransforms :" Pout<< "cyclicAMIPolyPatch::calcTransforms : patch:" << name()
<< "Specified translation : " << separationVector_ << " Specified translation : " << separationVector_
<< endl; << endl;
} }
@ -198,7 +208,8 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
{ {
if (debug) if (debug)
{ {
Pout<< "Assuming cyclic AMI pairs are colocated" << endl; Pout<< " patch:" << name()
<< " Assuming cyclic AMI pairs are colocated" << endl;
} }
const_cast<tensorField&>(forwardT()).clear(); const_cast<tensorField&>(forwardT()).clear();
@ -232,7 +243,8 @@ void Foam::cyclicAMIPolyPatch::resetAMI() const
if (debug) if (debug)
{ {
OFstream os(name() + "_neighbourPatch-org.obj"); const Time& t = boundaryMesh().mesh().time();
OFstream os(t.path()/name() + "_neighbourPatch-org.obj");
meshTools::writeOBJ(os, neighbPatch().localFaces(), nbrPoints); meshTools::writeOBJ(os, neighbPatch().localFaces(), nbrPoints);
} }
@ -250,10 +262,11 @@ void Foam::cyclicAMIPolyPatch::resetAMI() const
if (debug) if (debug)
{ {
OFstream osN(name() + "_neighbourPatch-trans.obj"); const Time& t = boundaryMesh().mesh().time();
OFstream osN(t.path()/name() + "_neighbourPatch-trans.obj");
meshTools::writeOBJ(osN, nbrPatch0.localFaces(), nbrPoints); meshTools::writeOBJ(osN, nbrPatch0.localFaces(), nbrPoints);
OFstream osO(name() + "_ownerPatch.obj"); OFstream osO(t.path()/name() + "_ownerPatch.obj");
meshTools::writeOBJ(osO, this->localFaces(), localPoints()); meshTools::writeOBJ(osO, this->localFaces(), localPoints());
} }
@ -451,10 +464,10 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
coupledPolyPatch(pp, bm), coupledPolyPatch(pp, bm),
nbrPatchName_(pp.nbrPatchName_), nbrPatchName_(pp.nbrPatchName_),
nbrPatchID_(-1), nbrPatchID_(-1),
transform_(UNKNOWN), transform_(pp.transform_),
rotationAxis_(vector::zero), rotationAxis_(pp.rotationAxis_),
rotationCentre_(point::zero), rotationCentre_(pp.rotationCentre_),
separationVector_(vector::zero), separationVector_(pp.separationVector_),
AMIPtr_(NULL), AMIPtr_(NULL),
surfPtr_(NULL), surfPtr_(NULL),
surfDict_(dictionary::null) surfDict_(dictionary::null)
@ -477,10 +490,10 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
coupledPolyPatch(pp, bm, index, newSize, newStart), coupledPolyPatch(pp, bm, index, newSize, newStart),
nbrPatchName_(nbrPatchName), nbrPatchName_(nbrPatchName),
nbrPatchID_(-1), nbrPatchID_(-1),
transform_(UNKNOWN), transform_(pp.transform_),
rotationAxis_(vector::zero), rotationAxis_(pp.rotationAxis_),
rotationCentre_(point::zero), rotationCentre_(pp.rotationCentre_),
separationVector_(vector::zero), separationVector_(pp.separationVector_),
AMIPtr_(NULL), AMIPtr_(NULL),
surfPtr_(NULL), surfPtr_(NULL),
surfDict_(dictionary::null) surfDict_(dictionary::null)
@ -535,6 +548,23 @@ Foam::cyclicAMIPolyPatch::~cyclicAMIPolyPatch()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::cyclicAMIPolyPatch::coupled() const
{
if
(
Pstream::parRun()
|| (size() && neighbPatch().size())
)
{
return true;
}
else
{
return false;
}
}
Foam::label Foam::cyclicAMIPolyPatch::neighbPatchID() const Foam::label Foam::cyclicAMIPolyPatch::neighbPatchID() const
{ {
if (nbrPatchID_ == -1) if (nbrPatchID_ == -1)
@ -769,6 +799,7 @@ void Foam::cyclicAMIPolyPatch::write(Ostream& os) const
coupledPolyPatch::write(os); coupledPolyPatch::write(os);
os.writeKeyword("neighbourPatch") << nbrPatchName_ os.writeKeyword("neighbourPatch") << nbrPatchName_
<< token::END_STATEMENT << nl; << token::END_STATEMENT << nl;
switch (transform_) switch (transform_)
{ {
case ROTATIONAL: case ROTATIONAL:
@ -801,8 +832,12 @@ void Foam::cyclicAMIPolyPatch::write(Ostream& os) const
} }
} }
os.writeKeyword(surfDict_.dictName());
os << surfDict_; if (surfDict_ != dictionary::null)
{
os.writeKeyword(surfDict_.dictName());
os << surfDict_;
}
} }

View File

@ -97,8 +97,8 @@ private:
// Private Member Functions // Private Member Functions
//- Return face ID of face at max distance from rotation axis //- Return normal of face at max distance from rotation axis
label findFaceMaxRadius(const pointField& faceCentres) const; vector findFaceMaxRadius(const pointField& faceCentres) const;
void calcTransforms void calcTransforms
( (
@ -254,6 +254,11 @@ public:
// Access // Access
//- Return true only if is coupled. Note that for non-parallel
// operation of a decomposed case this can return the wrong
// result
virtual bool coupled() const;
//- Neighbour patch name //- Neighbour patch name
inline const word& neighbPatchName() const; inline const word& neighbPatchName() const;