From 73b83058b191daefe60eb2170e49a756b0a614c0 Mon Sep 17 00:00:00 2001 From: henry Date: Thu, 16 Oct 2008 14:56:30 +0100 Subject: [PATCH] Added new extended stencil handling from Mattijs. --- src/OpenFOAM/meshes/MeshObject/MeshObject.C | 20 + src/OpenFOAM/meshes/MeshObject/MeshObject.H | 3 + .../mapPolyMesh/mapDistribute/mapDistribute.C | 140 +++++ .../mapPolyMesh/mapDistribute/mapDistribute.H | 22 + .../mapDistribute/mapDistributeTemplates.C | 353 +++++++++++- src/finiteVolume/Make/files | 24 +- .../extendedStencil/extendedCentredStencil.C | 65 +++ .../extendedStencil/extendedCentredStencil.H | 137 +++++ .../fvMesh/extendedStencil/extendedStencil.C | 312 ++++------- .../fvMesh/extendedStencil/extendedStencil.H | 69 +-- .../centredCECStencilObject.C | 38 ++ .../centredCECStencilObject.H | 88 +++ .../centredCFCStencilObject.C | 38 ++ .../centredCFCStencilObject.H | 88 +++ .../centredCPCStencilObject.C | 38 ++ .../centredCPCStencilObject.H | 88 +++ .../centredFECStencilObject.C | 38 ++ .../centredFECStencilObject.H | 88 +++ .../upwindCECStencilObject.C | 38 ++ .../upwindCECStencilObject.H | 89 +++ .../upwindCFCStencilObject.C | 38 ++ .../upwindCFCStencilObject.H | 89 +++ .../upwindCPCStencilObject.C | 38 ++ .../upwindCPCStencilObject.H | 89 +++ .../upwindFECStencilObject.C | 38 ++ .../upwindFECStencilObject.H | 89 +++ .../extendedStencilTemplates.C | 48 +- .../extendedStencil/extendedUpwindStencil.C | 429 +++++++++++++++ .../extendedStencil/extendedUpwindStencil.H | 177 ++++++ .../extendedUpwindStencilTemplates.C | 138 +++++ .../faceStencil/cellEdgeCellStencil.C | 209 +++++++ .../faceStencil/cellEdgeCellStencil.H | 95 ++++ .../faceStencil/cellFaceCellStencil.C | 167 ++++++ .../faceStencil/cellFaceCellStencil.H | 83 +++ .../faceStencil/cellPointCellStencil.C | 174 ++++++ .../faceStencil/cellPointCellStencil.H | 94 ++++ .../faceStencil/faceEdgeCellStencil.C | 445 +++++++++++++++ .../faceStencil/faceEdgeCellStencil.H | 96 ++++ .../extendedStencil/faceStencil/faceStencil.C | 512 ++++++++++++++++++ .../extendedStencil/faceStencil/faceStencil.H | 160 ++++++ .../schemes/linearFit/linearFit.C | 36 ++ .../schemes/linearFit/linearFit.H | 138 +++++ .../schemes/linearFit/linearFitData.C | 371 +++++++++++++ .../schemes/linearFit/linearFitData.H | 140 +++++ .../schemes/quadraticFit/quadraticFit.H | 12 +- .../schemes/quadraticFit/quadraticFitData.C | 8 +- .../schemes/quadraticFit/quadraticFitData.H | 13 +- src/sampling/probes/probes.C | 2 +- src/sampling/probes/probesTemplates.C | 2 +- 49 files changed, 5397 insertions(+), 309 deletions(-) create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedCentredStencil.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedCentredStencil.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCECStencilObject.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCECStencilObject.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCPCStencilObject.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCPCStencilObject.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredFECStencilObject.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredFECStencilObject.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCECStencilObject.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCECStencilObject.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCFCStencilObject.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCFCStencilObject.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCPCStencilObject.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCPCStencilObject.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindFECStencilObject.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindFECStencilObject.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencil.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencil.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellEdgeCellStencil.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellEdgeCellStencil.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellPointCellStencil.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellPointCellStencil.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/faceStencil/faceEdgeCellStencil.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/faceStencil/faceEdgeCellStencil.H create mode 100644 src/finiteVolume/fvMesh/extendedStencil/faceStencil/faceStencil.C create mode 100644 src/finiteVolume/fvMesh/extendedStencil/faceStencil/faceStencil.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.H create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.H diff --git a/src/OpenFOAM/meshes/MeshObject/MeshObject.C b/src/OpenFOAM/meshes/MeshObject/MeshObject.C index c089e0e531..fddfadfa59 100644 --- a/src/OpenFOAM/meshes/MeshObject/MeshObject.C +++ b/src/OpenFOAM/meshes/MeshObject/MeshObject.C @@ -83,6 +83,26 @@ const Type& Foam::MeshObject::New } +template +template +const Type& Foam::MeshObject::New +( + const Mesh& mesh, + const Data1& d1, + const Data2& d2 +) +{ + if (!mesh.db().objectRegistry::foundObject(Type::typeName)) + { + return store(new Type(mesh, d1, d2)); + } + else + { + return mesh.db().objectRegistry::lookupObject(Type::typeName); + } +} + + // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * // template diff --git a/src/OpenFOAM/meshes/MeshObject/MeshObject.H b/src/OpenFOAM/meshes/MeshObject/MeshObject.H index 595f2f1df4..b142b47f96 100644 --- a/src/OpenFOAM/meshes/MeshObject/MeshObject.H +++ b/src/OpenFOAM/meshes/MeshObject/MeshObject.H @@ -72,6 +72,9 @@ public: template static const Type& New(const Mesh& mesh, const Data& d); + template + static const Type& New(const Mesh& mesh, const Data1&, const Data2&); + // Destructor diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C index 5421b328d8..5ef7b693a5 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C @@ -271,4 +271,144 @@ Foam::mapDistribute::mapDistribute } +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::mapDistribute::compact(const boolList& elemIsUsed) +{ + // 1. send back to sender. Have him delete the corresponding element + // from the submap and do the same to the constructMap locally + // (and in same order). + + // Send elemIsUsed field to neighbour. Use nonblocking code from + // mapDistribute but in reverse order. + { + List sendFields(Pstream::nProcs()); + + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + const labelList& map = constructMap_[domain]; + + if (domain != Pstream::myProcNo() && map.size() > 0) + { + boolList& subField = sendFields[domain]; + subField.setSize(map.size()); + forAll(map, i) + { + subField[i] = elemIsUsed[map[i]]; + } + + OPstream::write + ( + Pstream::nonBlocking, + domain, + reinterpret_cast(subField.begin()), + subField.size()*sizeof(bool) + ); + } + } + + // Set up receives from neighbours + + List recvFields(Pstream::nProcs()); + + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + const labelList& map = subMap_[domain]; + + if (domain != Pstream::myProcNo() && map.size() > 0) + { + recvFields[domain].setSize(map.size()); + IPstream::read + ( + Pstream::nonBlocking, + domain, + reinterpret_cast(recvFields[domain].begin()), + recvFields[domain].size()*sizeof(bool) + ); + } + } + + + // Set up 'send' to myself - write directly into recvFields + + { + const labelList& map = constructMap_[Pstream::myProcNo()]; + + recvFields[Pstream::myProcNo()].setSize(map.size()); + forAll(map, i) + { + recvFields[Pstream::myProcNo()][i] = elemIsUsed[map[i]]; + } + } + + + // Wait for all to finish + + OPstream::waitRequests(); + IPstream::waitRequests(); + + + // Compact out all submap entries that are referring to unused elements + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + const labelList& map = subMap_[domain]; + + labelList newMap(map.size()); + label newI = 0; + + forAll(map, i) + { + if (recvFields[domain][i]) + { + // So element is used on destination side + newMap[newI++] = map[i]; + } + } + if (newI < map.size()) + { + newMap.setSize(newI); + subMap_[domain].transfer(newMap); + } + } + } + + + // 2. remove from construct map - since end-result (element in elemIsUsed) + // not used. + + label maxConstructIndex = -1; + + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + const labelList& map = constructMap_[domain]; + + labelList newMap(map.size()); + label newI = 0; + + forAll(map, i) + { + label destinationI = map[i]; + + // Is element is used on destination side + if (elemIsUsed[destinationI]) + { + maxConstructIndex = max(maxConstructIndex, destinationI); + + newMap[newI++] = destinationI; + } + } + if (newI < map.size()) + { + newMap.setSize(newI); + constructMap_[domain].transfer(newMap); + } + } + + constructSize_ = maxConstructIndex+1; + + // Clear the schedule (note:not necessary if nothing changed) + schedulePtr_.clear(); +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H index 748834ca4d..8427670a25 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H @@ -51,6 +51,7 @@ SourceFiles #include "labelList.H" #include "labelPair.H" #include "Pstream.H" +#include "boolList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -155,6 +156,12 @@ public: // Other + //- Compact maps. Gets per field a bool whether it is used (locally) + // and works out itself what this side and sender side can remove + // from maps. + void compact(const boolList& elemIsUsed); + + //- Distribute data. Note:schedule only used for Pstream::scheduled // for now, all others just use send-to-all, receive-from-all. template @@ -168,6 +175,21 @@ public: List& ); + //- Distribute data. If multiple processors writing to same + // position adds contributions using cop. + template + static void distribute + ( + const Pstream::commsTypes commsType, + const List& schedule, + const label constructSize, + const labelListList& subMap, + const labelListList& constructMap, + List&, + const CombineOp& cop, + const T& nullValue + ); + //- Distribute data using default commsType. template void distribute(List& fld) const diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C index 5f16ad1335..8b4598bfb8 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C @@ -268,25 +268,34 @@ void Foam::mapDistribute::distribute } - // Combine bits. Note that can reuse field storage + // Set up 'send' to myself - // Subset myself - const labelList& mySubMap = subMap[Pstream::myProcNo()]; - - List subField(mySubMap.size()); - forAll(mySubMap, i) { - subField[i] = field[mySubMap[i]]; + const labelList& map = subMap[Pstream::myProcNo()]; + + List& subField = sendFields[Pstream::myProcNo()]; + subField.setSize(map.size()); + forAll(map, i) + { + subField[i] = field[map[i]]; + } } + + // Combine bits. Note that can reuse field storage + field.setSize(constructSize); - // Receive sub field from myself (subField) - const labelList& map = constructMap[Pstream::myProcNo()]; - forAll(map, i) + // Receive sub field from myself (sendFields[Pstream::myProcNo()]) { - field[map[i]] = subField[i]; + const labelList& map = constructMap[Pstream::myProcNo()]; + const List& subField = sendFields[Pstream::myProcNo()]; + + forAll(map, i) + { + field[map[i]] = subField[i]; + } } @@ -339,4 +348,326 @@ void Foam::mapDistribute::distribute } +// Distribute list. +template +void Foam::mapDistribute::distribute +( + const Pstream::commsTypes commsType, + const List& schedule, + const label constructSize, + const labelListList& subMap, + const labelListList& constructMap, + List& field, + const CombineOp& cop, + const T& nullValue +) +{ + if (commsType == Pstream::blocking) + { + // Since buffered sending can reuse the field to collect the + // received data. + + // Send sub field to neighbour + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + const labelList& map = subMap[domain]; + + if (domain != Pstream::myProcNo() && map.size() > 0) + { + List subField(map.size()); + forAll(map, i) + { + subField[i] = field[map[i]]; + } + OPstream toNbr(Pstream::blocking, domain); + toNbr << subField; + } + } + + // Subset myself + const labelList& mySubMap = subMap[Pstream::myProcNo()]; + + List subField(mySubMap.size()); + forAll(mySubMap, i) + { + subField[i] = field[mySubMap[i]]; + } + + // Receive sub field from myself (subField) + const labelList& map = constructMap[Pstream::myProcNo()]; + + field.setSize(constructSize); + field = nullValue; + + forAll(map, i) + { + cop(field[map[i]], subField[i]); + } + + // Receive sub field from neighbour + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + const labelList& map = constructMap[domain]; + + if (domain != Pstream::myProcNo() && map.size() > 0) + { + IPstream fromNbr(Pstream::blocking, domain); + List subField(fromNbr); + + if (subField.size() != map.size()) + { + FatalErrorIn + ( + "template\n" + "void mapDistribute::distribute\n" + "(\n" + " const Pstream::commsTypes commsType,\n" + " const List& schedule,\n" + " const label constructSize,\n" + " const labelListList& subMap,\n" + " const labelListList& constructMap,\n" + " List& field\n" + ")\n" + ) << "Expected from processor " << domain + << " " << map.size() << " but received " + << subField.size() << " elements." + << abort(FatalError); + } + + forAll(map, i) + { + cop(field[map[i]], subField[i]); + } + } + } + } + else if (commsType == Pstream::scheduled) + { + // Need to make sure I don't overwrite field with received data + // since the data might need to be sent to another processor. So + // allocate a new field for the results. + List newField(constructSize, nullValue); + + // Subset myself + const labelList& mySubMap = subMap[Pstream::myProcNo()]; + + List subField(mySubMap.size()); + forAll(mySubMap, i) + { + subField[i] = field[mySubMap[i]]; + } + + // Receive sub field from myself (subField) + const labelList& map = constructMap[Pstream::myProcNo()]; + + forAll(map, i) + { + cop(newField[map[i]], subField[i]); + } + + // Schedule will already have pruned 0-sized comms + forAll(schedule, i) + { + const labelPair& twoProcs = schedule[i]; + label sendProc = twoProcs[0]; + label recvProc = twoProcs[1]; + + if (Pstream::myProcNo() == sendProc) + { + // I am sender. Send to recvProc. + const labelList& map = subMap[recvProc]; + + List subField(map.size()); + forAll(map, i) + { + subField[i] = field[map[i]]; + } + + OPstream toNbr(Pstream::scheduled, recvProc); + toNbr << subField; + } + else + { + // I am receiver. Receive from sendProc. + IPstream fromNbr(Pstream::scheduled, sendProc); + List subField(fromNbr); + + const labelList& map = constructMap[sendProc]; + + if (subField.size() != map.size()) + { + FatalErrorIn + ( + "template\n" + "void mapDistribute::distribute\n" + "(\n" + " const Pstream::commsTypes commsType,\n" + " const List& schedule,\n" + " const label constructSize,\n" + " const labelListList& subMap,\n" + " const labelListList& constructMap,\n" + " List& field\n" + ")\n" + ) << "Expected from processor " << sendProc + << " " << map.size() << " but received " + << subField.size() << " elements." + << abort(FatalError); + } + + forAll(map, i) + { + cop(newField[map[i]], subField[i]); + } + } + } + field.transfer(newField); + } + else if (commsType == Pstream::nonBlocking) + { + if (!contiguous()) + { + FatalErrorIn + ( + "template\n" + "void mapDistribute::distribute\n" + "(\n" + " const Pstream::commsTypes commsType,\n" + " const List& schedule,\n" + " const label constructSize,\n" + " const labelListList& subMap,\n" + " const labelListList& constructMap,\n" + " List& field\n" + ")\n" + ) << "Non-blocking only supported for contiguous data." + << exit(FatalError); + } + + // Set up sends to neighbours + + List > sendFields(Pstream::nProcs()); + + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + const labelList& map = subMap[domain]; + + if (domain != Pstream::myProcNo() && map.size() > 0) + { + List& subField = sendFields[domain]; + subField.setSize(map.size()); + forAll(map, i) + { + subField[i] = field[map[i]]; + } + + OPstream::write + ( + Pstream::nonBlocking, + domain, + reinterpret_cast(subField.begin()), + subField.size()*sizeof(T) + ); + } + } + + // Set up receives from neighbours + + List > recvFields(Pstream::nProcs()); + + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + const labelList& map = constructMap[domain]; + + if (domain != Pstream::myProcNo() && map.size() > 0) + { + recvFields[domain].setSize(map.size()); + IPstream::read + ( + Pstream::nonBlocking, + domain, + reinterpret_cast(recvFields[domain].begin()), + recvFields[domain].size()*sizeof(T) + ); + } + } + + // Set up 'send' to myself + + { + const labelList& map = subMap[Pstream::myProcNo()]; + + List& subField = sendFields[Pstream::myProcNo()]; + subField.setSize(map.size()); + forAll(map, i) + { + subField[i] = field[map[i]]; + } + } + + + // Combine bits. Note that can reuse field storage + + field.setSize(constructSize); + field = nullValue; + + // Receive sub field from myself (subField) + { + const labelList& map = constructMap[Pstream::myProcNo()]; + const List& subField = sendFields[Pstream::myProcNo()]; + + forAll(map, i) + { + cop(field[map[i]], subField[i]); + } + } + + + // Wait for all to finish + + OPstream::waitRequests(); + IPstream::waitRequests(); + + // Collect neighbour fields + + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + const labelList& map = constructMap[domain]; + + if (domain != Pstream::myProcNo() && map.size() > 0) + { + if (recvFields[domain].size() != map.size()) + { + FatalErrorIn + ( + "template\n" + "void mapDistribute::distribute\n" + "(\n" + " const Pstream::commsTypes commsType,\n" + " const List& schedule,\n" + " const label constructSize,\n" + " const labelListList& subMap,\n" + " const labelListList& constructMap,\n" + " List& field\n" + ")\n" + ) << "Expected from processor " << domain + << " " << map.size() << " but received " + << recvFields[domain].size() << " elements." + << abort(FatalError); + } + + forAll(map, i) + { + cop(field[map[i]], recvFields[domain][i]); + } + } + } + } + else + { + FatalErrorIn("mapDistribute::distribute(..)") + << "Unknown communication schedule " << commsType + << abort(FatalError); + } +} + + // ************************************************************************* // diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 923e4c50e6..7bf27115a5 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -40,6 +40,24 @@ $(fvMeshMapper)/fvSurfaceMapper.C extendedStencil = fvMesh/extendedStencil $(extendedStencil)/extendedStencil.C +$(extendedStencil)/extendedUpwindStencil.C +$(extendedStencil)/extendedCentredStencil.C + +$(extendedStencil)/faceStencil/faceStencil.C +$(extendedStencil)/faceStencil/faceEdgeCellStencil.C +$(extendedStencil)/faceStencil/cellFaceCellStencil.C +$(extendedStencil)/faceStencil/cellPointCellStencil.C +$(extendedStencil)/faceStencil/cellEdgeCellStencil.C + +$(extendedStencil)/extendedStencilMeshObjects/centredCECStencilObject.C +$(extendedStencil)/extendedStencilMeshObjects/centredCFCStencilObject.C +$(extendedStencil)/extendedStencilMeshObjects/centredCPCStencilObject.C +$(extendedStencil)/extendedStencilMeshObjects/centredFECStencilObject.C + +$(extendedStencil)/extendedStencilMeshObjects/upwindCECStencilObject.C +$(extendedStencil)/extendedStencilMeshObjects/upwindCFCStencilObject.C +$(extendedStencil)/extendedStencilMeshObjects/upwindCPCStencilObject.C +$(extendedStencil)/extendedStencilMeshObjects/upwindFECStencilObject.C fvPatchFields = fields/fvPatchFields $(fvPatchFields)/fvPatchField/fvPatchFields.C @@ -168,6 +186,8 @@ $(schemes)/harmonic/harmonic.C $(schemes)/localBlended/localBlended.C $(schemes)/localMax/localMax.C $(schemes)/localMin/localMin.C +//$(schemes)/linearFit/linearFit.C +//$(schemes)/linearFit/linearFitData.C $(schemes)/quadraticFit/quadraticFit.C $(schemes)/quadraticFit/quadraticFitData.C @@ -251,8 +271,8 @@ $(snGradSchemes)/snGradScheme/snGradSchemes.C $(snGradSchemes)/correctedSnGrad/correctedSnGrads.C $(snGradSchemes)/limitedSnGrad/limitedSnGrads.C $(snGradSchemes)/uncorrectedSnGrad/uncorrectedSnGrads.C -$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGradData.C -$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C +//$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGradData.C +//$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C convectionSchemes = finiteVolume/convectionSchemes $(convectionSchemes)/convectionScheme/convectionSchemes.C diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedCentredStencil.C b/src/finiteVolume/fvMesh/extendedStencil/extendedCentredStencil.C new file mode 100644 index 0000000000..956e796e18 --- /dev/null +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedCentredStencil.C @@ -0,0 +1,65 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "mapDistribute.H" +#include "extendedCentredStencil.H" +#include "faceStencil.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::extendedCentredStencil::extendedCentredStencil(const faceStencil& stencil) +: + extendedStencil(stencil.mesh()) +{ + stencil_ = stencil; + + // Calculate distribute map (also renumbers stencil) + mapPtr_ = calcDistributeMap(stencil.globalNumbering(), stencil_); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +// Per face which elements of the stencil to keep. +void Foam::extendedCentredStencil::compact() +{ + boolList isInStencil(map().constructSize(), false); + + forAll(stencil_, faceI) + { + const labelList& stencilCells = stencil_[faceI]; + + forAll(stencilCells, i) + { + isInStencil[stencilCells[i]] = true; + } + } + + mapPtr_().compact(isInStencil); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedCentredStencil.H b/src/finiteVolume/fvMesh/extendedStencil/extendedCentredStencil.H new file mode 100644 index 0000000000..0ae37d5eb6 --- /dev/null +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedCentredStencil.H @@ -0,0 +1,137 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::extendedCentredStencil + +Description + +SourceFiles + extendedCentredStencil.C + +\*---------------------------------------------------------------------------*/ + +#ifndef extendedCentredStencil_H +#define extendedCentredStencil_H + +#include "extendedStencil.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class faceStencil; + +/*---------------------------------------------------------------------------*\ + Class extendedCentredStencil Declaration +\*---------------------------------------------------------------------------*/ + +class extendedCentredStencil +: + public extendedStencil +{ + // Private data + + //- Swap map for getting neigbouring data + autoPtr mapPtr_; + + //- Per face the stencil. + labelListList stencil_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + extendedCentredStencil(const extendedCentredStencil&); + + //- Disallow default bitwise assignment + void operator=(const extendedCentredStencil&); + + +public: + + // Constructors + + //- Construct from uncompacted face stencil + explicit extendedCentredStencil(const faceStencil&); + + + // Member Functions + + //- Return reference to the parallel distribution map + const mapDistribute& map() const + { + return mapPtr_(); + } + + //- Return reference to the stencil + const labelListList& stencil() const + { + return stencil_; + } + + //- After removing elements from the stencil adapt the schedule (map). + void compact(); + + //- Use map to get the data into stencil order + template + void collectData + ( + const GeometricField& fld, + List >& stencilFld + ) const + { + extendedStencil::collectData(map(), stencil(), fld, stencilFld); + } + + //- Sum vol field contributions to create face values + template + tmp > weightedSum + ( + const GeometricField& fld, + const List >& stencilWeights + ) const + { + return extendedStencil::weightedSum + ( + map(), + stencil(), + fld, + stencilWeights + ); + } + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C index bb4c6c1534..c2273d4a52 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C @@ -32,154 +32,83 @@ License // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // Calculates per face a list of global cell/face indices. -void Foam::extendedStencil::calcFaceStencils +void Foam::extendedStencil::calcFaceStencil ( - const polyMesh& mesh, - const globalIndex& globalNumbering + const labelListList& globalCellCells, + labelListList& faceStencil ) { - const polyBoundaryMesh& patches = mesh.boundaryMesh(); - const label nBnd = mesh.nFaces()-mesh.nInternalFaces(); - const labelList& own = mesh.faceOwner(); - const labelList& nei = mesh.faceNeighbour(); - - - // Determine neighbouring global cell or boundary face - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - labelList neiGlobal(nBnd); - forAll(patches, patchI) - { - const polyPatch& pp = patches[patchI]; - label faceI = pp.start(); - - if (pp.coupled()) - { - // For coupled faces get the cell on the other side - forAll(pp, i) - { - label bFaceI = faceI-mesh.nInternalFaces(); - neiGlobal[bFaceI] = globalNumbering.toGlobal(own[faceI]); - faceI++; - } - } - else if (isA(pp)) - { - forAll(pp, i) - { - label bFaceI = faceI-mesh.nInternalFaces(); - neiGlobal[bFaceI] = -1; - faceI++; - } - } - else - { - // For noncoupled faces get the boundary face. - forAll(pp, i) - { - label bFaceI = faceI-mesh.nInternalFaces(); - neiGlobal[bFaceI] = - globalNumbering.toGlobal(mesh.nCells()+bFaceI); - faceI++; - } - } - } - syncTools::swapBoundaryFaceList(mesh, neiGlobal, false); - - - // Determine cellCells in global numbering - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - labelListList globalCellCells(mesh.nCells()); - forAll(globalCellCells, cellI) - { - const cell& cFaces = mesh.cells()[cellI]; - - labelList& cCells = globalCellCells[cellI]; - - cCells.setSize(cFaces.size()); - - // Collect neighbouring cells/faces - label nNbr = 0; - forAll(cFaces, i) - { - label faceI = cFaces[i]; - - if (mesh.isInternalFace(faceI)) - { - label nbrCellI = own[faceI]; - if (nbrCellI == cellI) - { - nbrCellI = nei[faceI]; - } - cCells[nNbr++] = globalNumbering.toGlobal(nbrCellI); - } - else - { - label nbrCellI = neiGlobal[faceI-mesh.nInternalFaces()]; - if (nbrCellI != -1) - { - cCells[nNbr++] = nbrCellI; - } - } - } - cCells.setSize(nNbr); - } + const polyBoundaryMesh& patches = mesh_.boundaryMesh(); + const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces(); + const labelList& own = mesh_.faceOwner(); + const labelList& nei = mesh_.faceNeighbour(); // Determine neighbouring global cell Cells // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ labelListList neiGlobalCellCells(nBnd); - for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++) + forAll(patches, patchI) { - neiGlobalCellCells[faceI-mesh.nInternalFaces()] = - globalCellCells[own[faceI]]; + const polyPatch& pp = patches[patchI]; + + if (pp.coupled()) + { + label faceI = pp.start(); + + forAll(pp, i) + { + neiGlobalCellCells[faceI-mesh_.nInternalFaces()] = + globalCellCells[own[faceI]]; + faceI++; + } + } } - syncTools::swapBoundaryFaceList(mesh, neiGlobalCellCells, false); + syncTools::swapBoundaryFaceList(mesh_, neiGlobalCellCells, false); // Construct stencil in global numbering // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - stencil_.setSize(mesh.nFaces()); + faceStencil.setSize(mesh_.nFaces()); - labelHashSet faceStencil; + labelHashSet faceStencilSet; - for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++) + for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) { - faceStencil.clear(); - label globalOwn = globalNumbering.toGlobal(own[faceI]); - faceStencil.insert(globalOwn); + faceStencilSet.clear(); + const labelList& ownCCells = globalCellCells[own[faceI]]; + label globalOwn = ownCCells[0]; + // Insert cellCells forAll(ownCCells, i) { - faceStencil.insert(ownCCells[i]); + faceStencilSet.insert(ownCCells[i]); } - label globalNei = globalNumbering.toGlobal(nei[faceI]); - faceStencil.insert(globalNei); const labelList& neiCCells = globalCellCells[nei[faceI]]; + label globalNei = neiCCells[0]; + // Insert cellCells forAll(neiCCells, i) { - faceStencil.insert(neiCCells[i]); + faceStencilSet.insert(neiCCells[i]); } // Guarantee owner first, neighbour second. - stencil_[faceI].setSize(faceStencil.size()); + faceStencil[faceI].setSize(faceStencilSet.size()); label n = 0; - stencil_[faceI][n++] = globalOwn; - stencil_[faceI][n++] = globalNei; - forAllConstIter(labelHashSet, faceStencil, iter) + faceStencil[faceI][n++] = globalOwn; + faceStencil[faceI][n++] = globalNei; + forAllConstIter(labelHashSet, faceStencilSet, iter) { if (iter.key() != globalOwn && iter.key() != globalNei) { - stencil_[faceI][n++] = iter.key(); + faceStencil[faceI][n++] = iter.key(); } } - //Pout<< "internalface:" << faceI << " toc:" << faceStencil.toc() - // << " stencil:" << stencil_[faceI] << endl; + //Pout<< "internalface:" << faceI << " toc:" << faceStencilSet.toc() + // << " faceStencil:" << faceStencil[faceI] << endl; } forAll(patches, patchI) { @@ -190,41 +119,40 @@ void Foam::extendedStencil::calcFaceStencils { forAll(pp, i) { - faceStencil.clear(); - label globalOwn = globalNumbering.toGlobal(own[faceI]); - faceStencil.insert(globalOwn); + faceStencilSet.clear(); + const labelList& ownCCells = globalCellCells[own[faceI]]; + label globalOwn = ownCCells[0]; forAll(ownCCells, i) { - faceStencil.insert(ownCCells[i]); + faceStencilSet.insert(ownCCells[i]); } - // Get the coupled cell - label globalNei = neiGlobal[faceI-mesh.nInternalFaces()]; - faceStencil.insert(globalNei); + // And the neighbours of the coupled cell const labelList& neiCCells = - neiGlobalCellCells[faceI-mesh.nInternalFaces()]; + neiGlobalCellCells[faceI-mesh_.nInternalFaces()]; + label globalNei = neiCCells[0]; forAll(neiCCells, i) { - faceStencil.insert(neiCCells[i]); + faceStencilSet.insert(neiCCells[i]); } // Guarantee owner first, neighbour second. - stencil_[faceI].setSize(faceStencil.size()); + faceStencil[faceI].setSize(faceStencilSet.size()); label n = 0; - stencil_[faceI][n++] = globalOwn; - stencil_[faceI][n++] = globalNei; - forAllConstIter(labelHashSet, faceStencil, iter) + faceStencil[faceI][n++] = globalOwn; + faceStencil[faceI][n++] = globalNei; + forAllConstIter(labelHashSet, faceStencilSet, iter) { if (iter.key() != globalOwn && iter.key() != globalNei) { - stencil_[faceI][n++] = iter.key(); + faceStencil[faceI][n++] = iter.key(); } } //Pout<< "coupledface:" << faceI - // << " toc:" << faceStencil.toc() - // << " stencil:" << stencil_[faceI] << endl; + // << " toc:" << faceStencilSet.toc() + // << " faceStencil:" << faceStencil[faceI] << endl; faceI++; } @@ -233,31 +161,30 @@ void Foam::extendedStencil::calcFaceStencils { forAll(pp, i) { - faceStencil.clear(); - label globalOwn = globalNumbering.toGlobal(own[faceI]); - faceStencil.insert(globalOwn); + faceStencilSet.clear(); + const labelList& ownCCells = globalCellCells[own[faceI]]; + label globalOwn = ownCCells[0]; forAll(ownCCells, i) { - faceStencil.insert(ownCCells[i]); + faceStencilSet.insert(ownCCells[i]); } - // Guarantee owner first, neighbour second. - stencil_[faceI].setSize(faceStencil.size()); + faceStencil[faceI].setSize(faceStencilSet.size()); label n = 0; - stencil_[faceI][n++] = globalOwn; - forAllConstIter(labelHashSet, faceStencil, iter) + faceStencil[faceI][n++] = globalOwn; + forAllConstIter(labelHashSet, faceStencilSet, iter) { if (iter.key() != globalOwn) { - stencil_[faceI][n++] = iter.key(); + faceStencil[faceI][n++] = iter.key(); } } //Pout<< "boundaryface:" << faceI - // << " toc:" << faceStencil.toc() - // << " stencil:" << stencil_[faceI] << endl; + // << " toc:" << faceStencilSet.toc() + // << " faceStencil:" << faceStencil[faceI] << endl; faceI++; } @@ -266,34 +193,13 @@ void Foam::extendedStencil::calcFaceStencils } -// Calculates extended stencil. This is per face -// - owner -// - cellCells of owner -// - neighbour -// - cellCells of neighbour -// It comes in two parts: -// - a map which collects/distributes all necessary data in a compact array -// - the stencil (a labelList per face) which is a set of indices into this -// compact array. -// The compact array is laid out as follows: -// - first data for current processor (Pstream::myProcNo()) -// - all cells -// - all boundary faces -// - then per processor -// - all used cells and boundary faces -void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh) +Foam::autoPtr Foam::extendedStencil::calcDistributeMap +( + const globalIndex& globalNumbering, + labelListList& faceStencil +) { - const label nBnd = mesh.nFaces()-mesh.nInternalFaces(); - - // Global numbering for cells and boundary faces - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - globalIndex globalNumbering(mesh.nCells()+nBnd); - - - // Calculate stencil in global cell indices - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - calcFaceStencils(mesh, globalNumbering); + const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces(); // Convert stencil to schedule @@ -309,8 +215,8 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh) // these are always all needed. List > globalToProc(Pstream::nProcs()); { - const labelList& procPatchMap = mesh.globalData().procPatchMap(); - const polyBoundaryMesh& patches = mesh.boundaryMesh(); + const labelList& procPatchMap = mesh_.globalData().procPatchMap(); + const polyBoundaryMesh& patches = mesh_.boundaryMesh(); // Presize with (as estimate) size of patch to neighbour. forAll(procPatchMap, procI) @@ -325,9 +231,9 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh) } // Collect all (non-local) globalcells/faces needed. - forAll(stencil_, faceI) + forAll(faceStencil, faceI) { - const labelList& stencilCells = stencil_[faceI]; + const labelList& stencilCells = faceStencil[faceI]; forAll(stencilCells, i) { @@ -359,17 +265,17 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh) } - // forAll(globalToProc, procI) - // { - // Pout<< "From processor:" << procI << " want cells/faces:" << endl; - // forAllConstIter(Map