translation vector on cyclicPolyPatch; adapted createPatch

This commit is contained in:
mattijs
2008-10-08 11:42:23 +01:00
parent 3feae1f3f6
commit 12bec33d04
5 changed files with 601 additions and 159 deletions

View File

@ -34,6 +34,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "cyclicPolyPatch.H"
#include "syncTools.H"
#include "argList.H"
#include "polyMesh.H"
@ -256,27 +257,6 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
);
}
// cycPatch.writeOBJ
// (
// prefix+cycPatch.name()+"_half0.obj",
// SubList<face>
// (
// cycPatch,
// halfSize
// ),
// cycPatch.points()
// );
// cycPatch.writeOBJ
// (
// prefix+cycPatch.name()+"_half1.obj",
// SubList<face>
// (
// cycPatch,
// halfSize,
// halfSize
// ),
// cycPatch.points()
// );
// Lines between corresponding face centres
OFstream str(prefix+cycPatch.name()+"_match.obj");
@ -289,7 +269,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
vertI++;
label nbrFaceI = halfSize + faceI;
const point& fc1 = mesh.faceCentres()[cycPatch.start()+nbrFaceI];
const point& fc1 =
mesh.faceCentres()[cycPatch.start()+nbrFaceI];
meshTools::writeOBJ(str, fc1);
vertI++;
@ -300,6 +281,247 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
}
void separateList
(
const vectorField& separation,
UList<vector>& field
)
{
if (separation.size() == 1)
{
// Single value for all.
forAll(field, i)
{
field[i] += separation[0];
}
}
else if (separation.size() == field.size())
{
forAll(field, i)
{
field[i] += separation[i];
}
}
else
{
FatalErrorIn
(
"separateList(const vectorField&, UList<vector>&)"
) << "Sizes of field and transformation not equal. field:"
<< field.size() << " transformation:" << separation.size()
<< abort(FatalError);
}
}
// Synchronise points on both sides of coupled boundaries.
template <class CombineOp>
void syncPoints
(
const polyMesh& mesh,
pointField& points,
const CombineOp& cop,
const point& nullValue
)
{
if (points.size() != mesh.nPoints())
{
FatalErrorIn
(
"syncPoints"
"(const polyMesh&, pointField&, const CombineOp&, const point&)"
) << "Number of values " << points.size()
<< " is not equal to the number of points in the mesh "
<< mesh.nPoints() << abort(FatalError);
}
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Is there any coupled patch with transformation?
bool hasTransformation = false;
if (Pstream::parRun())
{
// Send
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if
(
isA<processorPolyPatch>(pp)
&& pp.nPoints() > 0
&& refCast<const processorPolyPatch>(pp).owner()
)
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
// Get data per patchPoint in neighbouring point numbers.
pointField patchInfo(procPatch.nPoints(), nullValue);
const labelList& meshPts = procPatch.meshPoints();
const labelList& nbrPts = procPatch.neighbPoints();
forAll(nbrPts, pointI)
{
label nbrPointI = nbrPts[pointI];
if (nbrPointI >= 0 && nbrPointI < patchInfo.size())
{
patchInfo[nbrPointI] = points[meshPts[pointI]];
}
}
OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
toNbr << patchInfo;
}
}
// Receive and set.
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if
(
isA<processorPolyPatch>(pp)
&& pp.nPoints() > 0
&& !refCast<const processorPolyPatch>(pp).owner()
)
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
pointField nbrPatchInfo(procPatch.nPoints());
{
// We do not know the number of points on the other side
// so cannot use Pstream::read.
IPstream fromNbr
(
Pstream::blocking,
procPatch.neighbProcNo()
);
fromNbr >> nbrPatchInfo;
}
// Null any value which is not on neighbouring processor
nbrPatchInfo.setSize(procPatch.nPoints(), nullValue);
if (!procPatch.parallel())
{
hasTransformation = true;
transformList(procPatch.forwardT(), nbrPatchInfo);
}
else if (procPatch.separated())
{
hasTransformation = true;
separateList(-procPatch.separation(), nbrPatchInfo);
}
const labelList& meshPts = procPatch.meshPoints();
forAll(meshPts, pointI)
{
label meshPointI = meshPts[pointI];
points[meshPointI] = nbrPatchInfo[pointI];
}
}
}
}
// Do the cyclics.
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (isA<cyclicPolyPatch>(pp))
{
const cyclicPolyPatch& cycPatch =
refCast<const cyclicPolyPatch>(pp);
const edgeList& coupledPoints = cycPatch.coupledPoints();
const labelList& meshPts = cycPatch.meshPoints();
pointField half0Values(coupledPoints.size());
forAll(coupledPoints, i)
{
const edge& e = coupledPoints[i];
label point0 = meshPts[e[0]];
half0Values[i] = points[point0];
}
if (!cycPatch.parallel())
{
hasTransformation = true;
transformList(cycPatch.reverseT(), half0Values);
}
else if (cycPatch.separated())
{
hasTransformation = true;
const vectorField& v = cycPatch.coupledPolyPatch::separation();
separateList(v, half0Values);
}
forAll(coupledPoints, i)
{
const edge& e = coupledPoints[i];
label point1 = meshPts[e[1]];
points[point1] = half0Values[i];
}
}
}
//- Note: hasTransformation is only used for warning messages so
// reduction not strictly nessecary.
//reduce(hasTransformation, orOp<bool>());
// Synchronize multiple shared points.
const globalMeshData& pd = mesh.globalData();
if (pd.nGlobalPoints() > 0)
{
if (hasTransformation)
{
WarningIn
(
"syncPoints"
"(const polyMesh&, pointField&, const CombineOp&, const point&)"
) << "There are decomposed cyclics in this mesh with"
<< " transformations." << endl
<< "This is not supported. The result will be incorrect"
<< endl;
}
// Values on shared points.
pointField sharedPts(pd.nGlobalPoints(), nullValue);
forAll(pd.sharedPointLabels(), i)
{
label meshPointI = pd.sharedPointLabels()[i];
// Fill my entries in the shared points
sharedPts[pd.sharedPointAddr()[i]] = points[meshPointI];
}
// Combine on master.
Pstream::listCombineGather(sharedPts, cop);
Pstream::listCombineScatter(sharedPts);
// Now we will all have the same information. Merge it back with
// my local information.
forAll(pd.sharedPointLabels(), i)
{
label meshPointI = pd.sharedPointLabels()[i];
points[meshPointI] = sharedPts[pd.sharedPointAddr()[i]];
}
}
}
// Main program:
int main(int argc, char *argv[])
@ -393,25 +615,28 @@ int main(int argc, char *argv[])
label destPatchI = patches.findPatchID(patchName);
word patchType(dict.lookup("type"));
if (destPatchI == -1)
{
dictionary patchDict(dict.subDict("dictionary"));
destPatchI = allPatches.size();
Info<< "Adding new patch " << patchName
<< " of type " << patchType
<< " as patch " << destPatchI << endl;
<< " as patch " << destPatchI
<< " from " << patchDict << endl;
patchDict.remove("nFaces");
patchDict.add("nFaces", 0);
patchDict.remove("startFace");
patchDict.add("startFace", startFaceI);
// Add an empty patch.
allPatches.append
(
polyPatch::New
(
patchType,
patchName,
0, // size
startFaceI, // start
patchDict,
destPatchI,
patches
).ptr()
@ -557,16 +782,100 @@ int main(int argc, char *argv[])
}
else
{
Info<< "Synchronising points." << nl << endl;
// This is a bit tricky. Both normal and position might be out and
// current separation also includes the normal
// ( separation_ = (nf&(Cr - Cf))*nf ).
// For processor patches:
// - disallow multiple separation/transformation. This basically
// excludes decomposed cyclics. Use the (probably 0) separation
// to align the points.
// For cyclic patches:
// - for separated ones use our own recalculated offset vector
// - for rotational ones use current one.
forAll(mesh.boundaryMesh(), patchI)
{
const polyPatch& pp = mesh.boundaryMesh()[patchI];
if (pp.size() && isA<coupledPolyPatch>(pp))
{
const coupledPolyPatch& cpp =
refCast<const coupledPolyPatch>(pp);
if (cpp.separated())
{
Info<< "On coupled patch " << pp.name()
<< " separation[0] was "
<< cpp.separation()[0] << endl;
if (isA<cyclicPolyPatch>(pp))
{
const cyclicPolyPatch& cycpp =
refCast<const cyclicPolyPatch>(pp);
if (cycpp.transform() == cyclicPolyPatch::TRANSLATIONAL)
{
Info<< "On cyclic translation patch " << pp.name()
<< " forcing uniform separation of "
<< cycpp.separationVector() << endl;
const_cast<vectorField&>(cpp.separation()) =
pointField(1, cycpp.separationVector());
}
else
{
const_cast<vectorField&>(cpp.separation()) =
pointField
(
1,
pp[pp.size()/2].centre(mesh.points())
- pp[0].centre(mesh.points())
);
}
}
else
{
const_cast<vectorField&>(cpp.separation())
.setSize(1);
}
Info<< "On coupled patch " << pp.name()
<< " forcing uniform separation of "
<< cpp.separation() << endl;
}
else if (!cpp.parallel())
{
Info<< "On coupled patch " << pp.name()
<< " forcing uniform rotation of "
<< cpp.forwardT()[0] << endl;
const_cast<tensorField&>
(
cpp.forwardT()
).setSize(1);
const_cast<tensorField&>
(
cpp.reverseT()
).setSize(1);
Info<< "On coupled patch " << pp.name()
<< " forcing uniform rotation of "
<< cpp.forwardT() << endl;
}
}
}
Info<< "Synchronising points." << endl;
pointField newPoints(mesh.points());
syncTools::syncPointList
syncPoints
(
mesh,
newPoints,
nearestEqOp(), // cop
point(GREAT, GREAT, GREAT), // nullValue
true // applySeparation
nearestEqOp(),
point(GREAT, GREAT, GREAT)
);
scalarField diff(mag(newPoints-mesh.points()));

View File

@ -1,58 +1,77 @@
/*--------------------------------*- C++ -*----------------------------------*\
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\ / O peration | Version: 1.0 |
| \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
root "";
case "";
instance "system";
local "";
class dictionary;
object createPatchDict;
object createPatcheDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Tolerance used in matching faces. Absolute tolerance is span of
// face times this factor.
matchTolerance 1E-3;
matchTolerance 1E-6;
// Do a synchronisation of coupled points.
pointSync true;
// Patches to create.
// If no patches does a coupled point and face synchronisation anyway.
patches
(
{
// Name of new patch
name leftRight0;
name sidePatches;
// Type of new patch
// Dictionary for new patch
dictionary
{
type cyclic;
// Optional: used when matching and synchronising points.
//transform translational;
//separationVector (-2289 0 0);
}
// How to construct: either 'patches' or 'set'
constructFrom patches;
// If constructFrom = patches : names of patches
patches (half0 half1);
//patches (periodic-1 periodic-2);
patches (outlet-side1 outlet-side2);
// If constructFrom = set : name of faceSet
set f0;
}
{
name bottom;
type patch;
constructFrom set;
patches (half0 half1);
set bottomFaces;
}
//{
// name bottom;
// // Dictionary for new patch
// dictionary
// {
// type patch;
// }
//
// constructFrom set;
//
// patches (half0 half1);
//
// set bottomFaces;
//}
);

View File

@ -263,6 +263,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl
<< " (half)size:" << Cf.size() << nl
<< " absTol:" << absTol << nl
//<< " smallDist:" << smallDist << nl
<< " sum(mag(nf & nr)):" << sum(mag(nf & nr)) << endl;
}
@ -316,6 +317,13 @@ void Foam::coupledPolyPatch::calcTransformTensors
{
forwardT_.setSize(1);
reverseT_.setSize(1);
if (debug)
{
Pout<< " rotation " << sum(mag(forwardT_ - forwardT_[0]))
<< " more than local tolerance " << error
<< ". Assuming uniform rotation." << endl;
}
}
}
else
@ -384,7 +392,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
if (debug)
{
Pout<< " separation_:" << separation_ << nl
Pout<< " separation_:" << separation_.size() << nl
<< " forwardT size:" << forwardT_.size() << endl;
}
}

View File

@ -127,6 +127,22 @@ void Foam::cyclicPolyPatch::calcTransforms()
Pout<< "cyclicPolyPatch::calcTransforms : Writing half1"
<< " faces to OBJ file " << nm1 << endl;
writeOBJ(nm1, half1, half1.points());
OFstream str(casePath/name()+"_half0_to_half1.obj");
label vertI = 0;
Pout<< "cyclicPolyPatch::calcTransforms :"
<< " Writing coupled face centres as lines to " << str.name()
<< endl;
forAll(half0Ctrs, i)
{
const point& p0 = half0Ctrs[i];
str << "v " << p0.x() << ' ' << p0.y() << ' ' << p0.z() << nl;
vertI++;
const point& p1 = half1Ctrs[i];
str << "v " << p1.x() << ' ' << p1.y() << ' ' << p1.z() << nl;
vertI++;
str << "l " << vertI-1 << ' ' << vertI << nl;
}
}
vectorField half0Normals(half0.size());
@ -397,8 +413,6 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
anchors0 = getAnchorPoints(half0Faces, pp.points());
half1Ctrs = calcFaceCentres(half1Faces, pp.points());
vector n0 = vector::zero;
vector n1 = vector::zero;
switch (transform_)
{
case ROTATIONAL:
@ -406,34 +420,15 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
label face0 = getConsistentRotationFace(half0Ctrs);
label face1 = getConsistentRotationFace(half1Ctrs);
n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
n0 /= mag(n0) + VSMALL;
n1 /= mag(n1) + VSMALL;
break;
}
default:
{
// Assumes that cyclic is planar. This is also the initial
// condition for patches without faces.
// Determine the face with max area on both halves. These
// two faces are used to determine the transformation tensors
label max0I = findMaxArea(pp.points(), half0Faces);
n0 = half0Faces[max0I].normal(pp.points());
n0 /= mag(n0) + VSMALL;
label max1I = findMaxArea(pp.points(), half1Faces);
n1 = half1Faces[max1I].normal(pp.points());
n1 /= mag(n1) + VSMALL;
}
}
if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
{
if (debug)
{
Pout<< "cyclicPolyPatch::getCentresAndAnchors : Rotation :"
Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
<< " Specified rotation :"
<< " n0:" << n0 << " n1:" << n1 << endl;
}
@ -446,6 +441,66 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
}
break;
}
//- Problem: usually specified translation is not accurate enough
//- to get proper match so keep automatic determination over here.
//case TRANSLATIONAL:
//{
// // Transform 0 points.
//
// if (debug)
// {
// Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
// << "Specified translation : " << separationVector_ << endl;
// }
//
// half0Ctrs += separationVector_;
// anchors0 += separationVector_;
// break;
//}
default:
{
// Assumes that cyclic is planar. This is also the initial
// condition for patches without faces.
// Determine the face with max area on both halves. These
// two faces are used to determine the transformation tensors
label max0I = findMaxArea(pp.points(), half0Faces);
vector n0 = half0Faces[max0I].normal(pp.points());
n0 /= mag(n0) + VSMALL;
label max1I = findMaxArea(pp.points(), half1Faces);
vector n1 = half1Faces[max1I].normal(pp.points());
n1 /= mag(n1) + VSMALL;
if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
{
if (debug)
{
Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
<< " Detected rotation :"
<< " n0:" << n0 << " n1:" << n1 << endl;
}
// Rotation (around origin)
const tensor reverseT(rotationTensor(n0, -n1));
// Rotation
forAll(half0Ctrs, faceI)
{
half0Ctrs[faceI] = Foam::transform
(
reverseT,
half0Ctrs[faceI]
);
anchors0[faceI] = Foam::transform
(
reverseT,
anchors0[faceI]
);
}
}
else
{
@ -461,7 +516,8 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
if (debug)
{
Pout<< "cyclicPolyPatch::getCentresAndAnchors : Translation :"
Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
<< " Detected translation :"
<< " n0:" << n0 << " n1:" << n1
<< " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
}
@ -469,6 +525,10 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
half0Ctrs += ctr1 - ctr0;
anchors0 += ctr1 - ctr0;
}
break;
}
}
// Calculate typical distance per face
tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs);
@ -615,7 +675,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(0.9),
transform_(UNKNOWN),
rotationAxis_(vector::zero),
rotationCentre_(point::zero)
rotationCentre_(point::zero),
separationVector_(vector::zero)
{
calcTransforms();
}
@ -635,7 +696,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(0.9),
transform_(UNKNOWN),
rotationAxis_(vector::zero),
rotationCentre_(point::zero)
rotationCentre_(point::zero),
separationVector_(vector::zero)
{
dict.readIfPresent("featureCos", featureCos_);
@ -650,9 +712,14 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
dict.lookup("rotationCentre") >> rotationCentre_;
break;
}
case TRANSLATIONAL:
{
dict.lookup("separationVector") >> separationVector_;
break;
}
default:
{
// no additioanl info required
// no additional info required
}
}
}
@ -673,7 +740,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(pp.featureCos_),
transform_(pp.transform_),
rotationAxis_(pp.rotationAxis_),
rotationCentre_(pp.rotationCentre_)
rotationCentre_(pp.rotationCentre_),
separationVector_(pp.separationVector_)
{
calcTransforms();
}
@ -694,7 +762,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(pp.featureCos_),
transform_(pp.transform_),
rotationAxis_(pp.rotationAxis_),
rotationCentre_(pp.rotationCentre_)
rotationCentre_(pp.rotationCentre_),
separationVector_(pp.separationVector_)
{
calcTransforms();
}
@ -1322,6 +1391,8 @@ void Foam::cyclicPolyPatch::write(Ostream& os) const
{
os.writeKeyword("transform") << transformTypeNames[TRANSLATIONAL]
<< token::END_STATEMENT << nl;
os.writeKeyword("separationVector") << separationVector_
<< token::END_STATEMENT << nl;
break;
}
default:

View File

@ -101,12 +101,19 @@ private:
//- Type of transformation - rotational or translational
transformType transform_;
// For rotation
//- Axis of rotation for rotational cyclics
vector rotationAxis_;
//- point on axis of rotation for rotational cyclics
point rotationCentre_;
// For translation
//- Translation vector
vector separationVector_;
// Private member functions
@ -267,6 +274,9 @@ public:
const edgeList& coupledEdges() const;
// Transformation
vector separation(const label facei) const
{
if (facei < size()/2)
@ -328,6 +338,31 @@ public:
}
}
//- Type of transform
transformType transform() const
{
return transform_;
}
//- Axis of rotation for rotational cyclics
const vector& rotationAxis() const
{
return rotationAxis_;
}
//- point on axis of rotation for rotational cyclics
const point& rotationCentre() const
{
return rotationCentre_;
}
//- Translation vector for translational cyclics
const vector& separationVector() const
{
return separationVector_;
}
//- Initialize ordering for primitivePatch. Does not
// refer to *this (except for name() and type() etc.)