diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C index ece28f7b1..98f2c5ddb 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C +++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C @@ -41,13 +41,13 @@ Description int main(int argc, char *argv[]) { - #include "setRootCase.H" - #include "createTime.H" - #include "createMesh.H" - #include "readGravitationalAcceleration.H" - #include "createFields.H" - #include "createRadiationModel.H" - #include "initContinuityErrs.H" +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" +# include "readGravitationalAcceleration.H" +# include "createFields.H" +# include "createRadiationModel.H" +# include "initContinuityErrs.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -57,17 +57,17 @@ int main(int argc, char *argv[]) { Info<< "Time = " << runTime.timeName() << nl << endl; - #include "readSIMPLEControls.H" - #include "initConvergenceCheck.H" +# include "readSIMPLEControls.H" +# include "initConvergenceCheck.H" p.storePrevIter(); rho.storePrevIter(); // Pressure-velocity SIMPLE corrector { - #include "UEqn.H" - #include "hEqn.H" - #include "pEqn.H" +# include "UEqn.H" +# include "hEqn.H" +# include "pEqn.H" } turbulence->correct(); @@ -78,7 +78,7 @@ int main(int argc, char *argv[]) << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; - #include "convergenceCheck.H" +# include "convergenceCheck.H" } Info<< "End\n" << endl; diff --git a/doc/Doxygen/FoamFooter.html b/doc/Doxygen/FoamFooter.html index 352a32682..3f193433c 100644 --- a/doc/Doxygen/FoamFooter.html +++ b/doc/Doxygen/FoamFooter.html @@ -1,4 +1,6 @@ -Copyright © 2000-2009 OpenCFD Ltd +Copyright © 1993-2000 Henry Weller and Hrvoje Jasak +Copyright © 2000-2006 Programmers manual for FOAM, a product of Nabla Ltd. +Copyright © 2009-2011 OpenFOAM Extend Project diff --git a/doc/Doxygen/FoamHeader.html b/doc/Doxygen/FoamHeader.html index 3314d3dd6..e83fe4801 100644 --- a/doc/Doxygen/FoamHeader.html +++ b/doc/Doxygen/FoamHeader.html @@ -12,8 +12,8 @@ - - + + @@ -45,20 +45,8 @@ horizontal-align: left; "> class=menuLefton >Source Guide - OpenCFD - - - Solutions - - - Contact - - - OpenFOAM + OpenFOAM-Extend Project diff --git a/etc/controlDict b/etc/controlDict index b5a25f98b..9106caabc 100644 --- a/etc/controlDict +++ b/etc/controlDict @@ -904,8 +904,7 @@ Tolerances primitiveMeshFaceFlatnessThreshold 0.8; // Geometric matching tolerances - cyclicMatchTol 1e-4; - processorMatchTol 1e-4; + patchFaceMatchTol 1e-4; // Volumetric closed domain closedDomainTol 1e-10; diff --git a/src/OpenFOAM/containers/Lists/UList/UList.H b/src/OpenFOAM/containers/Lists/UList/UList.H index 184cee2df..5d3c41346 100644 --- a/src/OpenFOAM/containers/Lists/UList/UList.H +++ b/src/OpenFOAM/containers/Lists/UList/UList.H @@ -95,11 +95,13 @@ public: //- Declare friendship with the SubList class friend class SubList; + // Static Member Functions //- Return a null UList inline static const UList& null(); + // Public classes //- Less function class that can be used for sorting diff --git a/src/OpenFOAM/containers/NamedEnum/NamedEnum.C b/src/OpenFOAM/containers/NamedEnum/NamedEnum.C index 552d5593a..69c28ff5c 100644 --- a/src/OpenFOAM/containers/NamedEnum/NamedEnum.C +++ b/src/OpenFOAM/containers/NamedEnum/NamedEnum.C @@ -72,7 +72,7 @@ Enum Foam::NamedEnum::read(Istream& is) const ( "NamedEnum::read(Istream& is) const", is - ) << name << " is not in enumeration " << toc() + ) << name << " is not in enumeration " << toc() << exit(FatalIOError); } diff --git a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C index bfd09b417..f2e7beff4 100644 --- a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C +++ b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C @@ -61,45 +61,45 @@ void inv(Field& tf, const UList& tf1) } scalar scale = magSqr(tf1[0]); - Vector removeCmpts - ( - magSqr(tf1[0].xx())/scale < SMALL, - magSqr(tf1[0].yy())/scale < SMALL, - magSqr(tf1[0].zz())/scale < SMALL - ); - if (removeCmpts.x() || removeCmpts.y() || removeCmpts.z()) + // Fixed terrible hack. HJ, 20/Jan/2011 + boolList removeCmpts(3); + removeCmpts[0] = magSqr(tf1[0].xx())/scale < SMALL; + removeCmpts[1] = magSqr(tf1[0].yy())/scale < SMALL; + removeCmpts[2] = magSqr(tf1[0].zz())/scale < SMALL; + + if (removeCmpts[0] || removeCmpts[1] || removeCmpts[2]) { symmTensorField tf1Plus(tf1); - if (removeCmpts.x()) + if (removeCmpts[0]) { tf1Plus += symmTensor(1,0,0,0,0,0); } - if (removeCmpts.y()) + if (removeCmpts[1]) { tf1Plus += symmTensor(0,0,0,1,0,0); } - if (removeCmpts.z()) + if (removeCmpts[2]) { tf1Plus += symmTensor(0,0,0,0,0,1); } TFOR_ALL_F_OP_FUNC_F(symmTensor, tf, =, inv, symmTensor, tf1Plus) - if (removeCmpts.x()) + if (removeCmpts[0]) { tf -= symmTensor(1,0,0,0,0,0); } - if (removeCmpts.y()) + if (removeCmpts[1]) { tf -= symmTensor(0,0,0,1,0,0); } - if (removeCmpts.z()) + if (removeCmpts[2]) { tf -= symmTensor(0,0,0,0,0,1); } diff --git a/src/OpenFOAM/fields/Fields/tensorField/tensorField.C b/src/OpenFOAM/fields/Fields/tensorField/tensorField.C index 5b3a0c10b..2878746b3 100644 --- a/src/OpenFOAM/fields/Fields/tensorField/tensorField.C +++ b/src/OpenFOAM/fields/Fields/tensorField/tensorField.C @@ -26,6 +26,7 @@ License #include "tensorField.H" #include "transformField.H" +#include "boolList.H" #define TEMPLATE #include "FieldFunctionsM.C" @@ -61,45 +62,46 @@ void inv(Field& tf, const UList& tf1) } scalar scale = magSqr(tf1[0]); - Vector removeCmpts - ( - magSqr(tf1[0].xx())/scale < SMALL, - magSqr(tf1[0].yy())/scale < SMALL, - magSqr(tf1[0].zz())/scale < SMALL - ); - if (removeCmpts.x() || removeCmpts.y() || removeCmpts.z()) + // Fixed terrible hack. HJ, 20/Jan/2011 + boolList removeCmpts(3); + removeCmpts[0] = magSqr(tf1[0].xx())/scale < SMALL; + removeCmpts[1] = magSqr(tf1[0].yy())/scale < SMALL; + removeCmpts[2] = magSqr(tf1[0].zz())/scale < SMALL; + + + if (removeCmpts[0] || removeCmpts[1] || removeCmpts[2]) { tensorField tf1Plus(tf1); - if (removeCmpts.x()) + if (removeCmpts[0]) { tf1Plus += tensor(1,0,0,0,0,0,0,0,0); } - if (removeCmpts.y()) + if (removeCmpts[1]) { tf1Plus += tensor(0,0,0,0,1,0,0,0,0); } - if (removeCmpts.z()) + if (removeCmpts[2]) { tf1Plus += tensor(0,0,0,0,0,0,0,0,1); } TFOR_ALL_F_OP_FUNC_F(tensor, tf, =, inv, tensor, tf1Plus) - if (removeCmpts.x()) + if (removeCmpts[0]) { tf -= tensor(1,0,0,0,0,0,0,0,0); } - if (removeCmpts.y()) + if (removeCmpts[1]) { tf -= tensor(0,0,0,0,1,0,0,0,0); } - if (removeCmpts.z()) + if (removeCmpts[2]) { tf -= tensor(0,0,0,0,0,0,0,0,1); } diff --git a/src/OpenFOAM/fields/PointPatchFields/derived/global/GlobalPointPatchField.C b/src/OpenFOAM/fields/PointPatchFields/derived/global/GlobalPointPatchField.C index 0f4dac775..05dcfccfc 100644 --- a/src/OpenFOAM/fields/PointPatchFields/derived/global/GlobalPointPatchField.C +++ b/src/OpenFOAM/fields/PointPatchFields/derived/global/GlobalPointPatchField.C @@ -30,6 +30,8 @@ License #include "constraints.H" #include "PstreamCombineReduceOps.H" +#define OLD_COMBINE_REDUCE 1 + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -65,31 +67,59 @@ tmp > GlobalPointPatchField // Create the global list and insert local values if (globalPointPatch_.globalPointSize() > 0) { + // Get addressing + const labelList& sharedPointAddr = + globalPointPatch_.sharedPointAddr(); + + const Field& pField = tpField(); + + // Prepare result + tmp > tlpf(new Field(sharedPointAddr.size())); + Field& lpf = tlpf(); + +# ifdef OLD_COMBINE_REDUCE + Field gpf ( globalPointPatch_.globalPointSize(), pTraits::zero ); - const labelList& addr = globalPointPatch_.sharedPointAddr(); - const Field& pField = tpField(); - - forAll (addr, i) + forAll (sharedPointAddr, i) { - gpf[addr[i]] = pField[i]; + gpf[sharedPointAddr[i]] = pField[i]; } combineReduce(gpf, plusEqOp >()); // Extract local data - tmp > tlpf(new Field(addr.size())); - Field& lpf = tlpf(); - - forAll (addr, i) + forAll (sharedPointAddr, i) { - lpf[i] = gpf[addr[i]]; + lpf[i] = gpf[sharedPointAddr[i]]; } +# else + + // Pack data into a map + Map dataMap; + + forAll (sharedPointAddr, i) + { + dataMap.insert(sharedPointAddr[i], pField[i]); + } + + // Communicate map + Pstream::mapCombineGather(dataMap, plusEqOp()); + Pstream::mapCombineScatter(dataMap); + + // Extract local data + forAll (sharedPointAddr, i) + { + lpf[i] = dataMap[sharedPointAddr[i]]; + } + +# endif + return tlpf; } else @@ -126,6 +156,19 @@ tmp > GlobalPointPatchField { if (globalPointPatch_.globalEdgeSize() > 0) { + // Bug fix: use map-based communication. HJ, 18/Nov/2010 + + const labelList& sharedEdgeAddr = + globalPointPatch_.sharedEdgeAddr(); + + const Field& eField = teField(); + + // Prepare result + tmp > tlef(new Field(sharedEdgeAddr.size())); + Field& lef = tlef(); + +# ifdef OLD_COMBINE_REDUCE + // Create the global list and insert local values Field gef ( @@ -133,25 +176,41 @@ tmp > GlobalPointPatchField pTraits::zero ); - const labelList& addr = globalPointPatch_.sharedEdgeAddr(); - const Field& eField = teField(); - - forAll (addr, i) + forAll (sharedEdgeAddr, i) { - gef[addr[i]] = eField[i]; + gef[sharedEdgeAddr[i]] = eField[i]; } combineReduce(gef, plusEqOp >()); // Extract local data - tmp > tlef(new Field(addr.size())); - Field& lef = tlef(); - - forAll (addr, i) + forAll (sharedEdgeAddr, i) { - lef[i] = gef[addr[i]]; + lef[i] = gef[sharedEdgeAddr[i]]; } +# else + + // Pack data into a map + Map dataMap; + + forAll (sharedEdgeAddr, i) + { + dataMap.insert(sharedEdgeAddr[i], eField[i]); + } + + // Communicate map + Pstream::mapCombineGather(dataMap, plusEqOp()); + Pstream::mapCombineScatter(dataMap); + + // Extract local data + forAll (sharedEdgeAddr, i) + { + lef[i] = dataMap[sharedEdgeAddr[i]]; + } + +# endif + return tlef; } else @@ -189,6 +248,7 @@ void GlobalPointPatchField // Set the values from the global sum tmp > trpf = reduceExtractPoint(patchInternalField(pField)); + Field& rpf = trpf(); // Get addressing @@ -526,6 +586,8 @@ void GlobalPointPatchField } // Communicate map + // Note: Cannot use reduceExtract, because it uses plusEqOp + // HJ, 14/Jan/2011 Pstream::mapCombineGather(dataMap, eqOp()); Pstream::mapCombineScatter(dataMap); @@ -1063,7 +1125,8 @@ void GlobalPointPatchField { // Owner side localMult[doubleCutOwner[edgeI]] += - cutMask[coeffI]*coeffs[coeffI]*psiInternal[U[doubleCut[edgeI]]]; + cutMask[coeffI]*coeffs[coeffI]* + psiInternal[U[doubleCut[edgeI]]]; sumOffDiag[doubleCutOwner[edgeI]] += cutMask[coeffI]*coeffs[coeffI]; @@ -1072,7 +1135,8 @@ void GlobalPointPatchField // Neighbour side localMult[doubleCutNeighbour[edgeI]] += - cutMask[coeffI]*coeffs[coeffI]*psiInternal[L[doubleCut[edgeI]]]; + cutMask[coeffI]*coeffs[coeffI]* + psiInternal[L[doubleCut[edgeI]]]; sumOffDiag[doubleCutNeighbour[edgeI]] += cutMask[coeffI]*coeffs[coeffI]; @@ -1085,6 +1149,7 @@ void GlobalPointPatchField tmp > trpf = reduceExtractPoint(localMult); + Field& rpf = trpf(); // Get addressing diff --git a/src/OpenFOAM/include/OSspecific.H b/src/OpenFOAM/include/OSspecific.H index 136294f5f..045768e18 100644 --- a/src/OpenFOAM/include/OSspecific.H +++ b/src/OpenFOAM/include/OSspecific.H @@ -103,11 +103,11 @@ bool chDir(const fileName& dir); // // @return the full path name or fileName() if the name cannot be found // Optionally abort if the file cannot be found -fileName findEtcFile(const fileName&, bool mandatory=false); +fileName findEtcFile(const fileName&, bool mandatory = false); //- Make a directory and return an error if it could not be created // and does not already exist -bool mkDir(const fileName&, mode_t=0777); +bool mkDir(const fileName&, mode_t = 0777); //- Set the file mode bool chMod(const fileName&, const mode_t); @@ -120,14 +120,14 @@ fileName::Type type(const fileName&); //- Does the name exist (as DIRECTORY or FILE) in the file system? // Optionally enable/disable check for gzip file. -bool exists(const fileName&, const bool checkGzip=true); +bool exists(const fileName&, const bool checkGzip = true); //- Does the name exist as a DIRECTORY in the file system? bool isDir(const fileName&); //- Does the name exist as a FILE in the file system? // Optionally enable/disable check for gzip file. -bool isFile(const fileName&, const bool checkGzip=true); +bool isFile(const fileName&, const bool checkGzip = true); //- Return size of file off_t fileSize(const fileName&); @@ -139,8 +139,8 @@ time_t lastModified(const fileName&); fileNameList readDir ( const fileName&, - const fileName::Type=fileName::FILE, - const bool filtergz=true + const fileName::Type = fileName::FILE, + const bool filtergz = true ); //- Copy, recursively if necessary, the source to the destination @@ -172,7 +172,7 @@ void fdClose(const int); bool ping(const word&, const label port, const label timeOut); //- Check if machine is up by pinging port 22 (ssh) and 222 (rsh) -bool ping(const word&, const label timeOut=10); +bool ping(const word&, const label timeOut = 10); //- Execute the specified command int system(const string& command); diff --git a/src/OpenFOAM/meshes/meshShapes/cell/oppositeCellFace.C b/src/OpenFOAM/meshes/meshShapes/cell/oppositeCellFace.C index 9e6418dd5..5b6a23728 100644 --- a/src/OpenFOAM/meshes/meshShapes/cell/oppositeCellFace.C +++ b/src/OpenFOAM/meshes/meshShapes/cell/oppositeCellFace.C @@ -99,7 +99,7 @@ Foam::label Foam::cell::opposingFaceLabel { // There has already been an opposite face. // Non-prismatic cell - Info<< "Multiple faces not sharing vertex: " + Info<< "Multiple faces not sharing vertex: " << oppositeFaceLabel << " and " << curFaceLabels[faceI] << endl; return -1; diff --git a/src/OpenFOAM/meshes/meshShapes/face/oppositeFace.H b/src/OpenFOAM/meshes/meshShapes/face/oppositeFace.H index 892cb4b2f..241000237 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/oppositeFace.H +++ b/src/OpenFOAM/meshes/meshShapes/face/oppositeFace.H @@ -76,7 +76,6 @@ public: face(f), masterIndex_(masterIndex), oppositeIndex_(oppositeIndex) - {} diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C index a9084a23f..21907baea 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C @@ -260,26 +260,30 @@ void Foam::cyclicPolyPatch::calcTransforms() // Dump transformed first half if (debug) { - fileName fvPath(boundaryMesh().mesh().time().path()/"VTK"); - - pointField transformPoints = half0.localPoints(); - - forAll (transformPoints, pointI) + if (reverseT_.size() > 0) { - transformPoints[pointI] = - Foam::transform(reverseT_[0], transformPoints[pointI]); + fileName fvPath(boundaryMesh().mesh().time().path()/"VTK"); + + pointField transformPoints = half0.localPoints(); + + forAll (transformPoints, pointI) + { + transformPoints[pointI] = + Foam::transform(reverseT_[0], transformPoints[pointI]); + } + + standAlonePatch transformHalf0 + ( + half0.localFaces(), + transformPoints + ); + + fileName nm2(fvPath/name() + "_transform_half0_faces"); + Pout<< "cyclicPolyPatch::calcTransforms : Writing " + << "transform_half0 faces to file " << nm2 << endl; + + transformHalf0.writeVTK(nm2, transformHalf0, transformPoints); } - - standAlonePatch transformHalf0 - ( - half0.localFaces(), - transformPoints - ); - - fileName nm2(fvPath/name() + "_transform_half0_faces"); - Pout<< "cyclicPolyPatch::calcTransforms : Writing transform_half0" - << " faces to file " << nm2 << endl; - transformHalf0.writeVTK(nm2, transformHalf0, transformPoints); } // Check for error in face matching @@ -354,31 +358,33 @@ void Foam::cyclicPolyPatch::calcTransforms() } else { - maxDistance = - Foam::max - ( - maxDistance, - mag - ( - half0Ctrs[faceI] - - half1Ctrs[faceI] - ) - ); + // Disable checking for translational distance + // HJ, 13/Jan/2011 +// maxDistance = +// Foam::max +// ( +// maxDistance, +// mag +// ( +// half0Ctrs[faceI] +// - half1Ctrs[faceI] +// ) +// ); - maxRelDistance = - Foam::max - ( - maxRelDistance, - mag - ( - half0Ctrs[faceI] - - half1Ctrs[faceI] - ) - /( - mag(half1Ctrs[faceI] - half0Ctrs[faceI]) - + SMALL - ) - ); +// maxRelDistance = +// Foam::max +// ( +// maxRelDistance, +// mag +// ( +// half0Ctrs[faceI] +// - half1Ctrs[faceI] +// ) +// /( +// mag(half1Ctrs[faceI] - half0Ctrs[faceI]) +// + SMALL +// ) +// ); } } diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H index a4600f373..562acd41f 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H @@ -239,15 +239,6 @@ public: //- Attach regions void detach() const; - //- Is this the master side? - bool master() const; - - //- Is this the slave side? - bool slave() const - { - return !master(); - } - //- Return shadow patch const regionCouplePolyPatch& shadow() const; diff --git a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C index a8a68269f..9bec21322 100644 --- a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C +++ b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C @@ -217,8 +217,6 @@ void Foam::faceZone::calcCellLayers() const mc[faceI] = curMc; sc[faceI] = curSc; } - //Info << "masterCells: " << mc << endl; - //Info << "slaveCells: " << sc << endl; } } diff --git a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H index b32c43ce8..174d282db 100644 --- a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H +++ b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H @@ -263,8 +263,8 @@ public: } //- Map storing the local face index for every global face index. - // Used to find out the index of face in the zone from the known global - // face index. If the face is not in the zone, returns -1 + // Used to find out the index of face in the zone from the known + // global face index. If the face is not in the zone, returns -1 label whichFace(const label globalFaceID) const; //- Return reference to primitive patch diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/attachDetach.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/attachDetach.C index a4f270e2a..a02b01f85 100644 --- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/attachDetach.C +++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/attachDetach.C @@ -238,7 +238,7 @@ void Foam::attachDetach::checkDefinition() { FatalErrorIn("void attachDetach::checkDefinition()") << "Master and slave patch share " << nSharedPoints - << " point. This is not allowed." << nl + << " point. This is not allowed." << nl << "Please check mesh for topological errors." << abort(FatalError); } diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/attachInterface.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/attachInterface.C index c2e6073ac..2a0a75963 100644 --- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/attachInterface.C +++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/attachInterface.C @@ -109,7 +109,7 @@ void Foam::attachDetach::attachInterface polyModifyFace ( faces[masterPatchStart + faceI], // modified face - masterPatchStart + faceI, // label of face being modified + masterPatchStart + faceI, // label of face being modified masterFaceCells[faceI], // owner slaveFaceCells[faceI], // neighbour false, // face flip @@ -128,7 +128,7 @@ void Foam::attachDetach::attachInterface polyModifyFace ( faces[masterPatchStart + faceI].reverseFace(), // mod face - masterPatchStart + faceI, // label of face being modified + masterPatchStart + faceI, // label of face being modified slaveFaceCells[faceI], // owner masterFaceCells[faceI], // neighbour true, // face flip @@ -204,7 +204,7 @@ void Foam::attachDetach::attachInterface mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID) ]; } - + // Modify the face label curNeighbour = -1; @@ -222,7 +222,7 @@ void Foam::attachDetach::attachInterface own[curFaceID], // owner curNeighbour, // neighbour false, // face flip - mesh.boundaryMesh().whichPatch(curFaceID),// patch for face + mesh.boundaryMesh().whichPatch(curFaceID), // patch for face false, // remove from zone modifiedFaceZone, // zone for face modifiedFaceZoneFlip // face flip in zone @@ -251,7 +251,7 @@ void Foam::attachDetach::modifyMotionPoints if (debug) { - Pout<< "void attachDetach::modifyMotionPoints(" + Pout<< "void attachDetach::modifyMotionPoints(" << "pointField& motionPoints) const " << " for object " << name() << " : " << "Adjusting motion points." << endl; diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/detachInterface.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/detachInterface.C index e663508fc..aade9cc7b 100644 --- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/detachInterface.C +++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/attachDetach/detachInterface.C @@ -42,9 +42,9 @@ void Foam::attachDetach::detachInterface // 2. Modify all faces of the master zone, by putting them into the master // patch (look for orientation) and their renumbered mirror images // into the slave patch - // 3. Create a point renumbering list, giving a new point index for original - // points in the face patch - // 4. Grab all faces in cells on the master side and renumber them + // 3. Create a point renumbering list, giving a new point index for + // original points in the face patch + // 4. Grab all faces in cells on the master side and renumber them // using the point renumbering list. Exclude the ones that belong to // the master face zone // @@ -70,7 +70,9 @@ void Foam::attachDetach::detachInterface const polyMesh& mesh = topoChanger().mesh(); const faceZoneMesh& zoneMesh = mesh.faceZones(); - const primitiveFacePatch& masterFaceLayer = zoneMesh[faceZoneID_.index()](); + const primitiveFacePatch& masterFaceLayer = + zoneMesh[faceZoneID_.index()](); + const pointField& points = mesh.points(); const labelListList& meshEdgeFaces = mesh.edgeFaces(); @@ -88,7 +90,12 @@ void Foam::attachDetach::detachInterface // with their original labels to stop duplication label nIntEdges = masterFaceLayer.nInternalEdges(); - for (label curEdgeID = nIntEdges; curEdgeID < meshEdges.size(); curEdgeID++) + for + ( + label curEdgeID = nIntEdges; + curEdgeID < meshEdges.size(); + curEdgeID++ + ) { const labelList& curFaces = meshEdgeFaces[meshEdges[curEdgeID]]; @@ -374,6 +381,33 @@ void Foam::attachDetach::detachInterface // If the face has changed, create a modification entry if (changed) { + // Get zone ID and flipMap for the face + // Bug fix. Henrik Rusche, 20/Jan/2011 + const label oldZoneID = zoneMesh.whichZone(curFaceID); + bool oldFlip = false; + + if (oldZoneID > -1) + { + const label oldFaceInZoneID = + zoneMesh[oldZoneID].whichFace(curFaceID); + + if (oldFaceInZoneID > -1) + { + oldFlip = zoneMesh[oldZoneID].flipMap()[oldFaceInZoneID]; + } + else + { + FatalErrorIn + ( + "attachDetach::detachInterface\n" + "(\n" + " polyTopoChange& ref\n" + ") const\n" + ) << "Error in zone access." + << abort(FatalError); + } + } + if (mesh.isInternalFace(curFaceID)) { // No need to check for nei index: internal face. @@ -389,8 +423,8 @@ void Foam::attachDetach::detachInterface false, // flip flux -1, // patch for face false, // remove from zone - -1, // zone for face - false // face zone flip + oldZoneID, // zone for face + oldFlip // face zone flip ) ); // Pout << "modifying stick-out face. Internal Old face: " << oldFace << " new face: " << newFace << " own: " << own[curFaceID] << " nei: " << nei[curFaceID] << endl; @@ -408,12 +442,12 @@ void Foam::attachDetach::detachInterface false, // flip flux mesh.boundaryMesh().whichPatch(curFaceID), // patch false, // remove from zone - -1, // zone for face - false // face zone flip + oldZoneID, // zone for face + oldFlip // face zone flip ) - ); + ); // Pout << "modifying stick-out face. Boundary Old face: " << oldFace << " new face: " << newFace << " own: " << own[curFaceID] << " patch: " << mesh.boundaryMesh().whichPatch(curFaceID) << endl; - } + } } } diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/layerAdditionRemoval.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/layerAdditionRemoval.C index 436bd1d90..163f367c4 100644 --- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/layerAdditionRemoval.C +++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/layerAdditionRemoval.C @@ -217,6 +217,18 @@ Foam::layerAdditionRemoval::~layerAdditionRemoval() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +void Foam::layerAdditionRemoval::setRemoval() +{ + triggerRemoval_ = topoChanger().morphIndex(); +} + + +void Foam::layerAdditionRemoval::setAddition() +{ + triggerAddition_ = topoChanger().morphIndex(); +} + + bool Foam::layerAdditionRemoval::changeTopology() const { // Protect from multiple calculation in the same time-step diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/layerAdditionRemoval.H b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/layerAdditionRemoval.H index ce55d0d18..1e932ffff 100644 --- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/layerAdditionRemoval.H +++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/layerAdditionRemoval.H @@ -170,6 +170,12 @@ public: // Member Functions + //- Set layer removal + void setRemoval(); + + //- Set addition + void setAddition(); + //- Check for topology change virtual bool changeTopology() const; diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C index 08a8cec38..e306ccf23 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C @@ -168,7 +168,8 @@ void regionCoupleFvPatchField::evaluate const Pstream::commsTypes ) { - // Implement weights-based stabilised harmonic interpolation using magnitude + // Implement weights-based stabilised harmonic interpolation using + // magnitude of type // Algorithm: // 1) calculate magnitude of internal field and neighbour field // 2) calculate harmonic mean magnitude @@ -189,7 +190,7 @@ void regionCoupleFvPatchField::evaluate forAll (weights, faceI) { - den = (mOwn[faceI] - mNei[faceI]); + den = mOwn[faceI] - mNei[faceI]; if (mag(den) > SMALL) { diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C index 556f93363..86d303a7a 100644 --- a/src/finiteVolume/fvMesh/fvMesh.C +++ b/src/finiteVolume/fvMesh/fvMesh.C @@ -39,20 +39,12 @@ License #include "mapClouds.H" #include "volPointInterpolation.H" -#include "extendedLeastSquaresVectors.H" -#include "extendedLeastSquaresVectors.H" -#include "leastSquaresVectors.H" -#include "CentredFitData.H" -#include "linearFitPolynomial.H" -#include "quadraticFitPolynomial.H" -#include "quadraticLinearFitPolynomial.H" -//#include "quadraticFitSnGradData.H" -#include "skewCorrectionVectors.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // defineTypeNameAndDebug(Foam::fvMesh, 0); + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // void Foam::fvMesh::clearGeomNotOldVol() diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDTVD.H b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDTVD.H index b622447d0..af5fd4738 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDTVD.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDTVD.H @@ -113,9 +113,14 @@ public: } // Stabilise for division - gradf = stabilise(gradf, VSMALL); + // Changed to SMALL to prevent FPE. OB, 14/Jan/2011 + gradf = stabilise(gradf, SMALL); - return 2*(gradcf/gradf) - 1; + // New formulation. Oliver Borm and Aleks Jemcov + // HJ, 13/Jan/2011 + return max(2*(gradcf/gradf) - 1, 0); + +// return 2*(gradcf/gradf) - 1; } }; diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDVTVDV.H b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDVTVDV.H index e052b9019..287788a2c 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDVTVDV.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDVTVDV.H @@ -114,9 +114,14 @@ public: } // Stabilise for division - gradf = stabilise(gradf, VSMALL); + //Changed to SMALL to prevent FPE. OB, 14/Jan/2011 + gradf = stabilise(gradf, SMALL); - return 2*(gradcf/gradf) - 1; + // New formulation. Oliver Borm and Aleks Jemcov + // HJ, 13/Jan/2011 + return max(2*(gradcf/gradf) - 1, 0); + +// return 2*(gradcf/gradf) - 1; } }; diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/upwind/upwind.H b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/upwind/upwind.H index 1c4ea7048..7cb10c0fd 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/upwind/upwind.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/upwind/upwind.H @@ -121,7 +121,7 @@ public: this->mesh() ), this->mesh(), - dimensionedScalar("upwindLimiter", dimless, 0.0) + dimless ) ); } diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/vanAlbada/vanAlbada.H b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/vanAlbada/vanAlbada.H index e3d7556f1..14524c01c 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/vanAlbada/vanAlbada.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/vanAlbada/vanAlbada.H @@ -78,7 +78,11 @@ public: faceFlux, phiP, phiN, gradcP, gradcN, d ); - return r*(r + 1)/(sqr(r) + 1); + // New formulation. Oliver Borm and Aleks Jemcov + // HJ, 13/Jan/2011 + return (r + 1)/(r + 1/stabilise(r, VSMALL)); + +// return r*(r + 1)/(sqr(r) + 1); } }; diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.C index 16d36fdd6..364fbff66 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.C @@ -34,15 +34,7 @@ Description namespace Foam { - -defineTypeNameAndDebug(harmonic, 0); - -surfaceInterpolationScheme::addMeshFluxConstructorToTable - addharmonicScalarMeshFluxConstructorToTable_; - -surfaceInterpolationScheme::addMeshConstructorToTable - addharmonicScalarMeshConstructorToTable_; - + makeSurfaceInterpolationScheme(harmonic) } // ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H index ba6e591f8..b7d2e5df5 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H @@ -42,7 +42,6 @@ SourceFiles #include "surfaceInterpolationScheme.H" #include "volFields.H" #include "surfaceFields.H" -#include "reverseLinear.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -53,9 +52,10 @@ namespace Foam Class harmonic Declaration \*---------------------------------------------------------------------------*/ +template class harmonic : - public surfaceInterpolationScheme + public surfaceInterpolationScheme { // Private Member Functions @@ -74,10 +74,10 @@ public: //- Construct from mesh harmonic(const fvMesh& mesh) : - surfaceInterpolationScheme(mesh) + surfaceInterpolationScheme(mesh) {} - //- Construct from Istream. + //- Construct from Istream // The name of the flux field is read from the Istream and looked-up // from the mesh objectRegistry harmonic @@ -86,7 +86,7 @@ public: Istream& is ) : - surfaceInterpolationScheme(mesh) + surfaceInterpolationScheme(mesh) {} //- Construct from faceFlux and Istream @@ -97,7 +97,7 @@ public: Istream& is ) : - surfaceInterpolationScheme(mesh) + surfaceInterpolationScheme(mesh) {} @@ -106,26 +106,60 @@ public: //- Return the interpolation weighting factors virtual tmp weights ( - const GeometricField& + const GeometricField& phi ) const { - notImplemented + tmp tw ( - "harmonic::weights" - "(const GeometricField&)" + new surfaceScalarField + ( + IOobject + ( + "harmonicWeightingFactors" + phi.name(), + this->mesh().time().timeName(), + this->mesh() + ), + this->mesh() , + dimless + ) ); - return tmp(NULL); - } + surfaceScalarField& w = tw(); - //- Return the face-interpolate of the given cell field - virtual tmp > - interpolate - ( - const GeometricField& vf - ) const - { - return 1.0/(reverseLinear(vf.mesh()).interpolate(1.0/vf)); + const unallocLabelList& owner = this->mesh().owner(); + const unallocLabelList& neighbour = this->mesh().neighbour(); + + scalarField magPhi = mag(phi); + + // Initialise weights to 0.5 for uniform field (den = 0) + scalarField& wIn = w.internalField(); + wIn = 0.5; + + // Calculate internal weights using field magnitude + scalar mOwn, mNei, den, mean; + + forAll (owner, faceI) + { + mOwn = magPhi[owner[faceI]]; + mNei = magPhi[neighbour[faceI]]; + + mean = 2*(mOwn*mNei)/(mOwn + mNei + SMALL); + den = mOwn - mNei; + + if (mag(den) > SMALL) + { + wIn[faceI] = (mean - mNei)/den; + } + else + { + // Use 0.5 weights + } + } + + // Boundary weights are 1 + w.boundaryField() = 1; + + return tw; } }; diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linear/linear.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linear/linear.H index e458af55b..4ffabe77f 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linear/linear.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linear/linear.H @@ -120,7 +120,7 @@ template tmp > linearInterpolate(const tmp >& tvf) { - tmp > tinterp = + tmp > tinterp = linearInterpolate(tvf()); tvf.clear(); return tinterp; diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C index a3595bab6..721ca9f8f 100644 --- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C +++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C @@ -98,7 +98,9 @@ void Foam::directMappedPatchBase::collectSamples labelListList globalFaces(Pstream::nProcs()); globalFc[Pstream::myProcNo()] = patch_.faceCentres(); - globalSamples[Pstream::myProcNo()] = globalFc[Pstream::myProcNo()]+offsets_; + globalSamples[Pstream::myProcNo()] = + globalFc[Pstream::myProcNo()] + offsets_; + globalFaces[Pstream::myProcNo()] = identity(patch_.size()); // Distribute to all processors @@ -115,11 +117,13 @@ void Foam::directMappedPatchBase::collectSamples globalSamples, accessOp() ); + patchFaces = ListListOps::combine ( globalFaces, accessOp() ); + patchFc = ListListOps::combine ( globalFc, @@ -135,8 +139,9 @@ void Foam::directMappedPatchBase::collectSamples accessOp() ) ); + label sampleI = 0; - forAll(nPerProc, procI) + forAll (nPerProc, procI) { for (label i = 0; i < nPerProc[procI]; i++) { @@ -173,13 +178,14 @@ void Foam::directMappedPatchBase::findSamples "directMappedPatchBase::findSamples(const pointField&," " labelList&, labelList&, pointField&) const" ) << "No need to supply a patch name when in " - << sampleModeNames_[mode_] << " mode." << exit(FatalError); + << sampleModeNames_[mode_] << " mode." + << abort(FatalError); } // Octree based search engine meshSearch meshSearchEngine(mesh, false); - forAll(samples, sampleI) + forAll (samples, sampleI) { const point& sample = samples[sampleI]; @@ -200,7 +206,7 @@ void Foam::directMappedPatchBase::findSamples cc, cellI ); - nearest[sampleI].second().first() = magSqr(cc-sample); + nearest[sampleI].second().first() = magSqr(cc - sample); nearest[sampleI].second().second() = Pstream::myProcNo(); } } @@ -215,7 +221,7 @@ void Foam::directMappedPatchBase::findSamples if (pp.empty()) { - forAll(samples, sampleI) + forAll (samples, sampleI) { nearest[sampleI].second().first() = Foam::sqr(GREAT); nearest[sampleI].second().second() = Pstream::myProcNo(); @@ -251,7 +257,7 @@ void Foam::directMappedPatchBase::findSamples 3.0 // duplicity ); - forAll(samples, sampleI) + forAll (samples, sampleI) { const point& sample = samples[sampleI]; @@ -290,13 +296,14 @@ void Foam::directMappedPatchBase::findSamples "directMappedPatchBase::findSamples(const pointField&," " labelList&, labelList&, pointField&) const" ) << "No need to supply a patch name when in " - << sampleModeNames_[mode_] << " mode." << exit(FatalError); + << sampleModeNames_[mode_] << " mode." + << abort(FatalError); } // Octree based search engine meshSearch meshSearchEngine(mesh, false); - forAll(samples, sampleI) + forAll (samples, sampleI) { const point& sample = samples[sampleI]; @@ -340,20 +347,21 @@ void Foam::directMappedPatchBase::findSamples { Info<< "directMappedPatchBase::findSamples on mesh " << sampleRegion_ << " : " << endl; - forAll(nearest, sampleI) + + forAll (nearest, sampleI) { label procI = nearest[sampleI].second().second(); label localI = nearest[sampleI].first().index(); - Info<< " " << sampleI << " coord:"<< samples[sampleI] - << " found on processor:" << procI - << " in local cell/face:" << localI - << " with cc:" << nearest[sampleI].first().rawPoint() << endl; + Info<< " " << sampleI << " coord: "<< samples[sampleI] + << " found on processor: " << procI + << " in local cell/face: " << localI + << " with cc: " << nearest[sampleI].first().rawPoint() << endl; } } // Check for samples not being found - forAll(nearest, sampleI) + forAll (nearest, sampleI) { if (!nearest[sampleI].first().hit()) { @@ -363,8 +371,8 @@ void Foam::directMappedPatchBase::findSamples "(const pointField&, labelList&" ", labelList&, pointField&)" ) << "Did not find sample " << samples[sampleI] - << " on any processor of region" << sampleRegion_ - << exit(FatalError); + << " on any processor of region " << sampleRegion_ + << abort(FatalError); } } @@ -374,7 +382,7 @@ void Foam::directMappedPatchBase::findSamples sampleIndices.setSize(samples.size()); sampleLocations.setSize(samples.size()); - forAll(nearest, sampleI) + forAll (nearest, sampleI) { sampleProcs[sampleI] = nearest[sampleI].second().second(); sampleIndices[sampleI] = nearest[sampleI].first().index(); @@ -388,7 +396,8 @@ void Foam::directMappedPatchBase::calcMapping() const if (mapPtr_.valid()) { FatalErrorIn("directMappedPatchBase::calcMapping() const") - << "Mapping already calculated" << exit(FatalError); + << "Mapping already calculated" + << abort(FatalError); } if @@ -415,7 +424,7 @@ void Foam::directMappedPatchBase::calcMapping() const } - // Get global list of all samples and the processor and face they come from. + // Get global list of all samples and the processor and face they come from pointField samples; labelList patchFaceProcs; labelList patchFaces; @@ -435,7 +444,7 @@ void Foam::directMappedPatchBase::calcMapping() const // - cell/face sample is in (so source when mapping) // sampleIndices, sampleProcs. - //forAll(samples, i) + //forAll (samples, i) //{ // Info<< i << " need data in region " // << patch_.boundaryMesh().mesh().name() @@ -464,7 +473,7 @@ void Foam::directMappedPatchBase::calcMapping() const label vertI = 0; - forAll(patchFc, i) + forAll (patchFc, i) { meshTools::writeOBJ(str, patchFc[i]); vertI++; @@ -482,7 +491,7 @@ void Foam::directMappedPatchBase::calcMapping() const // const scalarField magOffset(mag(sampleLocations - patchFc)); // const scalar avgOffset(average(magOffset)); // - // forAll(magOffset, sampleI) + // forAll (magOffset, sampleI) // { // if // ( @@ -518,7 +527,7 @@ void Foam::directMappedPatchBase::calcMapping() const labelListList& subMap = mapPtr_().subMap(); labelListList& constructMap = mapPtr_().constructMap(); - forAll(subMap, procI) + forAll (subMap, procI) { subMap[procI] = UIndirectList