exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

MixingPlane with major mass conservation problem solved.

Browse source
This commit is contained in:
Martin Beaudoin 2011-12-20 15:44:43 -05:00
parent da4aef6f02
commit e905452e43
15 changed files with 692 additions and 169 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

219
src/OpenFOAM/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolate.C
View file

@ -132,8 +132,10 @@ interpolate
}
}
if (debug <= -200)
if (debug <= -500)
{
error::printStack(Info);
Info<< "srcF : " << srcF << nl
<< "srcAddr : " << srcAddr << nl
<< "srcWeights: " << srcWeights << nl
@ -214,10 +216,15 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToSlave
// MB: We need this back
Field<Type> profileFF = transform(masterPatchToProfileT(), patchFF);
Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToSlave: "
<< "patchFF: " << patchFF
<< "profileFF: " << profileFF
<< endl;
if (debug > 1)
{
Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToSlave: "
<< "patchFF: " << patchFF << endl
<< "profileFF: " << profileFF << endl
<< "masterPatchToProfileT(): " << masterPatchToProfileT() << endl
<< "slaveProfileToPatchT(): " << slaveProfileToPatchT() << endl
<< endl;
}
// Do interpolation
tmp<Field<Type> > tresult
@ -228,18 +235,19 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToSlave
pTraits<Type>::zero
)
);
Field<Type>& result = tresult();
interpolate
(
profileFF, // Master data in 'profile space'
masterPatchToProfileAddr(), // From master: compute the average
masterPatchToProfileWeights(),
slaveProfileToPatchAddr(), // To slave we distribute the average from
slaveProfileToPatchWeights(), // profile to patch
result
);
(
profileFF, // Master data in 'profile space'
masterPatchToProfileAddr(), // From master: compute the average
masterPatchToProfileWeights(),
slaveProfileToPatchAddr(), // To slave we distribute the average from
slaveProfileToPatchWeights(), // profile to patch
result
);
// Apply transform to bring the slave field back from 'profile space'
// to 'patch space'
transform(result, slaveProfileToPatchT(), result); // MB: We need this back
@ -284,16 +292,14 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToMaster
// Move slave data from 'patch space' to 'profile space'
Field<Type> profileFF = transform(slavePatchToProfileT(), patchFF);
Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToMaster: "
<< "patchFF: " << patchFF
<< "profileFF: " << profileFF
<< endl;
if (sizeof(Type)/sizeof(scalar) == 3)
if (debug > 1)
{
Field<Type> dummypatchFF(patchFF);
dummypatchFF.replace(vector::Z, scalar(0.0));
Info << "dummypatchFF: " << mag(dummypatchFF) << endl;
Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToMaster: "
<< "patchFF: " << patchFF << endl
<< "profileFF: " << profileFF << endl
<< "slavePatchToProfileT(): " << slavePatchToProfileT() << endl
<< "masterProfileToPatchT(): " << masterProfileToPatchT() << endl
<< endl;
}
// Do interpolation
@ -305,18 +311,19 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToMaster
pTraits<Type>::zero
)
);
Field<Type>& result = tresult();
interpolate
(
profileFF, // Slave data in 'profile space'
slavePatchToProfileAddr(), // From slave: compute the average
slavePatchToProfileWeights(),
masterProfileToPatchAddr(), // To master: distribute the average
masterProfileToPatchWeights(),
result
);
(
profileFF, // Slave data in 'profile space'
slavePatchToProfileAddr(), // From slave: compute the average
slavePatchToProfileWeights(),
masterProfileToPatchAddr(), // To master: distribute the average
masterProfileToPatchWeights(),
result
);
// Apply transform to bring the master field back from 'profile space'
// to 'patch space'
transform(result, masterProfileToPatchT(), result);
@ -338,6 +345,152 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToMaster
return tint;
}
template<class MasterPatch, class SlavePatch>
template<class Type>
tmp<Field<Type> >
MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToMaster
(
const Field<Type>& patchFF
) const
{
if (patchFF.size() != masterPatch_.size())
{
FatalErrorIn
(
"MixingPlaneInterpolation::masterToMaster("
"const Field<Type> ff)"
) << "given field does not correspond to patch. Patch size: "
<< masterPatch_.size() << " field size: " << patchFF.size()
<< abort(FatalError);
}
// Move master data from 'patch space' to 'profile space'
Field<Type> profileFF = transform(masterPatchToProfileT(), patchFF);
if (debug > 1)
{
Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToMaster: "
<< "patchFF: " << patchFF << endl
<< "profileFF: " << profileFF << endl
<< "masterPatchToProfileT(): " << masterPatchToProfileT() << endl
<< endl;
}
// Do interpolation
tmp<Field<Type> > tresult
(
new Field<Type>
(
masterPatch_.size(),
pTraits<Type>::zero
)
);
Field<Type>& result = tresult();
interpolate
(
profileFF, // Master data in 'profile space'
masterPatchToProfileAddr(), // From master: compute the average
masterPatchToProfileWeights(),
masterProfileToPatchAddr(), // To master: distribute the average
masterProfileToPatchWeights(),
result
);
// Apply transform to bring the master field back from 'profile space'
// to 'patch space'
transform(result, masterProfileToPatchT(), result);
return tresult;
}
template<class MasterPatch, class SlavePatch>
template<class Type>
tmp<Field<Type> >
MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToMaster
(
const tmp<Field<Type> >& tff
) const
{
tmp<Field<Type> > tint = masterToMaster(tff());
tff.clear();
return tint;
}
template<class MasterPatch, class SlavePatch>
template<class Type>
tmp<Field<Type> >
MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToSlave
(
const Field<Type>& patchFF
) const
{
if (patchFF.size() != slavePatch_.size())
{
FatalErrorIn
(
"MixingPlaneInterpolation::slaveToSlave("
"const Field<Type> ff)"
) << "given field does not correspond to patch. Patch size: "
<< slavePatch_.size() << " field size: " << patchFF.size()
<< abort(FatalError);
}
// Move slave data from 'patch space' to 'profile space'
Field<Type> profileFF = transform(slavePatchToProfileT(), patchFF);
if (debug > 1)
{
Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToSlave: "
<< "patchFF: " << patchFF << endl
<< "profileFF: " << profileFF << endl
<< "slavePatchToProfileT(): " << slavePatchToProfileT() << endl
<< endl;
}
// Do interpolation
tmp<Field<Type> > tresult
(
new Field<Type>
(
slavePatch_.size(),
pTraits<Type>::zero
)
);
Field<Type>& result = tresult();
interpolate
(
profileFF, // Slave data in 'profile space'
slavePatchToProfileAddr(), // From slave: compute the average
slavePatchToProfileWeights(),
slaveProfileToPatchAddr(), // To slave: distribute the average
slaveProfileToPatchWeights(),
result
);
// Apply transform to bring the slave field back from 'profile space'
// to 'patch space'
transform(result, slaveProfileToPatchT(), result);
return tresult;
}
template<class MasterPatch, class SlavePatch>
template<class Type>
tmp<Field<Type> >
MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToSlave
(
const tmp<Field<Type> >& tff
) const
{
tmp<Field<Type> > tint = slaveToSlave(tff());
tff.clear();
return tint;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

30
src/OpenFOAM/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolation.H
View file

@ -476,6 +476,36 @@ public:
) const;
//- Interpolate from master to master
//- (self circumferential averaging)
template<class Type>
tmp<Field<Type> > masterToMaster
(
const Field<Type>& pf
) const;
template<class Type>
tmp<Field<Type> > masterToMaster
(
const tmp<Field<Type> >& tpf
) const;
//- Interpolate from slave to slave
//- (self circumferential averaging)
template<class Type>
tmp<Field<Type> > slaveToSlave
(
const Field<Type>& pf
) const;
template<class Type>
tmp<Field<Type> > slaveToSlave
(
const tmp<Field<Type> >& tpf
) const;
// Edit
//- Correct weighting factors for moving mesh.

368
src/OpenFOAM/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationAddressing.C
View file

@ -39,6 +39,7 @@ Contributor
#include "PrimitivePatch.H"
#include "IOmanip.H"
#include "transform.H"
#include "RodriguesRotation.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -203,78 +204,154 @@ template<class MasterPatch, class SlavePatch>
void
MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcTransforms() const
{
// Collapse in spawise direction
// Collapse in spanwise direction
const direction spanDir = spanwiseSwitch();
// Master side
#define useRodrigues
#ifndef useRodrigues
const direction directionDir = directionalSwitch();
#endif
// Master side
masterPatchToProfileTPtr_ = new tensorField(masterPatch_.size());
tensorField& mPatchToProfileT = *masterPatchToProfileTPtr_;
masterProfileToPatchTPtr_ = new tensorField(masterPatch_.size());
tensorField& mProfileToPatchT = *masterProfileToPatchTPtr_;
// Get master patch face normals
const vectorField& globalMasterNormals = masterPatch_.faceNormals();
// Move face normals into the local coordinate system
vectorField localMasterNormals = cs_.localVector(globalMasterNormals);
localMasterNormals.replace(spanDir, 0);
// Transform back to global
vectorField transformMasterNormals = cs_.globalVector(localMasterNormals);
// Calculate transform tensors
mPatchToProfileT =
rotationTensor
(
globalMasterNormals,
transformMasterNormals
);
mProfileToPatchT =
rotationTensor
(
transformMasterNormals,
globalMasterNormals
);
if (debug != 0)
if (cs_.name() == "coordinateSystem")
{
// We can remove this later when agree this is ok.
const vectorField& masterFaceCntr = masterPatch_.faceCentres();
pointField transformedFaceCentre1 = transform
(
mPatchToProfileT,
masterFaceCntr
);
pointField transformedFaceCentre2 = transform
(
mProfileToPatchT,
transformedFaceCentre1
);
Info << "::calcTransforms(): globalMasterNormals: " << globalMasterNormals << endl;
Info << "::calcTransforms(): localMasterNormals: " << localMasterNormals << endl;
Info << "::calcTransforms(): transformMasterNormals: " << transformMasterNormals << endl;
Info << "::calcTransforms(): masterFaceCntr: " << masterFaceCntr << endl;
Info << "::calcTransforms(): transformedFaceCentre1: " << transformedFaceCentre1 << endl;
Info << "::calcTransforms(): transformedFaceCentre2: " << transformedFaceCentre2 << endl;
InfoIn
(
"MixingPlaneInterpolation"
"<MasterPatch, SlavePatch>::calcTransforms()"
) << "master face centre transformation check: "
<< "(should be zero!) = "
<< " mag : " << mag(transformedFaceCentre2 - masterFaceCntr) << endl
<< " sum mag: " << sum(mag(transformedFaceCentre2 - masterFaceCntr))
<< endl;
// Identity tensor for cartesian space
mPatchToProfileT = tensor(sphericalTensor::I);
mProfileToPatchT = tensor(sphericalTensor::I);
}
else
{
// Get master patch face centers.
// Warning: Face normals are not a good choice here. We need a vector
// that will intersect the rotation axis. Furthermore, if the surface
// normals are all parallel with the rotation axis, no valid rotation
// tensors can be computed
// Use face cell centers cause the fields are taken from cell centers
vectorField globalMasterVectors =
masterPatch_.faceCellCentres() - cs_.origin();
// We need unit vectors for computing rotation tensor
// We also need vectors lying into the plane normal to the rotation
// axis. This is a major limitation of the current implementation of
// rotationTensor() that is being used for computing rotation tensors:
// we cannot specify the rotation axis. Instead, when using
// rotationTensor(), the rotation axis will be aligned with the normal
// of the plane spanned by the two vectors.
// RodriguesRotation() was designed to overcome these limitations.
#ifndef useRodrigues
globalMasterVectors.replace(directionDir, 0);
Info << "Debug: " << masterPatch_.faceCellCentres() - globalMasterVectors << endl;
globalMasterVectors /= mag(globalMasterVectors);
#endif
// Move face vector into the local coordinate system
vectorField localMasterVectors = cs_.localVector(globalMasterVectors);
// Translate everything to theta=0
localMasterVectors.replace(spanDir, 0);
// Transform back to global
vectorField transformMasterVectors = cs_.globalVector(localMasterVectors);
if( debug )
{
Info << "Debug: globalMasterVectors: " << globalMasterVectors << endl;
Info << "Debug: transformMasterVectors: " << transformMasterVectors << endl;
}
// Calculate transform tensors. These are pure rotation tensors, aligned
// with the reference frame rotation axis
mPatchToProfileT =
#ifndef useRodrigues
rotationTensor
(
#else
RodriguesRotation
(
cs_.axis(),
#endif
globalMasterVectors,
transformMasterVectors
);
mProfileToPatchT =
#ifndef useRodrigues
rotationTensor
(
#else
RodriguesRotation
(
cs_.axis(),
#endif
transformMasterVectors,
globalMasterVectors
);
#if 0
{
vectorField globalMasterVectors2 = masterPatch_.faceCellCentres() - cs_.origin();
globalMasterVectors2 /= mag(globalMasterVectors);
vectorField localMasterVectors2 = cs_.localVector(globalMasterVectors2);
localMasterVectors2.replace(spanDir, 0);
vectorField transformMasterVectors2 = cs_.globalVector(localMasterVectors2);
Info << "master delta face transf: " << mag(mPatchToProfileT - RodriguesRotation(cs_.axis(), globalMasterVectors, transformMasterVectors)) << endl;
Info << "master delta face transf: " << mag(mProfileToPatchT - RodriguesRotation(cs_.axis(), transformMasterVectors, globalMasterVectors )) << endl;
Info << "master delta face: " << mag(mPatchToProfileT - RodriguesRotation(cs_.axis(), globalMasterVectors2, transformMasterVectors2)) << endl;
Info << "master delta face: " << mag(mProfileToPatchT - RodriguesRotation(cs_.axis(), transformMasterVectors2, globalMasterVectors2 )) << endl;
}
#endif
if (debug != 0)
{
// We can remove this later when agree this is ok.
const vectorField& masterFaceCntr = masterPatch_.faceCellCentres();
pointField transformedFaceCentre1 = transform
(
mPatchToProfileT,
masterFaceCntr
);
pointField transformedFaceCentre2 = transform
(
mProfileToPatchT,
transformedFaceCentre1
);
Info << "::calcTransforms(): mPatchToProfileT: " << mPatchToProfileT << endl;
Info << "::calcTransforms(): mProfileToPatchT: " << mProfileToPatchT << endl;
Info << "::calcTransforms(): globalMasterVectors: " << globalMasterVectors << endl;
Info << "::calcTransforms(): localMasterVectors: " << localMasterVectors << endl;
Info << "::calcTransforms(): transformMasterVectors: " << transformMasterVectors << endl;
Info << "::calcTransforms(): masterFaceCntr: " << masterFaceCntr << endl;
Info << "::calcTransforms(): transformedFaceCentre1: " << transformedFaceCentre1 << endl;
Info << "::calcTransforms(): transformedFaceCentre2: " << transformedFaceCentre2 << endl;
Info << "::calcTransforms(): deltaMag direct: " << mag(transformedFaceCentre1) - mag(masterFaceCntr) << endl;
Info << "::calcTransforms(): deltaMag reverse: " << mag(transformedFaceCentre2) - mag(masterFaceCntr) << endl;
InfoIn
(
"MixingPlaneInterpolation"
"<MasterPatch, SlavePatch>::calcTransforms()"
) << "master face centre transformation check: "
<< "(should be zero!) = "
<< " mag : " << mag(transformedFaceCentre2 - masterFaceCntr) << endl
<< " sum mag: " << sum(mag(transformedFaceCentre2 - masterFaceCntr))
<< endl;
}
}
// Slave side
@ -284,65 +361,128 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcTransforms() const
slaveProfileToPatchTPtr_ = new tensorField(slavePatch_.size());
tensorField& sProfileToPatchT = *slaveProfileToPatchTPtr_;
// Get slave patch face normals
const vectorField& globalSlaveNormals = slavePatch_.faceNormals();
// Move face normals into the local coordinate system
vectorField localSlaveNormals = cs_.localVector(globalSlaveNormals);
localSlaveNormals.replace(spanDir, 0);
// Transform back to global
vectorField transformSlaveNormals = cs_.globalVector(localSlaveNormals);
// Calculate transform tensors
sPatchToProfileT =
rotationTensor
(
globalSlaveNormals,
transformSlaveNormals
);
sProfileToPatchT =
rotationTensor
(
transformSlaveNormals,
globalSlaveNormals
);
if (debug != 0)
if (cs_.name() == "coordinateSystem")
{
// We can remove this later when agree this is ok.
const vectorField& slaveFaceCntr = slavePatch_.faceCentres();
// Identity tensor for cartesian space
sPatchToProfileT = tensor(sphericalTensor::I);
sProfileToPatchT = tensor(sphericalTensor::I);
}
else
{
// Get slave patch face centers.
// Warning: Face normals are not a good choice here. We need a vector
// that will intersect the rotation axis. Furthermore, if the surface
// normals are all parallel with the rotation axis, no valid rotation
// tensors can be computed
vectorField globalSlaveVectors =
slavePatch_.faceCellCentres() - cs_.origin();
pointField transformedFaceCentre1 = transform
(
sPatchToProfileT,
slaveFaceCntr
);
// We need unit vectors for computing rotation tensor
// We also need vectors lying into the plane normal to the rotation
// axis. This is a major limitation of the current implementation of
// rotationTensor() that is being used for computing rotation tensors:
// we cannot specify the rotation axis. Instead, when using
// rotationTensor(), the rotation axis will be aligned with the normal
// of the plane spanned by the two vectors.
// RodriguesRotation() was designed to overcome these limitations.
#ifndef useRodrigues
globalSlaveVectors.replace(directionDir, 0);
pointField transformedFaceCentre2 = transform
(
sProfileToPatchT,
transformedFaceCentre1
);
Info << "Debug: " << slavePatch_.faceCellCentres() - globalSlaveVectors << endl;
Info << "::calcTransforms(): globalSlaveNormals: " << globalSlaveNormals << endl;
Info << "::calcTransforms(): localSlaveNormals: " << localSlaveNormals << endl;
Info << "::calcTransforms(): transformSlaveNormals: " << transformSlaveNormals << endl;
Info << "::calcTransforms(): slaveFaceCntr: " << slaveFaceCntr << endl;
Info << "::calcTransforms(): transformedFaceCentre1: " << transformedFaceCentre1 << endl;
Info << "::calcTransforms(): transformedFaceCentre2: " << transformedFaceCentre2 << endl;
globalSlaveVectors /= mag(globalSlaveVectors);
#endif
InfoIn
(
"MixingPlaneInterpolation"
"<SlavePatch, SlavePatch>::calcTransforms()"
) << "slave face centre transformation check: "
<< "(should be zero!) = "
<< " mag : " << mag(transformedFaceCentre2 - slaveFaceCntr) << endl
<< " sum mag: " << sum(mag(transformedFaceCentre2 - slaveFaceCntr))
<< endl;
// Move face vector into the local coordinate system
vectorField localSlaveVectors = cs_.localVector(globalSlaveVectors);
// Translate everything to theta=0
localSlaveVectors.replace( spanDir, 0 );
// Transform back to global
vectorField transformSlaveVectors = cs_.globalVector(localSlaveVectors);
// Calculate transform tensors. These are pure rotation tensors, aligned
// with the reference frame rotation axis
sPatchToProfileT =
#ifndef useRodrigues
rotationTensor
(
#else
RodriguesRotation
(
cs_.axis(),
#endif
globalSlaveVectors,
transformSlaveVectors
);
sProfileToPatchT =
#ifndef useRodrigues
rotationTensor
(
#else
RodriguesRotation
(
cs_.axis(),
#endif
transformSlaveVectors,
globalSlaveVectors
);
#if 0
{
vectorField globalSlaveVectors2 = slavePatch_.faceCellCentres() - cs_.origin();
globalSlaveVectors2 /= mag(globalSlaveVectors);
vectorField localSlaveVectors2 = cs_.localVector(globalSlaveVectors2);
localSlaveVectors2.replace(spanDir, 0);
vectorField transformSlaveVectors2 = cs_.globalVector(localSlaveVectors2);
Info << "slave delta face transf: " << mag(sPatchToProfileT - RodriguesRotation(cs_.axis(), globalSlaveVectors, transformSlaveVectors)) << endl;
Info << "slave delta face transf: " << mag(sProfileToPatchT - RodriguesRotation(cs_.axis(), transformSlaveVectors, globalSlaveVectors )) << endl;
Info << "slave delta face: " << mag(sPatchToProfileT - RodriguesRotation(cs_.axis(), globalSlaveVectors2, transformSlaveVectors2)) << endl;
Info << "slave delta face: " << mag(sProfileToPatchT - RodriguesRotation(cs_.axis(), transformSlaveVectors2, globalSlaveVectors2 )) << endl;
}
#endif
if (debug != 0)
{
// We can remove this later when agree this is ok.
const vectorField& slaveFaceCntr = slavePatch_.faceCellCentres();
pointField transformedFaceCentre1 = transform
(
sPatchToProfileT,
slaveFaceCntr
);
pointField transformedFaceCentre2 = transform
(
sProfileToPatchT,
transformedFaceCentre1
);
Info << "::calcTransforms(): globalSlaveVectors: " << globalSlaveVectors << endl;
Info << "::calcTransforms(): localSlaveVectors: " << localSlaveVectors << endl;
Info << "::calcTransforms(): transformSlaveVectors: " << transformSlaveVectors << endl;
Info << "::calcTransforms(): slaveFaceCntr: " << slaveFaceCntr << endl;
Info << "::calcTransforms(): transformedFaceCentre1: " << transformedFaceCentre1 << endl;
Info << "::calcTransforms(): transformedFaceCentre2: " << transformedFaceCentre2 << endl;
Info << "::calcTransforms(): deltaMag direct: " << mag(transformedFaceCentre1) - mag(slaveFaceCntr) << endl;
Info << "::calcTransforms(): deltaMag reverse: " << mag(transformedFaceCentre2) - mag(slaveFaceCntr) << endl;
InfoIn
(
"MixingPlaneInterpolation"
"<SlavePatch, SlavePatch>::calcTransforms()"
) << "slave face centre transformation check: "
<< "(should be zero!) = "
<< " mag : " << mag(transformedFaceCentre2 - slaveFaceCntr) << endl
<< " sum mag: " << sum(mag(transformedFaceCentre2 - slaveFaceCntr))
<< endl;
}
}
}

16
src/OpenFOAM/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationPatches.C
View file

@ -243,6 +243,17 @@ calcTransformedPatches() const
correctStraddlingFaces(masterFaces, masterPointsLocalCoord);
correctStraddlingFaces(slaveFaces, slavePointsLocalCoord);
if(debug)
{
Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcTransformedPatches(): "
<< "masterPointsLocalCoord: "
<< masterPointsLocalCoord << endl;
Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcTransformedPatches(): "
<< "slavePointsLocalCoord: "
<< slavePointsLocalCoord << endl;
}
// Create the local coords patches
transformedMasterPatchPtr_ =
@ -372,6 +383,8 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcMixingPlanePatch() const
mixingPatchFaces[fI] = curFace;
#if 0
// Work in progress... MB 06_2011
// Compute face area
if (spanLimited[spanDir].first() || spanLimited[spanDir].first())
{
@ -381,6 +394,7 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcMixingPlanePatch() const
}
else
mixingPatchFacesArea[fI] = curFace.mag(mixingPatchPoints);
#endif
}
mixingPlanePatchPtr_ =
@ -390,10 +404,10 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcMixingPlanePatch() const
mixingPatchPoints
);
Info << "::calcMixingPlanePatch():mixingPatchFacesArea(): " << mixingPatchFacesArea << endl;
if (debug > 0)
{
Info << "::calcMixingPlanePatch():mixingPatchFacesArea(): " << mixingPatchFacesArea << endl;
InfoIn
(
"void MixingPlaneInterpolation::calcMixingPlanePatch()"

13
src/OpenFOAM/interpolations/MixingPlaneInterpolation/mixingPlaneInterpolation.H
View file

@ -44,22 +44,35 @@ Description
#include "face.H"
#include "SubList.H"
#include "pointField.H"
#include "polyPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
#if 0
typedef MixingPlaneInterpolation
<
PrimitivePatch<face, SubList, const pointField&>,
PrimitivePatch<face, SubList, const pointField&>
> mixingPlaneInterpolation;
#else
typedef MixingPlaneInterpolation
<
polyPatch,
polyPatch
> mixingPlaneInterpolation;
#endif
#if 0
// Activated later, when parallelization
typedef MixingPlaneInterpolation
<
PrimitivePatch<face, List, const pointField&>,
PrimitivePatch<face, List, const pointField&>
> mixingPlaneZoneInterpolation;
#endif
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

4
src/OpenFOAM/matrices/lduMatrix/preconditioners/DILUPreconditioner/DILUPreconditioner.C
View file

@ -32,6 +32,10 @@ namespace Foam
{
defineTypeNameAndDebug(DILUPreconditioner, 0);
lduPreconditioner::
addsymMatrixConstructorToTable<DILUPreconditioner>
addDILUPreconditionerSymMatrixConstructorToTable_;
lduPreconditioner::
addasymMatrixConstructorToTable<DILUPreconditioner>
addDILUPreconditionerAsymMatrixConstructorToTable_;

19
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/mixingPlane/mixingPlanePolyPatch.C
View file

@ -407,10 +407,14 @@ Foam::label Foam::mixingPlanePolyPatch::shadowIndex() const
if (!shadow.active())
{
FatalErrorIn("label mixingPlanePolyPatch::shadowIndex() const")
WarningIn("label mixingPlanePolyPatch::shadowIndex() const")
<< "Shadow patch name " << shadowName_
<< " not found. Please check your MixingPlane definition."
<< abort(FatalError);
<< " not found. Please check your MixingPlane definition. "
<< "This may be fine at mesh generation stage."
<< endl;
// Return a large label to indicate "undefined" or slave side
return 99999;
}
shadowIndex_ = shadow.index();
@ -507,7 +511,14 @@ void Foam::mixingPlanePolyPatch::calcGeometry()
// Reconstruct the cell face centres
if (patchToPatchPtr_ && master())
{
// Compute the neighbour face cell center
reconFaceCellCentres();
// Next, identify which cells are located at these locations
// Next, compute the weighting factors in order to properly interpolate
// the field values at those locations. We will be using an inverse
// distance interpolation scheme.
}
calcTransforms();
@ -579,7 +590,7 @@ void Foam::mixingPlanePolyPatch::write(Ostream& os) const
<< token::END_STATEMENT << nl;
// Note: only master writes the data
if (master())
if (master() || shadowIndex_ == -1)
{
// Write coordinate system dictionary. Check by hand. HJ, 26/Jan/2011
os.writeKeyword("coordinateSystem");

9
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/mixingPlane/mixingPlanePolyPatch.H
View file

@ -268,6 +268,15 @@ public:
template<class Type>
tmp<Field<Type> > interpolate(const tmp<Field<Type> >& tpf) const;
//- Interpolate face field from average: given field on a the
//- master/slave side, create an interpolated field on the same side
//- using averaged values
template<class Type>
tmp<Field<Type> > circumferentialAverage(const Field<Type>& pf) const;
template<class Type>
tmp<Field<Type> > circumferentialAverage(const tmp<Field<Type> >& tpf) const;
//- Return reconstructed cell centres
const vectorField& reconFaceCellCentres() const;

28
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/mixingPlane/mixingPlanePolyPatchTemplates.C
View file

@ -62,5 +62,33 @@ Foam::tmp<Foam::Field<Type> > Foam::mixingPlanePolyPatch::interpolate
return tint;
}
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::mixingPlanePolyPatch::circumferentialAverage
(
const Field<Type>& pf
) const
{
if (master())
{
return patchToPatch().masterToMaster(pf);
}
else
{
return patchToPatch().slaveToSlave(pf);
}
}
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::mixingPlanePolyPatch::circumferentialAverage
(
const tmp<Field<Type> >& tpf
) const
{
tmp<Field<Type> > tint = circumferentialAverage(tpf());
tpf.clear();
return tint;
}
// ************************************************************************* //

38
src/finiteVolume/fields/fvPatchFields/constraint/mixingPlane/mixingPlaneFvPatchField.C
View file

@ -31,6 +31,11 @@ Contributor
\*---------------------------------------------------------------------------*/
#include "mixingPlaneFvPatchField.H"
#include "volMesh.H"
#include "surfaceMesh.H"
#include "fvsPatchField.H"
#include "volFields.H"
#include "interpolationCellPoint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -138,13 +143,13 @@ mixingPlaneFvPatchField<Type>::mixingPlaneFvPatchField
template<class Type>
tmp<Field<Type> > mixingPlaneFvPatchField<Type>::patchNeighbourField() const
{
//if(debug == -99)
if(debug > 1)
{
Info << "mixingPlaneFvPatchField<Type>::patchNeighbourField(): for field: " << this->dimensionedInternalField().name();
Info << "mixingPlaneFvPatchField<Type>::patchNeighbourField(): for field: " << this->dimensionedInternalField().name();
Info << " on patch: " << this->patch().name() << endl;
Info << " surface Area: " << gSum(this->patch().magSf()) << endl;
}
const Field<Type>& iField = this->internalField();
// Get shadow face-cells and assemble shadow field
@ -175,6 +180,9 @@ void mixingPlaneFvPatchField<Type>::initEvaluate
const Pstream::commsTypes commsType
)
{
if(debug)
Info << "Inside mixingPlaneFvPatchField<Type>::initEvaluate: for field: " << this->dimensionedInternalField().name() << endl;
if(this->size() > 0)
{
if (!this->updated())
@ -257,6 +265,30 @@ void mixingPlaneFvPatchField<Type>::updateInterfaceMatrix
) const
{}
// Return averaged field on patch
template<class Type>
tmp<Field<Type> > mixingPlaneFvPatchField<Type>::patchCircumferentialAverageField() const
{
// Compute circum average of self
return this->mixingPlanePatch_.circumferentialAverage(*this);
}
// Return the averaged field for patchInternalField
template<class Type>
tmp<Field<Type> > mixingPlaneFvPatchField<Type>::patchInternalField() const
{
tmp<Field<Type> > tpnf(new Field<Type>(this->size()));
Field<Type>& pnf = tpnf();
if (this->size() > 0)
{
pnf = this->mixingPlanePatch_.circumferentialAverage(fvPatchField<Type>::patchInternalField());
}
return tpnf;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

15
src/finiteVolume/fields/fvPatchFields/constraint/mixingPlane/mixingPlaneFvPatchField.H
View file

@ -132,6 +132,9 @@ public:
// Evaluation functions
//- Return my average field given internal cell data
virtual tmp<Field<Type> > patchInternalField() const;
//- Return neighbour field given internal cell data
virtual tmp<Field<Type> > patchNeighbourField() const;
@ -166,10 +169,20 @@ public:
const Pstream::commsTypes commsType
) const;
//- Return patch averaged field given internal cell data
virtual tmp<Field<Type> > patchCircumferentialAverageField() const;
//- Return neighbour averaged field given internal cell data
virtual tmp<Field<Type> > patchNeighbourAveragedField() const
{
// Already averaged from patchNeighbourField
return patchNeighbourField();
}
// MixingPlane coupled interface functions
//- Does the patch field perform the transfromation
//- Does the patch field perform the transformation
virtual bool doTransform() const
{
return

20
src/finiteVolume/fvMesh/fvPatches/constraint/mixingPlane/mixingPlaneFvPatch.C
View file

@ -79,6 +79,9 @@ void Foam::mixingPlaneFvPatch::makeWeights(scalarField& w) const
w = interpolate(1 - masterWeights);
}
}
if (debug)
Info << ":mixingPlaneFvPatch::makeWeights: w: " << w << endl;
}
@ -98,6 +101,9 @@ void Foam::mixingPlaneFvPatch::makeDeltaCoeffs(scalarField& dc) const
dc = interpolate(masterDeltas);
}
}
if (debug)
Info << ":mixingPlaneFvPatch::makeDeltaCoeffs: dc: " << dc << endl;
}
@ -105,6 +111,8 @@ void Foam::mixingPlaneFvPatch::makeCorrVecs(vectorField& cv) const
{
if(size() > 0)
{
cv = vector::zero;
#if 0
// Full non-orthogonality treatment
// Calculate correction vectors on coupled patches
@ -113,7 +121,11 @@ void Foam::mixingPlaneFvPatch::makeCorrVecs(vectorField& cv) const
vectorField patchDeltas = delta();
vectorField n = nf();
cv = n - patchDeltas*patchDeltaCoeffs;
#endif
}
//if (debug)
Info << ":mixingPlaneFvPatch::makeCorrVecs: cv: " << cv << endl;
}
@ -134,6 +146,9 @@ Foam::tmp<Foam::vectorField> Foam::mixingPlaneFvPatch::delta() const
tmp<vectorField> tDelta =
mixingPlanePolyPatch_.reconFaceCellCentres() - Cn();
if(debug)
Info << "mixingPlaneFvPatch::delta: tDelta: " << tDelta() << endl;
return tDelta;
}
else
@ -144,6 +159,9 @@ Foam::tmp<Foam::vectorField> Foam::mixingPlaneFvPatch::delta() const
- mixingPlanePolyPatch_.shadow().reconFaceCellCentres()
);
if(debug)
Info << "mixingPlaneFvPatch::delta: tDelta: " << tDelta() << endl;
return tDelta;
}
}
@ -186,6 +204,8 @@ const Foam::labelListList& Foam::mixingPlaneFvPatch::addressing() const
const Foam::scalarListList& Foam::mixingPlaneFvPatch::weights() const
{
Info << "Foam::scalarListList& Foam::mixingPlaneFvPatch::weights()" << endl;
if (mixingPlanePolyPatch_.master())
{
return mixingPlanePolyPatch_.patchToPatch().

12
src/finiteVolume/fvMesh/fvPatches/constraint/mixingPlane/mixingPlaneFvPatch.H
View file

@ -127,6 +127,18 @@ public:
return mixingPlanePolyPatch_.interpolate(tpf);
}
//- Circumferential average face field
template<class Type>
tmp<Field<Type> > circumferentialAverage(const Field<Type>& pf) const
{
return mixingPlanePolyPatch_.circumferentialAverage(pf);
}
template<class Type>
tmp<Field<Type> > circumferentialAverage(const tmp<Field<Type> >& tpf) const
{
return mixingPlanePolyPatch_.circumferentialAverage(tpf);
}
// Interface transfer functions

17
src/finiteVolume/fvMesh/fvPatches/constraint/mixingPlane/traceMixingPlaneFlux.H
View file

@ -40,6 +40,7 @@ namespace Foam
// Temporary hack: useful for tracking balance of phi across interface
void
traceMixingPlaneFlux(
volVectorField& U,
surfaceScalarField& phi,
scalar& masterPatchScaleFactor,
scalar& shadowPatchScaleFactor,
@ -78,12 +79,22 @@ traceMixingPlaneFlux(
//scalar shadowFluxMag = shadowPatchScaleFactor_ * gSumMag(phi.boundaryField()[shadowPatchI]);
scalar shadowFluxMag = mag(shadowFlux);
scalar phiFromU(mag(gSum(U.boundaryField()[patchI] & U.mesh().Sf().boundaryField()[patchI])));
scalar shadowPhiFromU(mag(gSum(U.boundaryField()[shadowPatchI] & U.mesh().Sf().boundaryField()[shadowPatchI])));
Info<< "====> traceMixingPlaneFlux::: ID: " << traceId << ": mixingPlane pair (" << mixingPlanePatch.name() << ", " << mixingPlanePatch.shadow().name() << ") : "
<< " flux: " << localFlux << " " << shadowFlux
<< " : mag: " << localFluxMag << " " << shadowFluxMag
<< " Diff = " << localFlux + shadowFlux << " or "
<< " Diff mag = " << localFlux + shadowFlux << " or "
<< mag(localFlux + shadowFlux)/(localFluxMag + SMALL)*100
<< " %" << endl;
<< " %"
<< " phiFromU: " << phiFromU
<< " shadowPhiFromU: " << shadowPhiFromU
<< " diff: " << phiFromU - shadowPhiFromU << " or "
<< mag(phiFromU - shadowPhiFromU)/(phiFromU + SMALL)*100
<< " %"
<< endl;
}
}
}

53
src/postProcessing/functionObjects/check/mixingPlaneCheck/mixingPlaneCheckFunctionObject.C
View file

@ -84,14 +84,14 @@ bool Foam::mixingPlaneCheckFunctionObject::execute()
const polyMesh& mesh =
time_.lookupObject<polyMesh>(regionName_);
// const surfaceScalarField& phi =
// mesh.lookupObject<surfaceScalarField>(phiName_);
boolList visited(mesh.boundaryMesh().size(), false);
forAll (mesh.boundaryMesh(), patchI)
{
if (isA<mixingPlanePolyPatch>(mesh.boundaryMesh()[patchI]))
if (
isA<mixingPlanePolyPatch>(mesh.boundaryMesh()[patchI]) &&
mesh.boundaryMesh()[patchI].size()
)
{
if (!visited[patchI])
{
@ -128,8 +128,11 @@ bool Foam::mixingPlaneCheckFunctionObject::execute()
const scalarField masterAreas = mag(mixingMaster.faceAreas());
const scalarField shadowAreas = mag(mixingShadow.faceAreas());
scalar sumMasterAreas = gSum(masterAreas);
scalar sumShadowAreas = gSum(shadowAreas);
// Until the mixingPlane is fully parallelized, we stick with
// the serial version of sum. The interface is residing on a
// single processor when running in parallel
scalar sumMasterAreas = sum(masterAreas);
scalar sumShadowAreas = sum(shadowAreas);
scalar sumMixingAreas = sum(mixingPlanePatchAreas);
@ -165,8 +168,11 @@ bool Foam::mixingPlaneCheckFunctionObject::execute()
}
}
#if 0
// Does not work when in cylindrical coordinates
Info<< "Master scaling = "
<< mixingPlanePatchAreas/masterToStripsAreas << endl;
#endif
// Calculate shadow to strip sum
scalarField shadowToStripsAreas(mixingPlanePatch.size(), 0);
@ -189,8 +195,35 @@ bool Foam::mixingPlaneCheckFunctionObject::execute()
}
}
#if 0
// Does not work when in cylindrical coordinates
Info<< "Shadow scaling = "
<< mixingPlanePatchAreas/shadowToStripsAreas << endl;
#endif
// Old way of computing phi balance
if( mesh.foundObject<surfaceScalarField>(phiName_) )
{
const surfaceScalarField& phi =
mesh.lookupObject<surfaceScalarField>(phiName_);
// Calculate local and shadow flux
scalar masterPatchScaleFactor_ = 1.0;
scalar shadowPatchScaleFactor_ = sumMasterAreas/sumShadowAreas;
scalar localFlux = masterPatchScaleFactor_ * sum(phi.boundaryField()[patchI]);
scalar localFluxMag = mag(localFlux);
scalar shadowFlux = shadowPatchScaleFactor_ * sum(phi.boundaryField()[shadowPatchI]);
scalar shadowFluxMag = mag(shadowFlux);
Info<< "mixingPlane pair " << name() << " (" << mixingMaster.name() << ", " << mixingMaster.shadow().name() << ") : "
<< " flux: " << localFlux << " " << shadowFlux
<< " : mag: " << localFluxMag << " " << shadowFluxMag
<< " Diff = " << localFlux + shadowFlux << " or "
<< mag(localFlux + shadowFlux)/(localFluxMag + SMALL)*100
<< " %" << endl;
}
}
}
}
@ -204,8 +237,8 @@ bool Foam::mixingPlaneCheckFunctionObject::execute()
// visited[patchI] = true;
// // Calculate local and shadow flux
// scalar localFlux = masterPatchScaleFactor_ * gSum(phi.boundaryField()[patchI]);
// //scalar localFluxMag = masterPatchScaleFactor_ * gSumMag(phi.boundaryField()[patchI]);
// scalar localFlux = masterPatchScaleFactor_ * sum(phi.boundaryField()[patchI]);
// //scalar localFluxMag = masterPatchScaleFactor_ * sumMag(phi.boundaryField()[patchI]);
// scalar localFluxMag = mag(localFlux);
// const mixingPlanePolyPatch& mixingPlanePatch =
@ -218,8 +251,8 @@ bool Foam::mixingPlaneCheckFunctionObject::execute()
// visited[shadowPatchI] = true;
// scalar shadowFlux = shadowPatchScaleFactor_ * gSum(phi.boundaryField()[shadowPatchI]);
// //scalar shadowFluxMag = shadowPatchScaleFactor_ * gSumMag(phi.boundaryField()[shadowPatchI]);
// scalar shadowFlux = shadowPatchScaleFactor_ * sum(phi.boundaryField()[shadowPatchI]);
// //scalar shadowFluxMag = shadowPatchScaleFactor_ * sumMag(phi.boundaryField()[shadowPatchI]);
// scalar shadowFluxMag = mag(shadowFlux);
// Info<< "mixingPlane pair " << name_ << " (" << mixingPlanePatch.name() << ", " << mixingPlanePatch.shadow().name() << ") : "