diff --git a/applications/solvers/combustion/chemFoam/createFields.H b/applications/solvers/combustion/chemFoam/createFields.H index 4a971492c9c38395f0128bc193cef93582da1a2e..b37a8818e311f56d8576e44188aabe7866026f3d 100644 --- a/applications/solvers/combustion/chemFoam/createFields.H +++ b/applications/solvers/combustion/chemFoam/createFields.H @@ -45,9 +45,29 @@ ), thermo.rho() ); + volScalarField& p = thermo.p(); volScalarField& hs = thermo.hs(); + volScalarField Rspecific + ( + IOobject + ( + "Rspecific", + runTime.timeName(), + runTime, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar + ( + "zero", + dimensionSet(dimEnergy/dimMass/dimTemperature), + 0.0 + ) + ); + volVectorField U ( IOobject diff --git a/applications/solvers/combustion/chemFoam/pEqn.H b/applications/solvers/combustion/chemFoam/pEqn.H index 28d240940b6805c6ce5124bd2924e6c36c80a3fc..13f3d603ae7375fe6719e3fd0131f3fdad117270 100644 --- a/applications/solvers/combustion/chemFoam/pEqn.H +++ b/applications/solvers/combustion/chemFoam/pEqn.H @@ -3,7 +3,15 @@ rho = thermo.rho(); if (constProp == "volume") { - p[0] = rho0*R0*thermo.T()[0]; + scalar invW = 0.0; + forAll(Y, i) + { + invW += Y[i][0]/specieData[i].W(); + } + + Rspecific[0] = 1000.0*constant::physicoChemical::R.value()*invW; + + p[0] = rho0*Rspecific[0]*thermo.T()[0]; rho[0] = rho0; } -} \ No newline at end of file +} diff --git a/applications/solvers/combustion/chemFoam/readInitialConditions.H b/applications/solvers/combustion/chemFoam/readInitialConditions.H index 6b26e5fc37e1f8083600d9d01804c08df12b4c47..817ed264f9a18224280ddd54f24becd72472ff51 100644 --- a/applications/solvers/combustion/chemFoam/readInitialConditions.H +++ b/applications/solvers/combustion/chemFoam/readInitialConditions.H @@ -106,6 +106,7 @@ scalar rho0 = rho[0]; scalar u0 = hs0 - p0/rho0; scalar R0 = p0/(rho0*T0); - + Rspecific[0] = R0; + scalar integratedHeat = 0.0; diff --git a/applications/test/PatchEdgeFaceWave/Make/files b/applications/test/PatchEdgeFaceWave/Make/files new file mode 100644 index 0000000000000000000000000000000000000000..01051a355330a7b498c9dd948f10832fba51b1d7 --- /dev/null +++ b/applications/test/PatchEdgeFaceWave/Make/files @@ -0,0 +1,3 @@ +Test-PatchEdgeFaceWave.C + +EXE = $(FOAM_USER_APPBIN)/Test-PatchEdgeFaceWave diff --git a/applications/test/PatchEdgeFaceWave/Make/options b/applications/test/PatchEdgeFaceWave/Make/options new file mode 100644 index 0000000000000000000000000000000000000000..d27c95d033dd5d7b1995c8ff8dc406e35ca1f586 --- /dev/null +++ b/applications/test/PatchEdgeFaceWave/Make/options @@ -0,0 +1,7 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools diff --git a/applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C b/applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C new file mode 100644 index 0000000000000000000000000000000000000000..2dfb19118b2e86503119a8a72872980d6af13942 --- /dev/null +++ b/applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C @@ -0,0 +1,132 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Description + Test PatchEdgeFaceWave. + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "Time.H" +#include "fvMesh.H" +#include "volFields.H" +#include "PatchEdgeFaceWave.H" +#include "patchEdgeFaceInfo.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + argList::validArgs.append("patch"); + +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" + + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + + // Get name of patch + const word patchName = args[1]; + + const polyPatch& patch = patches[patchName]; + + // Data on all edges and faces + List<patchEdgeFaceInfo> allEdgeInfo(patch.nEdges()); + List<patchEdgeFaceInfo> allFaceInfo(patch.size()); + + // Initial seed + DynamicList<label> initialEdges; + DynamicList<patchEdgeFaceInfo> initialEdgesInfo; + + + // Just set an edge on the master + if (Pstream::master()) + { + label edgeI = 0; + Info<< "Starting walk on edge " << edgeI << endl; + + initialEdges.append(edgeI); + const edge& e = patch.edges()[edgeI]; + initialEdgesInfo.append + ( + patchEdgeFaceInfo + ( + e.centre(patch.localPoints()), + 0.0 + ) + ); + } + + + // Walk + PatchEdgeFaceWave + < + primitivePatch, + patchEdgeFaceInfo + > calc + ( + mesh, + patch, + initialEdges, + initialEdgesInfo, + allEdgeInfo, + allFaceInfo, + returnReduce(patch.nEdges(), sumOp<label>()) + ); + + + // Extract as patchField + volScalarField vsf + ( + IOobject + ( + "patchDist", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("patchDist", dimLength, 0.0) + ); + scalarField pf(vsf.boundaryField()[patch.index()].size()); + forAll(pf, faceI) + { + pf[faceI] = Foam::sqrt(allFaceInfo[faceI].distSqr()); + } + vsf.boundaryField()[patch.index()] = pf; + + Info<< "Writing patchDist volScalarField to " << runTime.value() + << endl; + + vsf.write(); + + + Info<< "\nEnd\n" << endl; + return 0; +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files index be2c8d6975524c2caf8d76df74dd2e358bb1757d..28ac7de544f74c5f83b46a863cf0647a8179800a 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files +++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files @@ -1,4 +1,3 @@ -patchPointEdgeCirculator.C createShellMesh.C extrudeToRegionMesh.C diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C index 935b30a5bd1a916fc4d476e07a34ef2c9b65a21f..6bfe6ca62e7c64c6960609011f282335aa4d653c 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C +++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C @@ -31,130 +31,381 @@ License #include "polyAddFace.H" #include "polyModifyFace.H" #include "polyAddCell.H" -#include "patchPointEdgeCirculator.H" +#include "labelPair.H" +#include "indirectPrimitivePatch.H" +#include "mapDistribute.H" +#include "globalMeshData.H" +#include "PatchTools.H" +#include "globalIndex.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // defineTypeNameAndDebug(Foam::createShellMesh, 0); +namespace Foam +{ +template<> +class minEqOp<labelPair> +{ +public: + void operator()(labelPair& x, const labelPair& y) const + { + x[0] = min(x[0], y[0]); + x[1] = min(x[1], y[1]); + } +}; +} + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +// Synchronise edges +void Foam::createShellMesh::syncEdges +( + const globalMeshData& globalData, + + const labelList& patchEdges, + const labelList& coupledEdges, + const PackedBoolList& sameEdgeOrientation, + + PackedBoolList& isChangedEdge, + DynamicList<label>& changedEdges, + labelPairList& allEdgeData +) +{ + const mapDistribute& map = globalData.globalEdgeSlavesMap(); + const PackedBoolList& cppOrientation = globalData.globalEdgeOrientation(); + + // Convert patch-edge data into cpp-edge data + labelPairList cppEdgeData + ( + map.constructSize(), + labelPair(labelMax, labelMax) + ); + + forAll(patchEdges, i) + { + label patchEdgeI = patchEdges[i]; + label coupledEdgeI = coupledEdges[i]; + + if (isChangedEdge[patchEdgeI]) + { + const labelPair& data = allEdgeData[patchEdgeI]; + + // Patch-edge data needs to be converted into coupled-edge data + // (optionally flipped) and consistent in orientation with + // other coupled edge (optionally flipped) + if (sameEdgeOrientation[i] == cppOrientation[coupledEdgeI]) + { + cppEdgeData[coupledEdgeI] = data; + } + else + { + cppEdgeData[coupledEdgeI] = labelPair(data[1], data[0]); + } + } + } + + // Synchronise + globalData.syncData + ( + cppEdgeData, + globalData.globalEdgeSlaves(), + globalData.globalEdgeTransformedSlaves(), + map, + minEqOp<labelPair>() + ); + + // Back from cpp-edge to patch-edge data + forAll(patchEdges, i) + { + label patchEdgeI = patchEdges[i]; + label coupledEdgeI = coupledEdges[i]; + + if (cppEdgeData[coupledEdgeI] != labelPair(labelMax, labelMax)) + { + const labelPair& data = cppEdgeData[coupledEdgeI]; + + if (sameEdgeOrientation[i] == cppOrientation[coupledEdgeI]) + { + allEdgeData[patchEdgeI] = data; + } + else + { + allEdgeData[patchEdgeI] = labelPair(data[1], data[0]); + } + + if (!isChangedEdge[patchEdgeI]) + { + changedEdges.append(patchEdgeI); + isChangedEdge[patchEdgeI] = true; + } + } + } +} + + void Foam::createShellMesh::calcPointRegions ( + const globalMeshData& globalData, const primitiveFacePatch& patch, const PackedBoolList& nonManifoldEdge, - faceList& pointRegions, - labelList& regionPoints + + faceList& pointGlobalRegions, + faceList& pointLocalRegions, + labelList& localToGlobalRegion ) { - pointRegions.setSize(patch.size()); - forAll(pointRegions, faceI) + const indirectPrimitivePatch& cpp = globalData.coupledPatch(); + + // Calculate correspondence between patch and globalData.coupledPatch. + labelList patchEdges; + labelList coupledEdges; + PackedBoolList sameEdgeOrientation; + PatchTools::matchEdges + ( + cpp, + patch, + + coupledEdges, + patchEdges, + sameEdgeOrientation + ); + + + // Initial unique regions + // ~~~~~~~~~~~~~~~~~~~~~~ + // These get merged later on across connected edges. + + // 1. Count + label nMaxRegions = 0; + forAll(patch.localFaces(), faceI) { const face& f = patch.localFaces()[faceI]; - pointRegions[faceI].setSize(f.size(), -1); + nMaxRegions += f.size(); } + const globalIndex globalRegions(nMaxRegions); + + // 2. Assign unique regions label nRegions = 0; - forAll(pointRegions, faceI) + pointGlobalRegions.setSize(patch.size()); + forAll(pointGlobalRegions, faceI) { const face& f = patch.localFaces()[faceI]; + labelList& pRegions = pointGlobalRegions[faceI]; + pRegions.setSize(f.size()); + forAll(pRegions, fp) + { + pRegions[fp] = globalRegions.toGlobal(nRegions++); + } + } - forAll(f, fp) + + DynamicList<label> changedEdges(patch.nEdges()); + labelPairList allEdgeData(patch.nEdges(), labelPair(labelMax, labelMax)); + PackedBoolList isChangedEdge(patch.nEdges()); + + + // Fill initial seed + // ~~~~~~~~~~~~~~~~~ + + forAll(patch.edgeFaces(), edgeI) + { + if (!nonManifoldEdge[edgeI]) + { + // Take over value from one face only. + const edge& e = patch.edges()[edgeI]; + label faceI = patch.edgeFaces()[edgeI][0]; + const face& f = patch.localFaces()[faceI]; + + label fp0 = findIndex(f, e[0]); + label fp1 = findIndex(f, e[1]); + allEdgeData[edgeI] = labelPair + ( + pointGlobalRegions[faceI][fp0], + pointGlobalRegions[faceI][fp1] + ); + if (!isChangedEdge[edgeI]) + { + changedEdges.append(edgeI); + isChangedEdge[edgeI] = true; + } + } + } + + + syncEdges + ( + globalData, + + patchEdges, + coupledEdges, + sameEdgeOrientation, + + isChangedEdge, + changedEdges, + allEdgeData + ); + + + // Edge-Face-Edge walk across patch + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Across edge minimum regions win + + while (true) + { + // From edge to face + // ~~~~~~~~~~~~~~~~~ + + DynamicList<label> changedFaces(patch.size()); + PackedBoolList isChangedFace(patch.size()); + + forAll(changedEdges, changedI) + { + label edgeI = changedEdges[changedI]; + const labelPair& edgeData = allEdgeData[edgeI]; + + const edge& e = patch.edges()[edgeI]; + const labelList& eFaces = patch.edgeFaces()[edgeI]; + + forAll(eFaces, i) + { + label faceI = eFaces[i]; + const face& f = patch.localFaces()[faceI]; + + // Combine edgeData with face data + label fp0 = findIndex(f, e[0]); + if (pointGlobalRegions[faceI][fp0] > edgeData[0]) + { + pointGlobalRegions[faceI][fp0] = edgeData[0]; + if (!isChangedFace[faceI]) + { + isChangedFace[faceI] = true; + changedFaces.append(faceI); + } + } + + label fp1 = findIndex(f, e[1]); + if (pointGlobalRegions[faceI][fp1] > edgeData[1]) + { + pointGlobalRegions[faceI][fp1] = edgeData[1]; + if (!isChangedFace[faceI]) + { + isChangedFace[faceI] = true; + changedFaces.append(faceI); + } + } + } + } + + + label nChangedFaces = returnReduce(changedFaces.size(), sumOp<label>()); + if (nChangedFaces == 0) { - if (pointRegions[faceI][fp] == -1) + break; + } + + + // From face to edge + // ~~~~~~~~~~~~~~~~~ + + isChangedEdge = false; + changedEdges.clear(); + + forAll(changedFaces, i) + { + label faceI = changedFaces[i]; + const face& f = patch.localFaces()[faceI]; + const labelList& fEdges = patch.faceEdges()[faceI]; + + forAll(fEdges, fp) { - // Found unassigned point. Distribute current region. - label pointI = f[fp]; - label edgeI = patch.faceEdges()[faceI][fp]; - - patchPointEdgeCirculator circ - ( - patch, - nonManifoldEdge, - edgeI, - findIndex(patch.edgeFaces()[edgeI], faceI), - pointI - ); - - for - ( - patchPointEdgeCirculator iter = circ.begin(); - iter != circ.end(); - ++iter - ) + label edgeI = fEdges[fp]; + + if (!nonManifoldEdge[edgeI]) { - label face2 = iter.faceID(); + const edge& e = patch.edges()[edgeI]; + label fp0 = findIndex(f, e[0]); + label region0 = pointGlobalRegions[faceI][fp0]; + label fp1 = findIndex(f, e[1]); + label region1 = pointGlobalRegions[faceI][fp1]; - if (face2 != -1) + if + ( + (allEdgeData[edgeI][0] > region0) + || (allEdgeData[edgeI][1] > region1) + ) { - const face& f2 = patch.localFaces()[face2]; - label fp2 = findIndex(f2, pointI); - label& region = pointRegions[face2][fp2]; - if (region != -1) + allEdgeData[edgeI] = labelPair(region0, region1); + if (!isChangedEdge[edgeI]) { - FatalErrorIn - ( - "createShellMesh::calcPointRegions(..)" - ) << "On point " << pointI - << " at:" << patch.localPoints()[pointI] - << " found region:" << region - << abort(FatalError); + changedEdges.append(edgeI); + isChangedEdge[edgeI] = true; } - region = nRegions; } } - - nRegions++; } } - } + syncEdges + ( + globalData, - // From region back to originating point (many to one, a point might - // have multiple regions though) - regionPoints.setSize(nRegions); - forAll(pointRegions, faceI) - { - const face& f = patch.localFaces()[faceI]; + patchEdges, + coupledEdges, + sameEdgeOrientation, - forAll(f, fp) + isChangedEdge, + changedEdges, + allEdgeData + ); + + + label nChangedEdges = returnReduce(changedEdges.size(), sumOp<label>()); + if (nChangedEdges == 0) { - regionPoints[pointRegions[faceI][fp]] = f[fp]; + break; } } - if (debug) + + // Assign local regions + // ~~~~~~~~~~~~~~~~~~~~ + + // Calculate addressing from global region back to local region + pointLocalRegions.setSize(patch.size()); + Map<label> globalToLocalRegion(globalRegions.localSize()/4); + DynamicList<label> dynLocalToGlobalRegion(globalToLocalRegion.size()); + forAll(patch.localFaces(), faceI) { - const labelListList& pointFaces = patch.pointFaces(); - forAll(pointFaces, pointI) + const face& f = patch.localFaces()[faceI]; + face& pRegions = pointLocalRegions[faceI]; + pRegions.setSize(f.size()); + forAll(f, fp) { - label region = -1; - const labelList& pFaces = pointFaces[pointI]; - forAll(pFaces, i) - { - label faceI = pFaces[i]; - const face& f = patch.localFaces()[faceI]; - label fp = findIndex(f, pointI); + label globalRegionI = pointGlobalRegions[faceI][fp]; - if (region == -1) - { - region = pointRegions[faceI][fp]; - } - else if (region != pointRegions[faceI][fp]) - { - Pout<< "Non-manifold point:" << pointI - << " at " << patch.localPoints()[pointI] - << " region:" << region - << " otherRegion:" << pointRegions[faceI][fp] - << endl; + Map<label>::iterator fnd = globalToLocalRegion.find(globalRegionI); - } + if (fnd != globalToLocalRegion.end()) + { + // Already encountered this global region. Assign same local one + pRegions[fp] = fnd(); + } + else + { + // Region not yet seen. Create new one + label localRegionI = globalToLocalRegion.size(); + pRegions[fp] = localRegionI; + globalToLocalRegion.insert(globalRegionI, localRegionI); + dynLocalToGlobalRegion.append(globalRegionI); } } } + localToGlobalRegion.transfer(dynLocalToGlobalRegion); } diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.H b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.H index 627191c6a6eed3bc8e93f3e8b092978b03cef7f0..7517a956df1dc16d241bb26cdc221a8c0281371f 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.H +++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.H @@ -41,6 +41,7 @@ SourceFiles #include "primitiveFacePatch.H" #include "PackedBoolList.H" +#include "labelPair.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -50,6 +51,7 @@ namespace Foam // Forward declaration of classes class mapPolyMesh; class polyTopoChange; +class globalMeshData; /*---------------------------------------------------------------------------*\ Class createShellMesh Declaration @@ -80,6 +82,18 @@ class createShellMesh // Private Member Functions + static void syncEdges + ( + const globalMeshData&, + const labelList&, + const labelList&, + const PackedBoolList& sameEdgeOrientation, + PackedBoolList& isChangedEdge, + DynamicList<label>& changedEdges, + labelPairList& allEdgeData + ); + + //- Disallow default bitwise copy construct createShellMesh(const createShellMesh&); @@ -142,13 +156,22 @@ public: // Other - //- Helper: calculate point regions + //- Helper: calculate point regions. The point region is the + // same on all faces connected to a point if they can be + // reached through a face-edge-face walk without crossing + // the nonManifoldEdge. + // pointGlobalRegions : non-compact. Guaranteed to be the same + // across processors. + // pointLocalRegions : compact. + // localToGlobalRegion : local to global region. static void calcPointRegions ( + const globalMeshData& globalData, const primitiveFacePatch& patch, const PackedBoolList& nonManifoldEdge, - faceList& pointRegions, - labelList& regionPoints + faceList& pointGlobalRegions, + faceList& pointLocalRegions, + labelList& localToGlobalRegion ); //- Play commands into polyTopoChange to create layer mesh. diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C index 0a085485f14206e5e8d113a1b54a7e1cc2efceb5..40e4021284244e8910b597bd62ea23abf830047d 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C +++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C @@ -120,58 +120,93 @@ Usage \*---------------------------------------------------------------------------*/ #include "argList.H" -#include "Time.H" #include "fvMesh.H" #include "polyTopoChange.H" -#include "patchPointEdgeCirculator.H" #include "OFstream.H" #include "meshTools.H" #include "mappedWallPolyPatch.H" #include "createShellMesh.H" -#include "volFields.H" -#include "surfaceFields.H" #include "syncTools.H" #include "cyclicPolyPatch.H" #include "wedgePolyPatch.H" #include "nonuniformTransformCyclicPolyPatch.H" #include "extrudeModel.H" +#include "globalIndex.H" +#include "addPatchCellLayer.H" -using namespace Foam; +#include "volFields.H" +#include "surfaceFields.H" +#include "pointFields.H" +//#include "ReadFields.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -template<class GeoField> -void addPatchFields(fvMesh& mesh, const word& patchFieldType) -{ - HashTable<const GeoField*> flds - ( - mesh.objectRegistry::lookupClass<GeoField>() - ); - - forAllConstIter(typename HashTable<const GeoField*>, flds, iter) - { - const GeoField& fld = *iter(); +using namespace Foam; - typename GeoField::GeometricBoundaryField& bfld = - const_cast<typename GeoField::GeometricBoundaryField&> - ( - fld.boundaryField() - ); +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - label sz = bfld.size(); - bfld.setSize(sz+1); - bfld.set - ( - sz, - GeoField::PatchFieldType::New - ( - patchFieldType, - mesh.boundary()[sz], - fld.dimensionedInternalField() - ) - ); - } -} +//template<class GeoField> +//void addPatchFields(const fvMesh& mesh, const word& patchFieldType) +//{ +// HashTable<const GeoField*> flds +// ( +// mesh.objectRegistry::lookupClass<GeoField>() +// ); +// +// forAllConstIter(typename HashTable<const GeoField*>, flds, iter) +// { +// const GeoField& fld = *iter(); +// +// typename GeoField::GeometricBoundaryField& bfld = +// const_cast<typename GeoField::GeometricBoundaryField&> +// ( +// fld.boundaryField() +// ); +// +// label sz = bfld.size(); +// +// for (label i = 0; i < sz; i++) +// { +// bfld.set +// ( +// i, +// bfld.clone(GeoField::PatchFieldType::New +// ( +// patchFieldType, +// fld.mesh().boundary()[sz], +// fld.dimensionedInternalField() +// ) +// ); +// +// +// +// Pout<< "fld:" << fld.name() << " had " << sz << " patches." << endl; +// Pout<< "fld before:" << fld << endl; +// Pout<< "adding on patch:" << fld.mesh().boundary()[sz].name() << endl; +// +// bfld.setSize(sz+1); +// bfld.set +// ( +// sz, +// GeoField::PatchFieldType::New +// ( +// patchFieldType, +// fld.mesh().boundary()[sz], +// fld.dimensionedInternalField() +// ) +// ); +// +// bfld[sz].operator=(pTraits<typename GeoField::value_type>::zero); +// +// Pout<< "fld:" << fld.name() << " now " << bfld.size() << " patches." +// << endl; +// +// const typename GeoField::PatchFieldType& pfld = bfld[sz]; +// Pout<< "pfld:" << pfld << endl; +// +// +// Pout<< "fld value:" << fld << endl; +// } +//} // Remove last patch field @@ -219,264 +254,467 @@ void reorderPatchFields(fvMesh& mesh, const labelList& oldToNew) } -void addAllPatchFields(fvMesh& mesh, const label insertPatchI) +//void addCalculatedPatchFields(const fvMesh& mesh) +//{ +// addPatchFields<volScalarField> +// ( +// mesh, +// calculatedFvPatchField<scalar>::typeName +// ); +// addPatchFields<volVectorField> +// ( +// mesh, +// calculatedFvPatchField<vector>::typeName +// ); +// addPatchFields<volSphericalTensorField> +// ( +// mesh, +// calculatedFvPatchField<sphericalTensor>::typeName +// ); +// addPatchFields<volSymmTensorField> +// ( +// mesh, +// calculatedFvPatchField<symmTensor>::typeName +// ); +// addPatchFields<volTensorField> +// ( +// mesh, +// calculatedFvPatchField<tensor>::typeName +// ); +// +// // Surface fields +// +// addPatchFields<surfaceScalarField> +// ( +// mesh, +// calculatedFvPatchField<scalar>::typeName +// ); +// addPatchFields<surfaceVectorField> +// ( +// mesh, +// calculatedFvPatchField<vector>::typeName +// ); +// addPatchFields<surfaceSphericalTensorField> +// ( +// mesh, +// calculatedFvPatchField<sphericalTensor>::typeName +// ); +// addPatchFields<surfaceSymmTensorField> +// ( +// mesh, +// calculatedFvPatchField<symmTensor>::typeName +// ); +// addPatchFields<surfaceTensorField> +// ( +// mesh, +// calculatedFvPatchField<tensor>::typeName +// ); +// +// // Point fields +// +// addPatchFields<pointScalarField> +// ( +// mesh, +// calculatedFvPatchField<scalar>::typeName +// ); +// addPatchFields<pointVectorField> +// ( +// mesh, +// calculatedFvPatchField<vector>::typeName +// ); +//} +// +// +//void addAllPatchFields(fvMesh& mesh, const label insertPatchI) +//{ +// polyBoundaryMesh& polyPatches = +// const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()); +// fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary()); +// +// label sz = polyPatches.size(); +// +// addPatchFields<volScalarField> +// ( +// mesh, +// calculatedFvPatchField<scalar>::typeName +// ); +// addPatchFields<volVectorField> +// ( +// mesh, +// calculatedFvPatchField<vector>::typeName +// ); +// addPatchFields<volSphericalTensorField> +// ( +// mesh, +// calculatedFvPatchField<sphericalTensor>::typeName +// ); +// addPatchFields<volSymmTensorField> +// ( +// mesh, +// calculatedFvPatchField<symmTensor>::typeName +// ); +// addPatchFields<volTensorField> +// ( +// mesh, +// calculatedFvPatchField<tensor>::typeName +// ); +// +// // Surface fields +// +// addPatchFields<surfaceScalarField> +// ( +// mesh, +// calculatedFvPatchField<scalar>::typeName +// ); +// addPatchFields<surfaceVectorField> +// ( +// mesh, +// calculatedFvPatchField<vector>::typeName +// ); +// addPatchFields<surfaceSphericalTensorField> +// ( +// mesh, +// calculatedFvPatchField<sphericalTensor>::typeName +// ); +// addPatchFields<surfaceSymmTensorField> +// ( +// mesh, +// calculatedFvPatchField<symmTensor>::typeName +// ); +// addPatchFields<surfaceTensorField> +// ( +// mesh, +// calculatedFvPatchField<tensor>::typeName +// ); +// +// // Create reordering list +// // patches before insert position stay as is +// labelList oldToNew(sz); +// for (label i = 0; i < insertPatchI; i++) +// { +// oldToNew[i] = i; +// } +// // patches after insert position move one up +// for (label i = insertPatchI; i < sz-1; i++) +// { +// oldToNew[i] = i+1; +// } +// // appended patch gets moved to insert position +// oldToNew[sz-1] = insertPatchI; +// +// // Shuffle into place +// polyPatches.reorder(oldToNew); +// fvPatches.reorder(oldToNew); +// +// reorderPatchFields<volScalarField>(mesh, oldToNew); +// reorderPatchFields<volVectorField>(mesh, oldToNew); +// reorderPatchFields<volSphericalTensorField>(mesh, oldToNew); +// reorderPatchFields<volSymmTensorField>(mesh, oldToNew); +// reorderPatchFields<volTensorField>(mesh, oldToNew); +// reorderPatchFields<surfaceScalarField>(mesh, oldToNew); +// reorderPatchFields<surfaceVectorField>(mesh, oldToNew); +// reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew); +// reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew); +// reorderPatchFields<surfaceTensorField>(mesh, oldToNew); +//} + + +//// Adds patch if not yet there. Returns patchID. +//template<class PatchType> +//label addPatch(fvMesh& mesh, const word& patchName, const dictionary& dict) +//{ +// polyBoundaryMesh& polyPatches = +// const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()); +// +// label patchI = polyPatches.findPatchID(patchName); +// if (patchI != -1) +// { +// if (isA<PatchType>(polyPatches[patchI])) +// { +// // Already there +// return patchI; +// } +// else +// { +// FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)") +// << "Already have patch " << patchName +// << " but of type " << PatchType::typeName +// << exit(FatalError); +// } +// } +// +// +// label insertPatchI = polyPatches.size(); +// label startFaceI = mesh.nFaces(); +// +// forAll(polyPatches, patchI) +// { +// const polyPatch& pp = polyPatches[patchI]; +// +// if (isA<processorPolyPatch>(pp)) +// { +// insertPatchI = patchI; +// startFaceI = pp.start(); +// break; +// } +// } +// +// dictionary patchDict(dict); +// patchDict.set("type", PatchType::typeName); +// patchDict.set("nFaces", 0); +// patchDict.set("startFace", startFaceI); +// +// +// // Below is all quite a hack. Feel free to change once there is a better +// // mechanism to insert and reorder patches. +// +// // Clear local fields and e.g. polyMesh parallelInfo. +// mesh.clearOut(); +// +// label sz = polyPatches.size(); +// +// fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary()); +// +// // Add polyPatch at the end +// polyPatches.setSize(sz+1); +// polyPatches.set +// ( +// sz, +// polyPatch::New +// ( +// patchName, +// patchDict, +// insertPatchI, +// polyPatches +// ) +// ); +// fvPatches.setSize(sz+1); +// fvPatches.set +// ( +// sz, +// fvPatch::New +// ( +// polyPatches[sz], // point to newly added polyPatch +// mesh.boundary() +// ) +// ); +// +// addAllPatchFields(mesh, insertPatchI); +// +// return insertPatchI; +//} +// +// +//template<class PatchType> +//label addPatch(fvMesh& mesh, const word& patchName) +//{ +//Pout<< "addPatch:" << patchName << endl; +// +// polyBoundaryMesh& polyPatches = +// const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()); +// +// label patchI = polyPatches.findPatchID(patchName); +// if (patchI != -1) +// { +// if (isA<PatchType>(polyPatches[patchI])) +// { +// // Already there +// return patchI; +// } +// else +// { +// FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)") +// << "Already have patch " << patchName +// << " but of type " << PatchType::typeName +// << exit(FatalError); +// } +// } +// +// +// label insertPatchI = polyPatches.size(); +// label startFaceI = mesh.nFaces(); +// +// forAll(polyPatches, patchI) +// { +// const polyPatch& pp = polyPatches[patchI]; +// +// if (isA<processorPolyPatch>(pp)) +// { +// insertPatchI = patchI; +// startFaceI = pp.start(); +// break; +// } +// } +// +// // Below is all quite a hack. Feel free to change once there is a better +// // mechanism to insert and reorder patches. +// +// // Clear local fields and e.g. polyMesh parallelInfo. +// mesh.clearOut(); +// +// label sz = polyPatches.size(); +// +// fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary()); +// +// // Add polyPatch at the end +// polyPatches.setSize(sz+1); +// polyPatches.set +// ( +// sz, +// polyPatch::New +// ( +// PatchType::typeName, +// patchName, +// 0, // size +// startFaceI, +// insertPatchI, +// polyPatches +// ) +// ); +// fvPatches.setSize(sz+1); +// fvPatches.set +// ( +// sz, +// fvPatch::New +// ( +// polyPatches[sz], // point to newly added polyPatch +// mesh.boundary() +// ) +// ); +// +// addAllPatchFields(mesh, insertPatchI); +// +// return insertPatchI; +//} + + +label findPatchID(const List<polyPatch*>& newPatches, const word& name) { - polyBoundaryMesh& polyPatches = - const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()); - fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary()); - - label sz = polyPatches.size(); - - addPatchFields<volScalarField> - ( - mesh, - calculatedFvPatchField<scalar>::typeName - ); - addPatchFields<volVectorField> - ( - mesh, - calculatedFvPatchField<vector>::typeName - ); - addPatchFields<volSphericalTensorField> - ( - mesh, - calculatedFvPatchField<sphericalTensor>::typeName - ); - addPatchFields<volSymmTensorField> - ( - mesh, - calculatedFvPatchField<symmTensor>::typeName - ); - addPatchFields<volTensorField> - ( - mesh, - calculatedFvPatchField<tensor>::typeName - ); - - // Surface fields - - addPatchFields<surfaceScalarField> - ( - mesh, - calculatedFvPatchField<scalar>::typeName - ); - addPatchFields<surfaceVectorField> - ( - mesh, - calculatedFvPatchField<vector>::typeName - ); - addPatchFields<surfaceSphericalTensorField> - ( - mesh, - calculatedFvPatchField<sphericalTensor>::typeName - ); - addPatchFields<surfaceSymmTensorField> - ( - mesh, - calculatedFvPatchField<symmTensor>::typeName - ); - addPatchFields<surfaceTensorField> - ( - mesh, - calculatedFvPatchField<tensor>::typeName - ); - - // Create reordering list - // patches before insert position stay as is - labelList oldToNew(sz); - for (label i = 0; i < insertPatchI; i++) - { - oldToNew[i] = i; - } - // patches after insert position move one up - for (label i = insertPatchI; i < sz-1; i++) + forAll(newPatches, i) { - oldToNew[i] = i+1; + if (newPatches[i]->name() == name) + { + return i; + } } - // appended patch gets moved to insert position - oldToNew[sz-1] = insertPatchI; - - // Shuffle into place - polyPatches.reorder(oldToNew); - fvPatches.reorder(oldToNew); - - reorderPatchFields<volScalarField>(mesh, oldToNew); - reorderPatchFields<volVectorField>(mesh, oldToNew); - reorderPatchFields<volSphericalTensorField>(mesh, oldToNew); - reorderPatchFields<volSymmTensorField>(mesh, oldToNew); - reorderPatchFields<volTensorField>(mesh, oldToNew); - reorderPatchFields<surfaceScalarField>(mesh, oldToNew); - reorderPatchFields<surfaceVectorField>(mesh, oldToNew); - reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew); - reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew); - reorderPatchFields<surfaceTensorField>(mesh, oldToNew); + return -1; } -// Adds patch if not yet there. Returns patchID. template<class PatchType> -label addPatch(fvMesh& mesh, const word& patchName, const dictionary& dict) +label addPatch +( + const polyBoundaryMesh& patches, + const word& patchName, + DynamicList<polyPatch*>& newPatches +) { - polyBoundaryMesh& polyPatches = - const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()); + label patchI = findPatchID(newPatches, patchName); - label patchI = polyPatches.findPatchID(patchName); if (patchI != -1) { - if (isA<PatchType>(polyPatches[patchI])) + if (isA<PatchType>(*newPatches[patchI])) { // Already there return patchI; } else { - FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)") - << "Already have patch " << patchName - << " but of type " << PatchType::typeName + FatalErrorIn + ( + "addPatch<PatchType>(const polyBoundaryMesh&," + " const word&, DynamicList<polyPatch*>)" + ) << "Already have patch " << patchName + << " but of type " << newPatches[patchI]->type() << exit(FatalError); } } - label insertPatchI = polyPatches.size(); - label startFaceI = mesh.nFaces(); + patchI = newPatches.size(); - forAll(polyPatches, patchI) + label startFaceI = 0; + if (patchI > 0) { - const polyPatch& pp = polyPatches[patchI]; - - if (isA<processorPolyPatch>(pp)) - { - insertPatchI = patchI; - startFaceI = pp.start(); - break; - } + const polyPatch& pp = *newPatches.last(); + startFaceI = pp.start()+pp.size(); } - dictionary patchDict(dict); - patchDict.set("type", PatchType::typeName); - patchDict.set("nFaces", 0); - patchDict.set("startFace", startFaceI); - - // Below is all quite a hack. Feel free to change once there is a better - // mechanism to insert and reorder patches. - - // Clear local fields and e.g. polyMesh parallelInfo. - mesh.clearOut(); - - label sz = polyPatches.size(); - - fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary()); - - // Add polyPatch at the end - polyPatches.setSize(sz+1); - polyPatches.set + newPatches.append ( - sz, polyPatch::New ( + PatchType::typeName, patchName, - patchDict, - insertPatchI, - polyPatches - ) + 0, // size + startFaceI, // nFaces + patchI, + patches + ).ptr() ); - fvPatches.setSize(sz+1); - fvPatches.set - ( - sz, - fvPatch::New - ( - polyPatches[sz], // point to newly added polyPatch - mesh.boundary() - ) - ); - - addAllPatchFields(mesh, insertPatchI); - return insertPatchI; + return patchI; } template<class PatchType> -label addPatch(fvMesh& mesh, const word& patchName) +label addPatch +( + const polyBoundaryMesh& patches, + const word& patchName, + const dictionary& dict, + DynamicList<polyPatch*>& newPatches +) { - polyBoundaryMesh& polyPatches = - const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()); + label patchI = findPatchID(newPatches, patchName); - label patchI = polyPatches.findPatchID(patchName); if (patchI != -1) { - if (isA<PatchType>(polyPatches[patchI])) + if (isA<PatchType>(*newPatches[patchI])) { // Already there return patchI; } else { - FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)") - << "Already have patch " << patchName - << " but of type " << PatchType::typeName + FatalErrorIn + ( + "addPatch<PatchType>(const polyBoundaryMesh&," + " const word&, DynamicList<polyPatch*>)" + ) << "Already have patch " << patchName + << " but of type " << newPatches[patchI]->type() << exit(FatalError); } } - label insertPatchI = polyPatches.size(); - label startFaceI = mesh.nFaces(); + patchI = newPatches.size(); - forAll(polyPatches, patchI) + label startFaceI = 0; + if (patchI > 0) { - const polyPatch& pp = polyPatches[patchI]; - - if (isA<processorPolyPatch>(pp)) - { - insertPatchI = patchI; - startFaceI = pp.start(); - break; - } + const polyPatch& pp = *newPatches.last(); + startFaceI = pp.start()+pp.size(); } - // Below is all quite a hack. Feel free to change once there is a better - // mechanism to insert and reorder patches. - - // Clear local fields and e.g. polyMesh parallelInfo. - mesh.clearOut(); - - label sz = polyPatches.size(); - - fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary()); + dictionary patchDict(dict); + patchDict.set("type", PatchType::typeName); + patchDict.set("nFaces", 0); + patchDict.set("startFace", startFaceI); - // Add polyPatch at the end - polyPatches.setSize(sz+1); - polyPatches.set + newPatches.append ( - sz, polyPatch::New ( - PatchType::typeName, patchName, - 0, // size - startFaceI, - insertPatchI, - polyPatches - ) - ); - fvPatches.setSize(sz+1); - fvPatches.set - ( - sz, - fvPatch::New - ( - polyPatches[sz], // point to newly added polyPatch - mesh.boundary() - ) + patchDict, + patchI, + patches + ).ptr() ); - addAllPatchFields(mesh, insertPatchI); - - return insertPatchI; + return patchI; } @@ -506,6 +744,8 @@ void reorderPatches reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew); reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew); reorderPatchFields<surfaceTensorField>(mesh, oldToNew); + reorderPatchFields<pointScalarField>(mesh, oldToNew); + reorderPatchFields<pointVectorField>(mesh, oldToNew); // Remove last. polyPatches.setSize(nNewPatches); @@ -520,6 +760,8 @@ void reorderPatches trimPatchFields<surfaceSphericalTensorField>(mesh, nNewPatches); trimPatchFields<surfaceSymmTensorField>(mesh, nNewPatches); trimPatchFields<surfaceTensorField>(mesh, nNewPatches); + trimPatchFields<pointScalarField>(mesh, nNewPatches); + trimPatchFields<pointVectorField>(mesh, nNewPatches); } @@ -531,15 +773,53 @@ void deleteEmptyPatches(fvMesh& mesh) labelList oldToNew(patches.size()); label usedI = 0; label notUsedI = patches.size(); + + //Pout<< "deleteEmptyPatches:" << endl; + //forAll(patches, patchI) + //{ + // Pout<< " patch:" << patchI << " name:" << patches[patchI].name() + // << " start:" << patches[patchI].start() + // << " nFaces:" << patches[patchI].size() + // << " index:" << patches[patchI].index() + // << endl; + //} + //Pout<< endl; + + + // Add all the non-empty, non-processor patches forAll(patches, patchI) { - if (returnReduce(patches[patchI].size(), sumOp<label>()) == 0) + if (isA<processorPolyPatch>(patches[patchI])) { - oldToNew[patchI] = --notUsedI; + // Processor patches are unique per processor so look at local + // size only + if (patches[patchI].size() == 0) + { + Pout<< "Deleting processor patch " << patchI + << " name:" << patches[patchI].name() + << endl; + oldToNew[patchI] = --notUsedI; + } + else + { + oldToNew[patchI] = usedI++; + } } else { - oldToNew[patchI] = usedI++; + // All non-processor patches are present everywhere to reduce + // size + if (returnReduce(patches[patchI].size(), sumOp<label>()) == 0) + { + Pout<< "Deleting patch " << patchI + << " name:" << patches[patchI].name() + << endl; + oldToNew[patchI] = --notUsedI; + } + else + { + oldToNew[patchI] = usedI++; + } } } @@ -628,7 +908,11 @@ label findUncoveredPatchFace } -// Count the number of faces in patches that need to be created +// Count the number of faces in patches that need to be created. Calculates: +// zoneSidePatch[zoneI] : the number of side faces to be created +// zoneZonePatch[zoneA,zoneB] : the number of faces inbetween zoneA and B +// Since this only counts we're not taking the processor patches into +// account. void countExtrudePatches ( const fvMesh& mesh, @@ -691,200 +975,410 @@ void countExtrudePatches } -// Constrain&sync normals on points that are on coupled patches to make sure -// the face extruded from the edge has a valid normal with its coupled -// equivalent. -// Note that only points on cyclic edges need to be constrained and not -// all points touching cyclics since only edges become faces. -void constrainCoupledNormals +void addCouplingPatches +( + const fvMesh& mesh, + const word& regionName, + const word& shellRegionName, + const wordList& zoneNames, + const wordList& zoneShadowNames, + const boolList& isInternal, + DynamicList<polyPatch*>& newPatches, + labelList& interRegionTopPatch, + labelList& interRegionBottomPatch +) +{ + Pout<< "Adding coupling patches:" << nl << nl + << "patchID\tpatch\ttype" << nl + << "-------\t-----\t----" + << endl; + + interRegionTopPatch.setSize(zoneNames.size()); + interRegionBottomPatch.setSize(zoneNames.size()); + + label nCoupled = 0; + forAll(zoneNames, i) + { + word interName(regionName+"_to_"+shellRegionName+'_'+zoneNames[i]); + + if (isInternal[i]) + { + interRegionTopPatch[i] = addPatch<mappedWallPolyPatch> + ( + mesh.boundaryMesh(), + interName + "_top", + newPatches + ); + nCoupled++; + Pout<< interRegionTopPatch[i] + << '\t' << newPatches[interRegionTopPatch[i]]->name() + << '\t' << newPatches[interRegionTopPatch[i]]->type() + << nl; + + interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch> + ( + mesh.boundaryMesh(), + interName + "_bottom", + newPatches + ); + nCoupled++; + Pout<< interRegionBottomPatch[i] + << '\t' << newPatches[interRegionBottomPatch[i]]->name() + << '\t' << newPatches[interRegionBottomPatch[i]]->type() + << nl; + } + else if (zoneShadowNames.size() == 0) + { + interRegionTopPatch[i] = addPatch<polyPatch> + ( + mesh.boundaryMesh(), + zoneNames[i] + "_top", + newPatches + ); + nCoupled++; + Pout<< interRegionTopPatch[i] + << '\t' << newPatches[interRegionTopPatch[i]]->name() + << '\t' << newPatches[interRegionTopPatch[i]]->type() + << nl; + + interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch> + ( + mesh.boundaryMesh(), + interName, + newPatches + ); + nCoupled++; + Pout<< interRegionBottomPatch[i] + << '\t' << newPatches[interRegionBottomPatch[i]]->name() + << '\t' << newPatches[interRegionBottomPatch[i]]->type() + << nl; + } + else //patch using shadow face zones. + { + interRegionTopPatch[i] = addPatch<mappedWallPolyPatch> + ( + mesh.boundaryMesh(), + zoneShadowNames[i] + "_top", + newPatches + ); + nCoupled++; + Pout<< interRegionTopPatch[i] + << '\t' << newPatches[interRegionTopPatch[i]]->name() + << '\t' << newPatches[interRegionTopPatch[i]]->type() + << nl; + + interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch> + ( + mesh.boundaryMesh(), + interName, + newPatches + ); + nCoupled++; + Pout<< interRegionBottomPatch[i] + << '\t' << newPatches[interRegionBottomPatch[i]]->name() + << '\t' << newPatches[interRegionBottomPatch[i]]->type() + << nl; + } + } + Pout<< "Added " << nCoupled << " inter-region patches." << nl + << endl; +} + + +void addProcPatches ( const fvMesh& mesh, const primitiveFacePatch& extrudePatch, - const labelList& meshEdges, - const faceList& pointRegions, // per face, per index the region + const labelList& extrudeMeshFaces, - vectorField& regionNormals + labelList& sidePatchID, + DynamicList<polyPatch*>& newPatches ) { - const polyBoundaryMesh& patches = mesh.boundaryMesh(); + // Global face indices engine + const globalIndex globalFaces(mesh.nFaces()); - // Mark edges that are on boundary of extrusion. - Map<label> meshToExtrudEdge + indirectPrimitivePatch pp ( - 2*(extrudePatch.nEdges()-extrudePatch.nInternalEdges()) + IndirectList<face>(mesh.faces(), extrudeMeshFaces), + mesh.points() ); - for - ( - label extrudeEdgeI = extrudePatch.nInternalEdges(); - extrudeEdgeI < extrudePatch.nEdges(); - extrudeEdgeI++ - ) + + forAll(extrudePatch.edges(), edgeI) { - meshToExtrudEdge.insert(meshEdges[extrudeEdgeI], extrudeEdgeI); + const edge& extrudeEdge = extrudePatch.edges()[edgeI]; + const edge& ppEdge = pp.edges()[edgeI]; + + if (extrudeEdge != ppEdge) + { + FatalErrorIn("addProcPatches()") + << "Problem: extrudeEdge:" << extrudeEdge + << " ppEdge:" << ppEdge + << exit(FatalError); + } } - // For owner: normal at first point of edge when walking through faces - // in order. - vectorField edgeNormals0(mesh.nEdges(), vector::zero); - vectorField edgeNormals1(mesh.nEdges(), vector::zero); + // Determine extrudePatch.edgeFaces in global numbering (so across + // coupled patches). + labelListList edgeGlobalFaces + ( + addPatchCellLayer::globalEdgeFaces + ( + mesh, + globalFaces, + pp + ) + ); + + // Calculate for every edge the patch to use. This will be an existing + // patch for uncoupled edge and a possibly new patch (in patchToNbrProc) + // for processor patches. + label nNewPatches; + Map<label> nbrProcToPatch; + Map<label> patchToNbrProc; + + addPatchCellLayer::calcSidePatch + ( + mesh, + globalFaces, + edgeGlobalFaces, + pp, + + sidePatchID, + nNewPatches, + nbrProcToPatch, + patchToNbrProc + ); - // Loop through all edges of patch. If they are to be extruded mark the - // point normals in order. - forAll(patches, patchI) - { - const polyPatch& pp = patches[patchI]; - if (isA<cyclicPolyPatch>(pp)) + // All patchIDs calcSidePatch are in mesh.boundaryMesh() numbering. + // Redo in newPatches numbering. + + Pout<< "Adding inter-processor patches:" << nl << nl + << "patchID\tpatch" << nl + << "-------\t-----" + << endl; + + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + forAll(sidePatchID, edgeI) + { + label meshPatchI = sidePatchID[edgeI]; + if (meshPatchI != -1) { - bool isOwner = refCast<const cyclicPolyPatch>(pp).owner(); + label newPatchI = -1; + if (meshPatchI < patches.size()) + { + // Existing mesh patch. See if I already have it. + newPatchI = findPatchID + ( + newPatches, + patches[meshPatchI].name() + ); + } - forAll(pp.faceEdges(), faceI) + if (newPatchI == -1) { - const labelList& fEdges = pp.faceEdges()[faceI]; - forAll(fEdges, fp) - { - label meshEdgeI = pp.meshEdges()[fEdges[fp]]; - if (meshToExtrudEdge.found(meshEdgeI)) - { - // Edge corresponds to a extrusion edge. Store extrusion - // normals on edge so we can syncTools it. - - //const edge& ppE = pp.edges()[fEdges[fp]]; - //Pout<< "ppedge:" << pp.localPoints()[ppE[0]] - // << pp.localPoints()[ppE[1]] - // << endl; - - const face& f = pp.localFaces()[faceI]; - label fp0 = fp; - label fp1 = f.fcIndex(fp0); - label mp0 = pp[faceI][fp0]; - label mp1 = pp[faceI][fp1]; - - // Find corresponding face and indices. - vector regionN0; - vector regionN1; - { - label exEdgeI = meshToExtrudEdge[meshEdgeI]; - const labelList& eFaces = - extrudePatch.edgeFaces()[exEdgeI]; - // Use 0th face. - label exFaceI = eFaces[0]; - const face& exF = extrudePatch[exFaceI]; - const face& exRegions = pointRegions[exFaceI]; - // Find points - label r0 = exRegions[findIndex(exF, mp0)]; - regionN0 = regionNormals[r0]; - label r1 = exRegions[findIndex(exF, mp1)]; - regionN1 = regionNormals[r1]; - } - - vector& nA = - ( - isOwner - ? edgeNormals0[meshEdgeI] - : edgeNormals1[meshEdgeI] - ); + // Newly added processor patch + label nbrProcI = patchToNbrProc[meshPatchI]; + word name = + "procBoundary" + + Foam::name(Pstream::myProcNo()) + + "to" + + Foam::name(nbrProcI); + + dictionary patchDict; + patchDict.add("myProcNo", Pstream::myProcNo()); + patchDict.add("neighbProcNo", nbrProcI); + + newPatchI = addPatch<processorPolyPatch> + ( + mesh.boundaryMesh(), + name, + patchDict, + newPatches + ); - nA = regionN0; - const vector& cyc0 = pp.pointNormals()[f[fp0]]; - nA -= (nA&cyc0)*cyc0; + Pout<< newPatchI << '\t' << name + << nl; + } - vector& nB = - ( - isOwner - ? edgeNormals1[meshEdgeI] - : edgeNormals0[meshEdgeI] - ); + sidePatchID[edgeI] = newPatchI; + } + } - nB = regionN1; - const vector& cyc1 = pp.pointNormals()[f[fp1]]; - nB -= (nB&cyc1)*cyc1; - } - } - } + + // Clear out sidePatchID for uncoupled edges. Just so we don't have + // to expose all the globalEdgeFaces info. + forAll(sidePatchID, edgeI) + { + if + ( + edgeGlobalFaces[edgeI].size() == 2 + && pp.edgeFaces()[edgeI].size() == 1 + ) + {} + else + { + sidePatchID[edgeI] = -1; } } +} - // Synchronise regionNormals - // ~~~~~~~~~~~~~~~~~~~~~~~~~ +void addZoneSidePatches +( + const fvMesh& mesh, + const word& oneDPolyPatchType, + const wordList& zoneNames, - // Synchronise - syncTools::syncEdgeList - ( - mesh, - edgeNormals0, - maxMagSqrEqOp<vector>(), - vector::zero // nullValue - ); - syncTools::syncEdgeList - ( - mesh, - edgeNormals1, - maxMagSqrEqOp<vector>(), - vector::zero // nullValue - ); + DynamicList<polyPatch*>& newPatches, + labelList& zoneSidePatch +) +{ + Pout<< "Adding patches for sides on zones:" << nl << nl + << "patchID\tpatch" << nl + << "-------\t-----" + << endl; + label nSide = 0; - // Re-work back into regionNormals - forAll(patches, patchI) + forAll(zoneNames, zoneI) { - const polyPatch& pp = patches[patchI]; + if (oneDPolyPatchType != word::null) + { + // Reuse single empty patch. + word patchName; + if (oneDPolyPatchType == "emptyPolyPatch") + { + patchName = "oneDEmptyPatch"; + zoneSidePatch[zoneI] = addPatch<emptyPolyPatch> + ( + mesh.boundaryMesh(), + patchName, + newPatches + ); + } + else if (oneDPolyPatchType == "wedgePolyPatch") + { + patchName = "oneDWedgePatch"; + zoneSidePatch[zoneI] = addPatch<wedgePolyPatch> + ( + mesh.boundaryMesh(), + patchName, + newPatches + ); + } + else + { + FatalErrorIn("addZoneSidePatches()") + << "Type " << oneDPolyPatchType << " does not exist " + << exit(FatalError); + } - if (isA<cyclicPolyPatch>(pp)) + Pout<< zoneSidePatch[zoneI] << '\t' << patchName << nl; + + nSide++; + } + else if (zoneSidePatch[zoneI] > 0) { - bool isOwner = refCast<const cyclicPolyPatch>(pp).owner(); + word patchName = zoneNames[zoneI] + "_" + "side"; + + zoneSidePatch[zoneI] = addPatch<polyPatch> + ( + mesh.boundaryMesh(), + patchName, + newPatches + ); + + Pout<< zoneSidePatch[zoneI] << '\t' << patchName << nl; + + nSide++; + } + } + Pout<< "Added " << nSide << " zone-side patches." << nl + << endl; +} + + +void addInterZonePatches +( + const fvMesh& mesh, + const wordList& zoneNames, + const bool oneD, + + labelList& zoneZonePatch_min, + labelList& zoneZonePatch_max, + DynamicList<polyPatch*>& newPatches +) +{ + Pout<< "Adding inter-zone patches:" << nl << nl + << "patchID\tpatch" << nl + << "-------\t-----" + << endl; + + dictionary transformDict; + transformDict.add + ( + "transform", + cyclicPolyPatch::transformTypeNames[cyclicPolyPatch::NOORDERING] + ); - forAll(pp.faceEdges(), faceI) + label nInter = 0; + if (!oneD) + { + forAll(zoneZonePatch_min, minZone) + { + for (label maxZone = minZone; maxZone < zoneNames.size(); maxZone++) { - const labelList& fEdges = pp.faceEdges()[faceI]; - forAll(fEdges, fp) - { - label meshEdgeI = pp.meshEdges()[fEdges[fp]]; - if (meshToExtrudEdge.found(meshEdgeI)) - { - const face& f = pp.localFaces()[faceI]; - label fp0 = fp; - label fp1 = f.fcIndex(fp0); - label mp0 = pp[faceI][fp0]; - label mp1 = pp[faceI][fp1]; + label index = minZone*zoneNames.size()+maxZone; + if (zoneZonePatch_min[index] > 0) + { + word minToMax = + zoneNames[minZone] + + "_to_" + + zoneNames[maxZone]; + word maxToMin = + zoneNames[maxZone] + + "_to_" + + zoneNames[minZone]; - const vector& nA = + { + transformDict.set("neighbourPatch", maxToMin); + zoneZonePatch_min[index] = + addPatch<nonuniformTransformCyclicPolyPatch> ( - isOwner - ? edgeNormals0[meshEdgeI] - : edgeNormals1[meshEdgeI] + mesh.boundaryMesh(), + minToMax, + transformDict, + newPatches ); - - const vector& nB = + Pout<< zoneZonePatch_min[index] << '\t' << minToMax + << nl; + nInter++; + } + { + transformDict.set("neighbourPatch", minToMax); + zoneZonePatch_max[index] = + addPatch<nonuniformTransformCyclicPolyPatch> ( - isOwner - ? edgeNormals1[meshEdgeI] - : edgeNormals0[meshEdgeI] + mesh.boundaryMesh(), + maxToMin, + transformDict, + newPatches ); - - // Find corresponding face and indices. - { - label exEdgeI = meshToExtrudEdge[meshEdgeI]; - const labelList& eFaces = - extrudePatch.edgeFaces()[exEdgeI]; - // Use 0th face. - label exFaceI = eFaces[0]; - const face& exF = extrudePatch[exFaceI]; - const face& exRegions = pointRegions[exFaceI]; - // Find points - label r0 = exRegions[findIndex(exF, mp0)]; - regionNormals[r0] = nA; - label r1 = exRegions[findIndex(exF, mp1)]; - regionNormals[r1] = nB; - } + Pout<< zoneZonePatch_max[index] << '\t' << maxToMin + << nl; + nInter++; } + } } } } + Pout<< "Added " << nInter << " inter-zone patches." << nl + << endl; } @@ -914,23 +1408,83 @@ tmp<pointField> calcOffset } +void setCouplingInfo +( + fvMesh& mesh, + const labelList& zoneToPatch, + const word& sampleRegion, + const mappedWallPolyPatch::sampleMode mode, + const List<pointField>& offsets +) +{ + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + + List<polyPatch*> newPatches(patches.size(), NULL); + + forAll(zoneToPatch, zoneI) + { + label patchI = zoneToPatch[zoneI]; + + const polyPatch& pp = patches[patchI]; + + if (isA<mappedWallPolyPatch>(pp)) + { + newPatches[patchI] = new mappedWallPolyPatch + ( + pp.name(), + pp.size(), + pp.start(), + patchI, + sampleRegion, // sampleRegion + mode, // sampleMode + pp.name(), // samplePatch + offsets[zoneI], // offset + patches + ); + } + } + + forAll(newPatches, patchI) + { + if (!newPatches[patchI]) + { + newPatches[patchI] = patches[patchI].clone(patches).ptr(); + } + } + + mesh.removeFvBoundary(); + mesh.addFvPatches(newPatches, true); +} + // Main program: int main(int argc, char *argv[]) { - argList::noParallel(); argList::addNote("Create region mesh by extruding a faceZone"); #include "addRegionOption.H" #include "addOverwriteOption.H" - argList::addOption("dict", "name", "specify alternative dictionary"); - argList::addBoolOption("AMI", "apply mapped AMI boundary type"); + argList::addNote("Create region mesh by extruding a faceZone"); #include "setRootCase.H" #include "createTime.H" #include "createNamedMesh.H" + if (mesh.boundaryMesh().checkParallelSync(true)) + { + List<wordList> allNames(Pstream::nProcs()); + allNames[Pstream::myProcNo()] = mesh.boundaryMesh().names(); + Pstream::gatherList(allNames); + Pstream::scatterList(allNames); + + FatalErrorIn(args.executable()) + << "Patches are not synchronised on all processors." + << " Per processor patches " << allNames + << exit(FatalError); + } + + const word oldInstance = mesh.pointsInstance(); bool overwrite = args.optionFound("overwrite"); const word dictName @@ -966,7 +1520,7 @@ int main(int argc, char *argv[]) const Switch oneD(dict.lookup("oneD")); const Switch adaptMesh(dict.lookup("adaptMesh")); - Info<< "Extruding zones " << zoneNames + Pout<< "Extruding zones " << zoneNames << " on mesh " << regionName << " into shell mesh " << shellRegionName << endl; @@ -982,6 +1536,53 @@ int main(int argc, char *argv[]) + //// Read objects in time directory + //IOobjectList objects(mesh, runTime.timeName()); + // + //// Read vol fields. + // + //PtrList<volScalarField> vsFlds; + //ReadFields(mesh, objects, vsFlds); + // + //PtrList<volVectorField> vvFlds; + //ReadFields(mesh, objects, vvFlds); + // + //PtrList<volSphericalTensorField> vstFlds; + //ReadFields(mesh, objects, vstFlds); + // + //PtrList<volSymmTensorField> vsymtFlds; + //ReadFields(mesh, objects, vsymtFlds); + // + //PtrList<volTensorField> vtFlds; + //ReadFields(mesh, objects, vtFlds); + // + //// Read surface fields. + // + //PtrList<surfaceScalarField> ssFlds; + //ReadFields(mesh, objects, ssFlds); + // + //PtrList<surfaceVectorField> svFlds; + //ReadFields(mesh, objects, svFlds); + // + //PtrList<surfaceSphericalTensorField> sstFlds; + //ReadFields(mesh, objects, sstFlds); + // + //PtrList<surfaceSymmTensorField> ssymtFlds; + //ReadFields(mesh, objects, ssymtFlds); + // + //PtrList<surfaceTensorField> stFlds; + //ReadFields(mesh, objects, stFlds); + // + //// Read point fields. + // + //PtrList<pointScalarField> psFlds; + //ReadFields(pointMesh::New(mesh), objects, psFlds); + // + //PtrList<pointVectorField> pvFlds; + //ReadFields(pointMesh::New(mesh), objects, pvFlds); + + + // Create dummy fv* files createDummyFvMeshFiles(mesh, shellRegionName); @@ -996,7 +1597,7 @@ int main(int argc, char *argv[]) { meshInstance = oldInstance; } - Info<< "Writing meshes to " << meshInstance << nl << endl; + Pout<< "Writing meshes to " << meshInstance << nl << endl; const polyBoundaryMesh& patches = mesh.boundaryMesh(); @@ -1093,7 +1694,7 @@ int main(int argc, char *argv[]) const labelListList& edgeFaces = extrudePatch.edgeFaces(); - Info<< "extrudePatch :" + Pout<< "extrudePatch :" << " faces:" << extrudePatch.size() << " points:" << extrudePatch.nPoints() << " edges:" << extrudePatch.nEdges() @@ -1149,112 +1750,92 @@ int main(int argc, char *argv[]) const faceZone& fz = faceZones[zoneIDs[i]]; if (isInternal[i]) { - Info<< "FaceZone " << fz.name() << " has internal faces" << endl; + Pout<< "FaceZone " << fz.name() << " has internal faces" << endl; } else { - Info<< "FaceZone " << fz.name() << " has boundary faces" << endl; + Pout<< "FaceZone " << fz.name() << " has boundary faces" << endl; } } - Info<< endl; + Pout<< endl; - // Add interface patches - // ~~~~~~~~~~~~~~~~~~~~~ - // Note that these actually get added to the original mesh - // so the shell mesh creation copies them. They then get removed - // from the original mesh. - - Info<< "Adding coupling patches:" << nl << nl - << "patchID\tpatch\ttype" << nl - << "-------\t-----\t----" - << endl; - labelList interRegionTopPatch(zoneNames.size()); - labelList interRegionBottomPatch(zoneNames.size()); - label nCoupled = 0; - forAll(zoneIDs, i) + DynamicList<polyPatch*> regionPatches(patches.size()); + // Copy all non-local patches since these are used on boundary edges of + // the extrusion + forAll(patches, patchI) { - word interName(regionName+"_to_"+shellRegionName+'_'+zoneNames[i]); - - if (isInternal[i]) + if (!isA<processorPolyPatch>(patches[patchI])) { - interRegionTopPatch[i] = addPatch<mappedWallPolyPatch> - ( - mesh, - interName + "_top" - ); - nCoupled++; - Info<< interRegionTopPatch[i] - << '\t' << patches[interRegionTopPatch[i]].name() - << '\t' << patches[interRegionTopPatch[i]].type() - << nl; - - interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch> + label newPatchI = regionPatches.size(); + regionPatches.append ( - mesh, - interName + "_bottom" + patches[patchI].clone + ( + patches, + newPatchI, + 0, // size + 0 // start + ).ptr() ); - nCoupled++; - Info<< interRegionBottomPatch[i] - << '\t' << patches[interRegionBottomPatch[i]].name() - << '\t' << patches[interRegionBottomPatch[i]].type() - << nl; } - else if (zoneShadowNames.size() == 0) - { - interRegionTopPatch[i] = addPatch<polyPatch> - ( - mesh, - zoneNames[i] + "_top" - ); - nCoupled++; - Info<< interRegionTopPatch[i] - << '\t' << patches[interRegionTopPatch[i]].name() - << '\t' << patches[interRegionTopPatch[i]].type() - << nl; + } - interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch> - ( - mesh, - interName - ); - nCoupled++; - Info<< interRegionBottomPatch[i] - << '\t' << patches[interRegionBottomPatch[i]].name() - << '\t' << patches[interRegionBottomPatch[i]].type() - << nl; - } - else if (zoneShadowNames.size() > 0) //patch using shadow face zones. - { - interRegionTopPatch[i] = addPatch<mappedWallPolyPatch> - ( - mesh, - zoneShadowNames[i] + "_top" - ); - nCoupled++; - Info<< interRegionTopPatch[i] - << '\t' << patches[interRegionTopPatch[i]].name() - << '\t' << patches[interRegionTopPatch[i]].type() - << nl; - interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch> - ( - mesh, - interName - ); - nCoupled++; - Info<< interRegionBottomPatch[i] - << '\t' << patches[interRegionBottomPatch[i]].name() - << '\t' << patches[interRegionBottomPatch[i]].type() - << nl; + // Add interface patches + // ~~~~~~~~~~~~~~~~~~~~~ + + // From zone to interface patch (region side) + labelList interRegionTopPatch(zoneNames.size()); + labelList interRegionBottomPatch(zoneNames.size()); + + addCouplingPatches + ( + mesh, + regionName, + shellRegionName, + zoneNames, + zoneShadowNames, + isInternal, + + regionPatches, + interRegionTopPatch, + interRegionBottomPatch + ); + + // From zone to interface patch (mesh side) + labelList interMeshTopPatch(zoneNames.size()); + labelList interMeshBottomPatch(zoneNames.size()); + + if (adaptMesh) + { + DynamicList<polyPatch*> newPatches(patches.size()); + forAll(patches, patchI) + { + newPatches.append(patches[patchI].clone(patches).ptr()); } + addCouplingPatches + ( + mesh, + regionName, + shellRegionName, + zoneNames, + zoneShadowNames, + isInternal, + + newPatches, + interMeshTopPatch, + interMeshBottomPatch + ); + mesh.clearOut(); + mesh.removeFvBoundary(); + mesh.addFvPatches(newPatches, true); } - Info<< "Added " << nCoupled << " inter-region patches." << nl - << endl; + // Patch per extruded face labelList extrudeTopPatchID(extrudePatch.size()); labelList extrudeBottomPatchID(extrudePatch.size()); @@ -1294,136 +1875,79 @@ int main(int argc, char *argv[]) zoneZonePatch_min // reuse for counting ); - // Now check which patches to add. - Info<< "Adding patches for edges on zones:" << nl << nl - << "patchID\tpatch" << nl - << "-------\t-----" - << endl; - - label nSide = 0; - - forAll(zoneSidePatch, zoneI) - { - if (oneD) - { - // Reuse single empty patch. - word patchType = dict.lookup("oneDPolyPatchType"); - word patchName; - if (patchType == "emptyPolyPatch") - { - patchName = "oneDEmptyPatch"; - zoneSidePatch[zoneI] = addPatch<emptyPolyPatch> - ( - mesh, - patchName - ); - } - else if (patchType == "wedgePolyPatch") - { - patchName = "oneDWedgePatch"; - zoneSidePatch[zoneI] = addPatch<wedgePolyPatch> - ( - mesh, - patchName - ); - } - else - { - FatalErrorIn(args.executable()) - << "Type " << patchType << " does not exist " - << exit(FatalError); - } - - Info<< zoneSidePatch[zoneI] << '\t' << patchName << nl; - - nSide++; - } - else if (zoneSidePatch[zoneI] > 0) - { - word patchName = faceZones[zoneI].name() + "_" + "side"; - - zoneSidePatch[zoneI] = addPatch<polyPatch> - ( - mesh, - patchName - ); - - Info<< zoneSidePatch[zoneI] << '\t' << patchName << nl; + // Now we'll have: + // zoneSidePatch[zoneA] : number of faces needed on the side of zoneA + // zoneZonePatch_min[zoneA,zoneB] : number of faces needed inbetween A,B - nSide++; - } - } - Info<< "Added " << nSide << " zone-edge patches." << nl - << endl; + // Add the zone-side patches. + addZoneSidePatches + ( + mesh, + (oneD ? dict.lookup("oneDPolyPatchType") : word::null), + zoneNames, + regionPatches, + zoneSidePatch + ); - Info<< "Adding inter-zone patches:" << nl << nl - << "patchID\tpatch" << nl - << "-------\t-----" - << endl; - dictionary transformDict; - transformDict.add + // Add the patches inbetween zones + addInterZonePatches ( - "transform", - cyclicPolyPatch::transformTypeNames[cyclicPolyPatch::NOORDERING] + mesh, + zoneNames, + oneD, + + zoneZonePatch_min, + zoneZonePatch_max, + regionPatches ); - label nInter = 0; - if (!oneD) - { - forAll(zoneZonePatch_min, minZone) - { - for (label maxZone = minZone; maxZone < faceZones.size(); maxZone++) - { - label index = minZone*faceZones.size()+maxZone; - if (zoneZonePatch_min[index] > 0) - { - word minToMax = - faceZones[minZone].name() - + "_to_" - + faceZones[maxZone].name(); - word maxToMin = - faceZones[maxZone].name() - + "_to_" - + faceZones[minZone].name(); + // Additionally check if any interprocessor patches need to be added. + // Reuses addPatchCellLayer functionality. + // Note: does not handle edges with > 2 faces? + labelList sidePatchID; + addProcPatches + ( + mesh, + extrudePatch, + extrudeMeshFaces, - { - transformDict.set("neighbourPatch", maxToMin); - zoneZonePatch_min[index] = - addPatch<nonuniformTransformCyclicPolyPatch> - ( - mesh, - minToMax, - transformDict - ); - Info<< zoneZonePatch_min[index] << '\t' << minToMax - << nl; - nInter++; - } - { - transformDict.set("neighbourPatch", minToMax); - zoneZonePatch_max[index] = - addPatch<nonuniformTransformCyclicPolyPatch> - ( - mesh, - maxToMin, - transformDict - ); - Info<< zoneZonePatch_max[index] << '\t' << maxToMin - << nl; - nInter++; - } + sidePatchID, + regionPatches + ); - } - } - } - } - Info<< "Added " << nInter << " inter-zone patches." << nl - << endl; +// // Add all the newPatches to the mesh and fields +// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// { +// forAll(newPatches, patchI) +// { +// Pout<< "Adding patch " << patchI +// << " name:" << newPatches[patchI]->name() +// << endl; +// } +// //label nOldPatches = mesh.boundary().size(); +// mesh.clearOut(); +// mesh.removeFvBoundary(); +// mesh.addFvPatches(newPatches, true); +// //// Add calculated fvPatchFields for the added patches +// //for +// //( +// // label patchI = nOldPatches; +// // patchI < mesh.boundary().size(); +// // patchI++ +// //) +// //{ +// // Pout<< "ADDing calculated to patch " << patchI +// // << endl; +// // addCalculatedPatchFields(mesh); +// //} +// //Pout<< "** Added " << mesh.boundary().size()-nOldPatches +// // << " patches." << endl; +// } // Set patches to use for edges to be extruded into boundary faces @@ -1433,7 +1957,7 @@ int main(int argc, char *argv[]) // If empty size create an internal face. labelListList extrudeEdgePatches(extrudePatch.nEdges()); - // Is edge an non-manifold edge + // Is edge a non-manifold edge PackedBoolList nonManifoldEdge(extrudePatch.nEdges()); // Note: logic has to be same as in countExtrudePatches. @@ -1483,6 +2007,15 @@ int main(int argc, char *argv[]) nonManifoldEdge[edgeI] = 1; } } + else if (sidePatchID[edgeI] != -1) + { + // Coupled extrusion + ePatches.setSize(eFaces.size()); + forAll(eFaces, i) + { + ePatches[i] = sidePatchID[edgeI]; + } + } else { label faceI = findUncoveredPatchFace @@ -1494,8 +2027,12 @@ int main(int argc, char *argv[]) if (faceI != -1) { - label patchI = mesh.boundaryMesh().whichPatch(faceI); - ePatches.setSize(eFaces.size(), patchI); + label newPatchI = findPatchID + ( + regionPatches, + patches[patches.whichPatch(faceI)].name() + ); + ePatches.setSize(eFaces.size(), newPatchI); } else { @@ -1515,76 +2052,102 @@ int main(int argc, char *argv[]) // ~~~~~~~~~~~~~~~~~~~~ // Per face, per point the region number. - faceList pointRegions(extrudePatch.size()); - // Per region the originating point - labelList regionPoints; + faceList pointGlobalRegions; + faceList pointLocalRegions; + labelList localToGlobalRegion; createShellMesh::calcPointRegions ( + mesh.globalData(), extrudePatch, nonManifoldEdge, - pointRegions, - regionPoints + pointGlobalRegions, + pointLocalRegions, + localToGlobalRegion ); - - // Calculate a normal per region - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - vectorField regionNormals(regionPoints.size(), vector::zero); - vectorField regionCentres(regionPoints.size(), vector::zero); - labelList nRegionFaces(regionPoints.size(), 0); - - forAll(pointRegions, faceI) + // Per local region an originating point + labelList localRegionPoints(localToGlobalRegion.size()); + forAll(pointLocalRegions, faceI) { - const face& f = extrudeFaces[faceI]; - - forAll(f, fp) + const face& f = extrudePatch.localFaces()[faceI]; + const face& pRegions = pointLocalRegions[faceI]; + forAll(pRegions, fp) { - label region = pointRegions[faceI][fp]; - regionNormals[region] += extrudePatch.faceNormals()[faceI]; - regionCentres[region] += extrudePatch.faceCentres()[faceI]; - nRegionFaces[region]++; + localRegionPoints[pRegions[fp]] = f[fp]; } } - regionNormals /= mag(regionNormals); - forAll(regionCentres, regionI) + + // Calculate region normals by reducing local region normals + pointField localRegionNormals(localToGlobalRegion.size()); { - regionCentres[regionI] /= nRegionFaces[regionI]; - } + pointField localSum(localToGlobalRegion.size(), vector::zero); + labelList localNum(localToGlobalRegion.size(), 0); + forAll(pointLocalRegions, faceI) + { + const face& pRegions = pointLocalRegions[faceI]; + forAll(pRegions, fp) + { + label localRegionI = pRegions[fp]; + localSum[localRegionI] += extrudePatch.faceNormals()[faceI]; + localNum[localRegionI]++; + } + } + + Map<point> globalSum(2*localToGlobalRegion.size()); + Map<label> globalNum(2*localToGlobalRegion.size()); + + forAll(localSum, localRegionI) + { + label globalRegionI = localToGlobalRegion[localRegionI]; + globalSum.insert(globalRegionI, localSum[localRegionI]); + globalNum.insert(globalRegionI, localNum[localRegionI]); + } + + // Reduce + Pstream::mapCombineGather(globalSum, plusEqOp<point>()); + Pstream::mapCombineScatter(globalSum); + Pstream::mapCombineGather(globalNum, plusEqOp<label>()); + Pstream::mapCombineScatter(globalNum); + + forAll(localToGlobalRegion, localRegionI) + { + label globalRegionI = localToGlobalRegion[localRegionI]; + localRegionNormals[localRegionI] = + globalSum[globalRegionI] + / globalNum[globalRegionI]; + } + localRegionNormals /= mag(localRegionNormals); + } - // Constrain&sync normals on points that are on coupled patches. - constrainCoupledNormals - ( - mesh, - extrudePatch, - extrudeMeshEdges, - pointRegions, - regionNormals - ); // For debugging: dump hedgehog plot of normals if (false) { - OFstream str(runTime.path()/"regionNormals.obj"); + OFstream str(runTime.path()/"localRegionNormals.obj"); label vertI = 0; scalar thickness = model().sumThickness(1); - forAll(pointRegions, faceI) + forAll(pointLocalRegions, faceI) { const face& f = extrudeFaces[faceI]; forAll(f, fp) { - label region = pointRegions[faceI][fp]; + label region = pointLocalRegions[faceI][fp]; const point& pt = extrudePoints[f[fp]]; meshTools::writeOBJ(str, pt); vertI++; - meshTools::writeOBJ(str, pt+thickness*regionNormals[region]); + meshTools::writeOBJ + ( + str, + pt+thickness*localRegionNormals[region] + ); vertI++; str << "l " << vertI-1 << ' ' << vertI << nl; } @@ -1593,11 +2156,18 @@ int main(int argc, char *argv[]) // Use model to create displacements of first layer - vectorField firstDisp(regionNormals.size()); + vectorField firstDisp(localRegionNormals.size()); forAll(firstDisp, regionI) { - const point& regionPt = regionCentres[regionI]; - const vector& n = regionNormals[regionI]; + //const point& regionPt = regionCentres[regionI]; + const point& regionPt = extrudePatch.points() + [ + extrudePatch.meshPoints() + [ + localRegionPoints[regionI] + ] + ]; + const vector& n = localRegionNormals[regionI]; firstDisp[regionI] = model()(regionPt, n, 1) - regionPt; } @@ -1606,13 +2176,47 @@ int main(int argc, char *argv[]) // Create a new mesh // ~~~~~~~~~~~~~~~~~ - createShellMesh extruder(extrudePatch, pointRegions, regionPoints); + createShellMesh extruder + ( + extrudePatch, + pointLocalRegions, + localRegionPoints + ); - autoPtr<fvMesh> regionMeshPtr; autoPtr<mapPolyMesh> shellMap; + fvMesh regionMesh + ( + IOobject + ( + shellRegionName, + meshInstance, + runTime, + IOobject::NO_READ, + IOobject::AUTO_WRITE, + false + ), + xferCopy(pointField()), + xferCopy(faceList()), + xferCopy(labelList()), + xferCopy(labelList()), + false + ); + + // Add the new patches + forAll(regionPatches, patchI) + { + regionPatches[patchI] = regionPatches[patchI]->clone + ( + regionMesh.boundaryMesh() + ).ptr(); + } + regionMesh.clearOut(); + regionMesh.removeFvBoundary(); + regionMesh.addFvPatches(regionPatches, true); + { - polyTopoChange meshMod(mesh.boundaryMesh().size()); + polyTopoChange meshMod(regionPatches.size()); extruder.setRefinement ( @@ -1625,22 +2229,12 @@ int main(int argc, char *argv[]) meshMod ); - shellMap = meshMod.makeMesh + shellMap = meshMod.changeMesh ( - regionMeshPtr, // mesh to create - IOobject - ( - shellRegionName, - meshInstance, - runTime, - IOobject::NO_READ, - IOobject::AUTO_WRITE, - false - ), - mesh // mesh to get original patch info from + regionMesh, // mesh to change + false // inflate ); } - fvMesh& regionMesh = regionMeshPtr(); // Necessary? regionMesh.setInstance(meshInstance); @@ -1650,82 +2244,92 @@ int main(int argc, char *argv[]) extruder.updateMesh(shellMap); - // Change top and bottom boundary conditions on regionMesh - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Calculate offsets from shell mesh back to original mesh + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - // Save offsets from shell mesh back to original mesh List<pointField> topOffsets(zoneIDs.size()); List<pointField> bottomOffsets(zoneIDs.size()); + forAll(regionMesh.boundaryMesh(), patchI) { - const polyBoundaryMesh& regionPatches = regionMesh.boundaryMesh(); - List<polyPatch*> newPatches(regionPatches.size()); - forAll(regionPatches, patchI) + const polyPatch& pp = regionMesh.boundaryMesh()[patchI]; + + if + ( + isA<mappedWallPolyPatch>(pp) + && (findIndex(interRegionTopPatch, patchI) != -1) + ) { - const polyPatch& pp = regionPatches[patchI]; + label zoneI = findIndex(interRegionTopPatch, patchI); + topOffsets[zoneI] = calcOffset(extrudePatch, extruder, pp); + } + else if + ( + isA<mappedWallPolyPatch>(pp) + && (findIndex(interRegionBottomPatch, patchI) != -1) + ) + { + label zoneI = findIndex(interRegionBottomPatch, patchI); - if - ( - isA<mappedWallPolyPatch>(pp) - && (findIndex(interRegionTopPatch, patchI) != -1) - ) - { - label index = findIndex(interRegionTopPatch, patchI); + bottomOffsets[zoneI] = calcOffset(extrudePatch, extruder, pp); + } + } - topOffsets[index] = calcOffset(extrudePatch, extruder, pp); - newPatches[patchI] = new mappedWallPolyPatch - ( - pp.name(), - pp.size(), - pp.start(), - patchI, - regionName, // sampleRegion - sampleMode, // sampleMode - pp.name(), // samplePatch - topOffsets[index], // offset - patches - ); - } - else if - ( - isA<mappedWallPolyPatch>(pp) - && (findIndex(interRegionBottomPatch, patchI) != -1) - ) - { - label index = findIndex(interRegionBottomPatch, patchI); + // Change top and bottom boundary conditions on regionMesh + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + { + // Correct top patches for offset + setCouplingInfo + ( + regionMesh, + interRegionTopPatch, + regionName, // name of main mesh + sampleMode, // sampleMode + topOffsets + ); - bottomOffsets[index] = calcOffset(extrudePatch, extruder, pp); + // Correct bottom patches for offset + setCouplingInfo + ( + regionMesh, + interRegionBottomPatch, + regionName, + sampleMode, // sampleMode + bottomOffsets + ); - newPatches[patchI] = new mappedWallPolyPatch - ( - pp.name(), - pp.size(), - pp.start(), - patchI, - regionName, // sampleRegion - sampleMode, // sampleMode - pp.name(), // samplePatch - bottomOffsets[index], // offset - patches - ); - } - else - { - newPatches[patchI] = pp.clone - ( - regionPatches, - patchI, - pp.size(), - pp.start() - ).ptr(); - } - } - regionMesh.removeFvBoundary(); - regionMesh.addFvPatches(newPatches, true); + // Remove any unused patches deleteEmptyPatches(regionMesh); } + // Change top and bottom boundary conditions on main mesh + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + if (adaptMesh) + { + // Correct top patches for offset + setCouplingInfo + ( + mesh, + interMeshTopPatch, + shellRegionName, // name of shell mesh + sampleMode, // sampleMode + -topOffsets + ); + + // Correct bottom patches for offset + setCouplingInfo + ( + mesh, + interMeshBottomPatch, + shellRegionName, + sampleMode, + -bottomOffsets + ); + } + // Write addressing from region mesh back to originating patch @@ -1799,7 +2403,7 @@ int main(int argc, char *argv[]) "point to patch point addressing"; - Info<< "Writing mesh " << regionMesh.name() + Pout<< "Writing mesh " << regionMesh.name() << " to " << regionMesh.facesInstance() << nl << endl; @@ -1847,7 +2451,7 @@ int main(int argc, char *argv[]) mesh.faceOwner()[meshFaceI],// owner -1, // neighbour false, // face flip - extrudeBottomPatchID[zoneFaceI],// patch for face + interMeshBottomPatch[zoneI],// patch for face zoneI, // zone for face flip // face flip in zone ); @@ -1861,7 +2465,7 @@ int main(int argc, char *argv[]) mesh.faceNeighbour()[meshFaceI],// owner -1, // neighbour true, // face flip - extrudeBottomPatchID[zoneFaceI],// patch for face + interMeshBottomPatch[zoneI], // patch for face zoneI, // zone for face !flip // face flip in zone ); @@ -1886,7 +2490,7 @@ int main(int argc, char *argv[]) mesh.faceOwner()[meshFaceI],// owner -1, // neighbour false, // face flip - extrudeTopPatchID[zoneFaceI],// patch for face + interMeshTopPatch[zoneI], // patch for face zoneI, // zone for face flip // face flip in zone ); @@ -1900,7 +2504,7 @@ int main(int argc, char *argv[]) mesh.faceNeighbour()[meshFaceI],// owner -1, // neighbour true, // face flip - extrudeTopPatchID[zoneFaceI], // patch for face + interMeshTopPatch[zoneI], // patch for face zoneI, // zone for face !flip // face flip in zone ); @@ -1914,6 +2518,7 @@ int main(int argc, char *argv[]) forAll(extrudeMeshFaces, zoneFaceI) { label meshFaceI = extrudeMeshFaces[zoneFaceI]; + label zoneI = zoneID[zoneFaceI]; bool flip = zoneFlipMap[zoneFaceI]; const face& f = mesh.faces()[meshFaceI]; @@ -1930,7 +2535,7 @@ int main(int argc, char *argv[]) -1, // master edge meshFaceI, // master face true, // flip flux - extrudeTopPatchID[zoneFaceI], // patch for face + interMeshTopPatch[zoneI], // patch for face -1, // zone for face false //face flip in zone ); @@ -1947,7 +2552,7 @@ int main(int argc, char *argv[]) -1, // master edge meshFaceI, // master face false, // flip flux - extrudeTopPatchID[zoneFaceI], // patch for face + interMeshTopPatch[zoneI], // patch for face -1, // zone for face false // zone flip ); @@ -1973,84 +2578,9 @@ int main(int argc, char *argv[]) } mesh.setInstance(meshInstance); - } - - - - // Change master and slave boundary conditions on originating mesh - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - if (adaptMesh) - { - const polyBoundaryMesh& patches = mesh.boundaryMesh(); - List<polyPatch*> newPatches(patches.size()); - forAll(patches, patchI) - { - const polyPatch& pp = patches[patchI]; - - if - ( - isA<mappedWallPolyPatch>(pp) - && (findIndex(interRegionTopPatch, patchI) != -1) - ) - { - label index = findIndex(interRegionTopPatch, patchI); - newPatches[patchI] = new mappedWallPolyPatch - ( - pp.name(), - pp.size(), - pp.start(), - patchI, - shellRegionName, // sampleRegion - sampleMode, // sampleMode - pp.name(), // samplePatch - -topOffsets[index], // offset - patches - ); - } - else if - ( - isA<mappedWallPolyPatch>(pp) - && (findIndex(interRegionBottomPatch, patchI) != -1) - ) - { - label index = findIndex(interRegionBottomPatch, patchI); - - newPatches[patchI] = new mappedWallPolyPatch - ( - pp.name(), - pp.size(), - pp.start(), - patchI, - shellRegionName, // sampleRegion - sampleMode, // sampleMode - pp.name(), // samplePatch - -bottomOffsets[index], // offset - patches - ); - } - else - { - newPatches[patchI] = pp.clone - ( - patches, - patchI, - pp.size(), - pp.start() - ).ptr(); - } - } - mesh.removeFvBoundary(); - mesh.addFvPatches(newPatches, true); - - // Remove any unused patches - deleteEmptyPatches(mesh); - } - if (adaptMesh) - { - Info<< "Writing mesh " << mesh.name() + Pout<< "Writing mesh " << mesh.name() << " to " << mesh.facesInstance() << nl << endl; diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.H b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.H deleted file mode 100644 index 4bc46bfb245b27eab5ac3a76280238ef79958ab8..0000000000000000000000000000000000000000 --- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.H +++ /dev/null @@ -1,205 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation - \\/ 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 3 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, see <http://www.gnu.org/licenses/>. - -Class - Foam::patchPointEdgeCirculator - -Description - Walks from starting edge/face around point on patch. - - -# Use in-place: \n - \code - patchPointEdgeCirculator circ(..); - - // Walk - do - { - Info<< "edge:" << circ.edgeID() << endl; - ++circ; - } - while (circ != circ.end()); - \endcode - - -# Use like STL iterator: \n - \code - patchPointEdgeCirculator circ(..); - for - ( - patchPointEdgeCirculator iter = circ.begin(); - iter != circ.end(); - ++iter - ) - { - Info<< "edge:" << iter.edgeID() << endl; - } - \endcode - - -SourceFiles - patchPointEdgeCirculator.C - -\*---------------------------------------------------------------------------*/ - -#ifndef patchPointEdgeCirculator_H -#define patchPointEdgeCirculator_H - -#include "face.H" -#include "ListOps.H" -#include "primitiveFacePatch.H" -#include "PackedBoolList.H" -#include "InfoProxy.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes - -/*---------------------------------------------------------------------------*\ - Class patchPointEdgeCirculator Declaration -\*---------------------------------------------------------------------------*/ - -class patchPointEdgeCirculator -{ - // Static data members - - //- end iterator - static const patchPointEdgeCirculator endConstIter; - - - // Private data - - //- Patch - const primitiveFacePatch& patch_; - - const PackedBoolList& nonManifoldEdge_; - - //- Current edge - label edgeID_; - - //- Current direction (face, expressed as index into edgeFaces()[edgeID] - label index_; - - //- Point - label pointID_; - - //- Starting edge so we know when to stop. - label startEdgeID_; - - - // Private Member Functions - - //- Set to end() iterator - inline void setEnd(); - - //- Set edgeID_ to be the other edge on the face that uses the - // point. - inline void otherEdge(const label cellI); - -public: - - // Constructors - - //- Construct from components - inline patchPointEdgeCirculator - ( - const primitiveFacePatch&, - const PackedBoolList& nonManifoldEdge, - const label edgeID, - const label index, - const label pointID - ); - - //- Construct as copy - inline patchPointEdgeCirculator(const patchPointEdgeCirculator&); - - - // Member Functions - - inline label edgeID() const; - - inline label index() const; - - inline label pointID() const; - - //- Helper: get face according to the index. - // Returns -1 if on end. - inline label faceID() const; - - //- Walk back until non-manifold edge (if any) or minimum edge. - inline void setCanonical(); - - - // Member Operators - - inline void operator=(const patchPointEdgeCirculator& iter); - - inline bool operator==(const patchPointEdgeCirculator& iter) const; - - inline bool operator!=(const patchPointEdgeCirculator& iter) const; - - //- Step to next face. - inline patchPointEdgeCirculator& operator++(); - - //- iterator set to the beginning face. For internal edges this is - // the current face. For boundary edges this is the first boundary face - // reached from walking back (i.e. in opposite direction to ++) - inline patchPointEdgeCirculator begin() const; - inline patchPointEdgeCirculator cbegin() const; - - //- iterator set to beyond the end of the walk. - inline const patchPointEdgeCirculator& end() const; - inline const patchPointEdgeCirculator& cend() const; - - // Info - - //- Return info proxy. - // Used to print token information to a stream - InfoProxy<patchPointEdgeCirculator> info() const - { - return *this; - } - - friend Ostream& operator<< - ( - Ostream&, - const InfoProxy<patchPointEdgeCirculator>& - ); -}; - -Ostream& operator<<(Ostream&, const InfoProxy<patchPointEdgeCirculator>&); - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#include "patchPointEdgeCirculatorI.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculatorI.H b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculatorI.H deleted file mode 100644 index 55844065a7aea0406f9ba966ac3ed8fed765f363..0000000000000000000000000000000000000000 --- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculatorI.H +++ /dev/null @@ -1,308 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation - \\/ 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 3 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, see <http://www.gnu.org/licenses/>. - -\*---------------------------------------------------------------------------*/ - - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::patchPointEdgeCirculator::setEnd() -{ - edgeID_ = -1; - pointID_ = -1; -} - - -// Cross face to other edge on point -void Foam::patchPointEdgeCirculator::otherEdge(const label faceI) -{ - const labelList& fEdges = patch_.faceEdges()[faceI]; - const face& f = patch_.localFaces()[faceI]; - label fp = findIndex(f, pointID_); - - if (fEdges[fp] == edgeID_) - { - edgeID_ = fEdges[f.rcIndex(fp)]; - } - else - { - // Check for now - if (fEdges[f.rcIndex(fp)] != edgeID_) - { - FatalErrorIn("patchPointEdgeCirculator::otherEdge(const label)") - << "face:" << faceI - << " verts:" << f - << " edges:" << fEdges - << " looking for other edge than " << edgeID_ - << abort(FatalError); - } - edgeID_ = fEdges[fp]; - } -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -//- Construct from components -Foam::patchPointEdgeCirculator::patchPointEdgeCirculator -( - const primitiveFacePatch& patch, - const PackedBoolList& nonManifoldEdge, - const label edgeID, - const label index, - - const label pointID -) -: - patch_(patch), - nonManifoldEdge_(nonManifoldEdge), - edgeID_(edgeID), - index_(index), - pointID_(pointID), - startEdgeID_(edgeID_) -{ - if (edgeID_ != -1) - { - const edge& e = patch_.edges()[edgeID_]; - - if (pointID_ != e[0] && pointID_ != e[1]) - { - FatalErrorIn - ( - "patchPointEdgeCirculator::patchPointEdgeCirculator(..)" - ) << "edge " << edgeID_ << " verts " << e - << " does not contain point " << pointID_ << abort(FatalError); - } - } -} - - -//- Construct copy -Foam::patchPointEdgeCirculator::patchPointEdgeCirculator -( - const patchPointEdgeCirculator& circ -) -: - patch_(circ.patch_), - nonManifoldEdge_(circ.nonManifoldEdge_), - edgeID_(circ.edgeID_), - index_(circ.index_), - pointID_(circ.pointID_), - startEdgeID_(circ.startEdgeID_) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -Foam::label Foam::patchPointEdgeCirculator::edgeID() const -{ - return edgeID_; -} - - -Foam::label Foam::patchPointEdgeCirculator::index() const -{ - return index_; -} - - -Foam::label Foam::patchPointEdgeCirculator::pointID() const -{ - return pointID_; -} - - -Foam::label Foam::patchPointEdgeCirculator::faceID() const -{ - if (edgeID_ != -1 && index_ != -1) - { - return patch_.edgeFaces()[edgeID_][index_]; - } - else - { - return -1; - } -} - - -void Foam::patchPointEdgeCirculator::operator= -( - const patchPointEdgeCirculator& circ -) -{ - edgeID_ = circ.edgeID_; - index_ = circ.index_; - pointID_ = circ.pointID_; - startEdgeID_ = circ.startEdgeID_; -} - - -bool Foam::patchPointEdgeCirculator::operator== -( - const patchPointEdgeCirculator& circ -) const -{ - // Do just enough to have setEnd() produce an iterator equal to end(). - // Could include the direction(i.e. index_) to make sure that two - // circulators around same point but in different direction are not equal. - return edgeID_ == circ.edgeID_ && pointID_ == circ.pointID_; -} - - -bool Foam::patchPointEdgeCirculator::operator!= -( - const patchPointEdgeCirculator& circ -) const -{ - return !(*this == circ); -} - - -void Foam::patchPointEdgeCirculator::setCanonical() -{ - if (edgeID_ == -1) - { - FatalErrorIn("patchPointEdgeCirculator::setCanonical()") - << "Already reached end(). Cannot walk any further." - << abort(FatalError); - } - - label minEdgeID = edgeID_; - label minIndex = index_; - - while (true) - { - if (nonManifoldEdge_[edgeID_]) - { - break; - } - - // Step back - const labelList& eFaces = patch_.edgeFaces()[edgeID_]; - - if (eFaces.size() != 2) - { - FatalErrorIn("patchPointEdgeCirculator::setCanonical()") - << "problem eFaces:" << eFaces << abort(FatalError); - } - - label faceI = (index_ == 0 ? eFaces[1] : eFaces[0]); - - // Step to other edge on face - otherEdge(faceI); - - // Update index - index_ = findIndex(patch_.edgeFaces()[edgeID_], faceI); - - if (edgeID_ < minEdgeID) - { - minEdgeID = edgeID_; - minIndex = index_; - } - - if (edgeID_ == startEdgeID_) - { - edgeID_ = minEdgeID; - index_ = minIndex; - break; - } - } - - startEdgeID_ = edgeID_; -} - - -//- Step to next edge. -Foam::patchPointEdgeCirculator& -Foam::patchPointEdgeCirculator::operator++() -{ - if (index_ == -1) - { - setEnd(); - } - else - { - // Step to other edge on face - label faceI = patch_.edgeFaces()[edgeID_][index_]; - otherEdge(faceI); - - if (edgeID_ == startEdgeID_) - { - setEnd(); - } - else if (nonManifoldEdge_[edgeID_]) - { - // Reached non-manifold edge. Cannot walk further. - // Mark so it gets set to end next time. - index_ = -1; - } - else - { - const labelList& eFaces = patch_.edgeFaces()[edgeID_]; - - if (eFaces.size() != 2) - { - FatalErrorIn("patchPointEdgeCirculator:::operator++()") - << "problem eFaces:" << eFaces << abort(FatalError); - } - // Point to the face that is not faceI - index_ = (eFaces[0] != faceI ? 0 : 1); - } - } - - return *this; -} - - -Foam::patchPointEdgeCirculator Foam::patchPointEdgeCirculator::begin() const -{ - patchPointEdgeCirculator iter(*this); - iter.setCanonical(); - - return iter; -} - - -Foam::patchPointEdgeCirculator Foam::patchPointEdgeCirculator::cbegin() const -{ - patchPointEdgeCirculator iter(*this); - iter.setCanonical(); - - return iter; -} - - -const Foam::patchPointEdgeCirculator& Foam::patchPointEdgeCirculator::end() - const -{ - return endConstIter; -} - -const Foam::patchPointEdgeCirculator& Foam::patchPointEdgeCirculator::cend() - const -{ - return endConstIter; -} - - -// ************************************************************************* // diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C index daf71ac7e97d4b83906d40508cd56de3aa3d029b..a6ae8b76817ae687f0635a670bf737bbc385b389 100644 --- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C +++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C @@ -36,7 +36,7 @@ Foam::label Foam::checkTopology label nTotCells = returnReduce(mesh.cells().size(), sumOp<label>()); // These are actually warnings, not errors. - if (nEmpty % nTotCells) + if (nTotCells && (nEmpty % nTotCells)) { Info<< " ***Total number of faces on empty patches" << " is not divisible by the number of cells in the mesh." @@ -255,7 +255,7 @@ Foam::label Foam::checkTopology { regionSplit rs(mesh); - if (rs.nRegions() == 1) + if (rs.nRegions() <= 1) { Info<< " Number of regions: " << rs.nRegions() << " (OK)." << endl; diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSetDict b/applications/utilities/mesh/manipulation/topoSet/topoSetDict index 7dd6cf4e935936f87e6b20299917dfb6b6329106..86e5b98f1c1cc0d2e7e9187458354428d176fbe0 100644 --- a/applications/utilities/mesh/manipulation/topoSet/topoSetDict +++ b/applications/utilities/mesh/manipulation/topoSet/topoSetDict @@ -319,7 +319,7 @@ FoamFile // source setToFaceZone; // sourceInfo // { -// set f0; // name of faceSet +// faceSet f0; // name of faceSet // } // // // Select based on faceSet, using cellSet to determine orientation diff --git a/src/OSspecific/POSIX/signals/sigStopAtWriteNow.C b/src/OSspecific/POSIX/signals/sigStopAtWriteNow.C index 09d862ca5d78199aaa08545959edaf0a737a9509..e973bc1bd8fc90039bb43542cf5b3a7dce8000b3 100644 --- a/src/OSspecific/POSIX/signals/sigStopAtWriteNow.C +++ b/src/OSspecific/POSIX/signals/sigStopAtWriteNow.C @@ -83,6 +83,21 @@ Foam::sigStopAtWriteNow::sigStopAtWriteNow { if (signal_ > 0) { + // Check that the signal is different from the writeNowSignal + if (sigWriteNow::signal_ == signal_) + { + FatalErrorIn + ( + "Foam::sigStopAtWriteNow::sigStopAtWriteNow" + "(const bool, const Time&)" + ) << "stopAtWriteNowSignal : " << signal_ + << " cannot be the same as the writeNowSignal." + << " Please change this in the controlDict (" + << findEtcFile("controlDict", false) << ")." + << exit(FatalError); + } + + // Store runTime runTimePtr_ = &runTime; diff --git a/src/OSspecific/POSIX/signals/sigWriteNow.H b/src/OSspecific/POSIX/signals/sigWriteNow.H index 871b980edd99da64f7c44583692b8685acf21aed..477bae825f4b2a48ffd06d3b9203b7af57f40896 100644 --- a/src/OSspecific/POSIX/signals/sigWriteNow.H +++ b/src/OSspecific/POSIX/signals/sigWriteNow.H @@ -67,6 +67,8 @@ class sigWriteNow public: + friend class sigStopAtWriteNow; + // Constructors //- Construct null diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C index 9800ba6006d28f4f40ed9cc9418d78194a79c67d..d1b4926e763910f2f24449d1d4ba6606734354cf 100644 --- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C +++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C @@ -46,6 +46,20 @@ defineTypeNameAndDebug(Foam::globalMeshData, 0); // Geometric matching tolerance. Factor of mesh bounding box. const Foam::scalar Foam::globalMeshData::matchTol_ = 1E-8; +namespace Foam +{ +template<> +class minEqOp<labelPair> +{ +public: + void operator()(labelPair& x, const labelPair& y) const + { + x[0] = min(x[0], y[0]); + x[1] = min(x[1], y[1]); + } +}; +} + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -1063,6 +1077,128 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const } +void Foam::globalMeshData::calcGlobalEdgeOrientation() const +{ + if (debug) + { + Pout<< "globalMeshData::calcGlobalEdgeOrientation() :" + << " calculating edge orientation w.r.t. master edge." << endl; + } + + const globalIndex& globalPoints = globalPointNumbering(); + + // 1. Determine master point + labelList masterPoint; + { + const mapDistribute& map = globalPointSlavesMap(); + + masterPoint.setSize(map.constructSize()); + masterPoint = labelMax; + + for (label pointI = 0; pointI < coupledPatch().nPoints(); pointI++) + { + masterPoint[pointI] = globalPoints.toGlobal(pointI); + } + syncData + ( + masterPoint, + globalPointSlaves(), + globalPointTransformedSlaves(), + map, + minEqOp<label>() + ); + } + + // Now all points should know who is master by comparing their global + // pointID with the masterPointID. We now can use this information + // to find the orientation of the master edge. + + { + const mapDistribute& map = globalEdgeSlavesMap(); + const labelListList& slaves = globalEdgeSlaves(); + const labelListList& transformedSlaves = globalEdgeTransformedSlaves(); + + // Distribute orientation of master edge (in masterPoint numbering) + labelPairList masterEdgeVerts(map.constructSize()); + masterEdgeVerts = labelPair(labelMax, labelMax); + + for (label edgeI = 0; edgeI < coupledPatch().nEdges(); edgeI++) + { + if + ( + ( + slaves[edgeI].size() + + transformedSlaves[edgeI].size() + ) + > 0 + ) + { + // I am master. Fill in my masterPoint equivalent. + + const edge& e = coupledPatch().edges()[edgeI]; + masterEdgeVerts[edgeI] = labelPair + ( + masterPoint[e[0]], + masterPoint[e[1]] + ); + } + } + syncData + ( + masterEdgeVerts, + slaves, + transformedSlaves, + map, + minEqOp<labelPair>() + ); + + // Now check my edges on how they relate to the master's edgeVerts + globalEdgeOrientationPtr_.reset + ( + new PackedBoolList(coupledPatch().nEdges()) + ); + PackedBoolList& globalEdgeOrientation = globalEdgeOrientationPtr_(); + + forAll(coupledPatch().edges(), edgeI) + { + const edge& e = coupledPatch().edges()[edgeI]; + const labelPair masterE + ( + masterPoint[e[0]], + masterPoint[e[1]] + ); + + label stat = labelPair::compare + ( + masterE, + masterEdgeVerts[edgeI] + ); + if (stat == 0) + { + FatalErrorIn + ( + "globalMeshData::calcGlobalEdgeOrientation() const" + ) << "problem : my edge:" << e + << " in master points:" << masterE + << " v.s. masterEdgeVerts:" << masterEdgeVerts[edgeI] + << exit(FatalError); + } + else + { + globalEdgeOrientation[edgeI] = (stat == 1); + } + } + } + + if (debug) + { + Pout<< "globalMeshData::calcGlobalEdgeOrientation() :" + << " finished calculating edge orientation." + << endl; + } +} + + // Calculate uncoupled boundary faces (without calculating // primitiveMesh::pointFaces()) void Foam::globalMeshData::calcPointBoundaryFaces @@ -1660,6 +1796,7 @@ void Foam::globalMeshData::clearOut() globalEdgeNumberingPtr_.clear(); globalEdgeSlavesPtr_.clear(); globalEdgeTransformedSlavesPtr_.clear(); + globalEdgeOrientationPtr_.clear(); globalEdgeSlavesMapPtr_.clear(); // Face @@ -2095,6 +2232,16 @@ const } +const Foam::PackedBoolList& Foam::globalMeshData::globalEdgeOrientation() const +{ + if (!globalEdgeOrientationPtr_.valid()) + { + calcGlobalEdgeOrientation(); + } + return globalEdgeOrientationPtr_(); +} + + const Foam::mapDistribute& Foam::globalMeshData::globalEdgeSlavesMap() const { if (!globalEdgeSlavesMapPtr_.valid()) diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H index 3fca52ae9761cca1a5e5be938b45e8643234b99b..a3530dc86c39700c7edc7da10dfc54e9f4ef0d36 100644 --- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H +++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H @@ -95,6 +95,7 @@ class mapDistribute; template<class T> class EdgeMap; class globalIndex; class globalIndexAndTransform; +class PackedBoolList; /*---------------------------------------------------------------------------*\ Class globalMeshData Declaration @@ -191,6 +192,7 @@ class globalMeshData mutable autoPtr<globalIndex> globalEdgeNumberingPtr_; mutable autoPtr<labelListList> globalEdgeSlavesPtr_; mutable autoPtr<labelListList> globalEdgeTransformedSlavesPtr_; + mutable autoPtr<PackedBoolList> globalEdgeOrientationPtr_; mutable autoPtr<mapDistribute> globalEdgeSlavesMapPtr_; @@ -297,6 +299,9 @@ class globalMeshData //- Calculate global edge addressing. void calcGlobalEdgeSlaves() const; + //- Calculate orientation w.r.t. edge master. + void calcGlobalEdgeOrientation() const; + // Global boundary face/cell addressing @@ -539,6 +544,8 @@ public: const labelListList& globalEdgeSlaves() const; const labelListList& globalEdgeTransformedSlaves() const; const mapDistribute& globalEdgeSlavesMap() const; + //- Is my edge same orientation master edge + const PackedBoolList& globalEdgeOrientation() const; // Coupled point to boundary faces. These are uncoupled boundary // faces only but include empty patches. diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.C index 4a40eac11aa598dab0315f22b4203a2e083283d1..793411c162ab2efa6450b4c9c0a739d175e394a5 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.C +++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.C @@ -30,6 +30,7 @@ License #include "PatchToolsSearch.C" #include "PatchToolsSortEdges.C" #include "PatchToolsNormals.C" +#include "PatchToolsMatch.C" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H index 80e50333bdfc95c78326dcf8e1c2e273be7dffa7..e4ade383ef366829d8b701579a8da83bc32fd4c2 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H +++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H @@ -36,6 +36,7 @@ SourceFiles PatchToolsSearch.C PatchToolsSortEdges.C PatchToolsNormals.C + PatchToolsMatch.C \*---------------------------------------------------------------------------*/ @@ -51,6 +52,7 @@ namespace Foam { class polyMesh; +class PackedBoolList; /*---------------------------------------------------------------------------*\ Class PatchTools Declaration @@ -169,6 +171,55 @@ public: ); + //- Find corresponding points on patches sharing the same points + // p1PointLabels : points on p1 that were matched + // p2PointLabels : corresponding points on p2 + template + < + class Face1, + template<class> class FaceList1, + class PointField1, + class PointType1, + class Face2, + template<class> class FaceList2, + class PointField2, + class PointType2 + > + static void matchPoints + ( + const PrimitivePatch<Face1, FaceList1, PointField1, PointType1>& p1, + const PrimitivePatch<Face2, FaceList2, PointField2, PointType2>& p2, + + labelList& p1PointLabels, + labelList& p2PointLabels + ); + + //- Find corresponding edges on patches sharing the same points + // p1EdgeLabels : edges on p1 that were matched + // p2EdgeLabels : corresponding edges on p2 + // sameOrientation : same orientation? + template + < + class Face1, + template<class> class FaceList1, + class PointField1, + class PointType1, + class Face2, + template<class> class FaceList2, + class PointField2, + class PointType2 + > + static void matchEdges + ( + const PrimitivePatch<Face1, FaceList1, PointField1, PointType1>& p1, + const PrimitivePatch<Face2, FaceList2, PointField2, PointType2>& p2, + + labelList& p1EdgeLabels, + labelList& p2EdgeLabels, + PackedBoolList& sameOrientation + ); + + //- Return parallel consistent point normals for patches (on boundary faces) // using mesh points. template diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsMatch.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsMatch.C new file mode 100644 index 0000000000000000000000000000000000000000..1af68e8991236c2f9d88f50826243e7adbef6853 --- /dev/null +++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsMatch.C @@ -0,0 +1,137 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "PatchTools.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +< + class Face1, + template<class> class FaceList1, + class PointField1, + class PointType1, + class Face2, + template<class> class FaceList2, + class PointField2, + class PointType2 +> +void Foam::PatchTools::matchPoints +( + const PrimitivePatch<Face1, FaceList1, PointField1, PointType1>& p1, + const PrimitivePatch<Face2, FaceList2, PointField2, PointType2>& p2, + + labelList& p1PointLabels, + labelList& p2PointLabels +) +{ + p1PointLabels.setSize(p1.nPoints()); + p2PointLabels.setSize(p1.nPoints()); + + label nMatches = 0; + + forAll(p1.meshPoints(), pointI) + { + label meshPointI = p1.meshPoints()[pointI]; + + Map<label>::const_iterator iter = p2.meshPointMap().find + ( + meshPointI + ); + + if (iter != p2.meshPointMap().end()) + { + p1PointLabels[nMatches] = pointI; + p2PointLabels[nMatches] = iter(); + nMatches++; + } + } + p1PointLabels.setSize(nMatches); + p2PointLabels.setSize(nMatches); +} + + +template +< + class Face1, + template<class> class FaceList1, + class PointField1, + class PointType1, + class Face2, + template<class> class FaceList2, + class PointField2, + class PointType2 +> +void Foam::PatchTools::matchEdges +( + const PrimitivePatch<Face1, FaceList1, PointField1, PointType1>& p1, + const PrimitivePatch<Face2, FaceList2, PointField2, PointType2>& p2, + + labelList& p1EdgeLabels, + labelList& p2EdgeLabels, + PackedBoolList& sameOrientation +) +{ + p1EdgeLabels.setSize(p1.nEdges()); + p2EdgeLabels.setSize(p1.nEdges()); + sameOrientation.setSize(p1.nEdges()); + sameOrientation = 0; + + label nMatches = 0; + + EdgeMap<label> edgeToIndex(2*p1.nEdges()); + forAll(p1.edges(), edgeI) + { + const edge& e = p1.edges()[edgeI]; + const edge meshE + ( + p1.meshPoints()[e[0]], + p1.meshPoints()[e[1]] + ); + edgeToIndex.insert(meshE, edgeI); + } + + forAll(p2.edges(), edgeI) + { + const edge& e = p2.edges()[edgeI]; + const edge meshE(p2.meshPoints()[e[0]], p2.meshPoints()[e[1]]); + + EdgeMap<label>::const_iterator iter = edgeToIndex.find(meshE); + + if (iter != edgeToIndex.end()) + { + p1EdgeLabels[nMatches] = iter(); + p2EdgeLabels[nMatches] = edgeI; + sameOrientation[nMatches] = (meshE[0] == iter.key()[0]); + nMatches++; + } + } + p1EdgeLabels.setSize(nMatches); + p2EdgeLabels.setSize(nMatches); + sameOrientation.setSize(nMatches); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 3bcd7b80d55dfec16641962a6dbf6ff0e12c12fd..cc2f3e4225a0e329c88e772acb0748631e5042a7 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -390,6 +390,12 @@ $(basicSource)/actuationDiskSource/actuationDiskSource.C $(basicSource)/radialActuationDiskSource/radialActuationDiskSource.C $(basicSource)/explicitSource/explicitSource.C $(basicSource)/explicitSetValue/explicitSetValue.C +$(basicSource)/rotorDiskSource/rotorDiskSource.C +$(basicSource)/rotorDiskSource/bladeModel/bladeModel.C +$(basicSource)/rotorDiskSource/profileModel/profileModel.C +$(basicSource)/rotorDiskSource/profileModel/profileModelList.C +$(basicSource)/rotorDiskSource/profileModel/lookup/lookupProfile.C +$(basicSource)/rotorDiskSource/profileModel/series/seriesProfile.C LIB = $(FOAM_LIBBIN)/libfiniteVolume diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.C index 32a025d993bd4010f8c5a3e4aee8656ce9d04b80..e7e996d1db9784c8a07821fdf425c03c25cca17f 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.C @@ -75,11 +75,11 @@ Foam::actuationDiskSource::actuationDiskSource ) : basicSource(name, modelType, dict, mesh), - dict_(dict.subDict(modelType + "Coeffs")), - diskDir_(dict_.lookup("diskDir")), - Cp_(readScalar(dict_.lookup("Cp"))), - Ct_(readScalar(dict_.lookup("Ct"))), - diskArea_(readScalar(dict_.lookup("diskArea"))) + coeffs_(dict.subDict(modelType + "Coeffs")), + diskDir_(coeffs_.lookup("diskDir")), + Cp_(readScalar(coeffs_.lookup("Cp"))), + Ct_(readScalar(coeffs_.lookup("Ct"))), + diskArea_(readScalar(coeffs_.lookup("diskArea"))) { Info<< " - creating actuation disk zone: " << this->name() << endl; @@ -132,20 +132,8 @@ void Foam::actuationDiskSource::addSu(fvMatrix<vector>& UEqn) void Foam::actuationDiskSource::writeData(Ostream& os) const { - os << indent << token::BEGIN_BLOCK << incrIndent << nl; - os.writeKeyword("name") << this->name() << token::END_STATEMENT << nl; - - if (dict_.found("note")) - { - os.writeKeyword("note") << string(dict_.lookup("note")) - << token::END_STATEMENT << nl; - } - - os << indent << "actuationDisk"; - + os << indent << name_ << endl; dict_.write(os); - - os << decrIndent << indent << token::END_BLOCK << endl; } @@ -153,13 +141,11 @@ bool Foam::actuationDiskSource::read(const dictionary& dict) { if (basicSource::read(dict)) { - const dictionary& sourceDict = dict.subDict(name()); - const dictionary& subDictCoeffs = - sourceDict.subDict(typeName + "Coeffs"); - subDictCoeffs.readIfPresent("diskDir", diskDir_); - subDictCoeffs.readIfPresent("Cp", Cp_); - subDictCoeffs.readIfPresent("Ct", Ct_); - subDictCoeffs.readIfPresent("diskArea", diskArea_); + coeffs_ = dict.subDict(typeName + "Coeffs"); + coeffs_.readIfPresent("diskDir", diskDir_); + coeffs_.readIfPresent("Cp", Cp_); + coeffs_.readIfPresent("Ct", Ct_); + coeffs_.readIfPresent("diskArea", diskArea_); checkData(); diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.H index 8b54ac08710686b5f43e3555abbd07685f62ddb7..cbe98ea4c8d24e2489a9691c24931256bbe11138 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.H +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.H @@ -74,8 +74,8 @@ protected: // Protected data - //- Sub dictionary with actuationDisk information - const dictionary& dict_; + //- Coefficients dictionary + dictionary coeffs_; //- Disk area normal vector diskDir_; diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C index 114e6e58cbe692549e742f4d515bf2cbbe962c74..cd7394a173bc9029c9517e2e411b32113868b262 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C @@ -277,28 +277,20 @@ bool Foam::basicSource::isActive() void Foam::basicSource::addSu(Foam::fvMatrix<vector>& Eqn) { - notImplemented - ( - "Foam::basicSource addSu(Foam::fvMatrix<vector>& Eqn)" - ); + notImplemented("Foam::basicSource addSu(Foam::fvMatrix<vector>& Eqn)"); } void Foam::basicSource::addSu(Foam::fvMatrix<scalar>& Eqn) { - notImplemented - ( - "Foam::basicSource addSu(Foam::fvMatrix<scalar>& Eqn)" - ); + notImplemented("Foam::basicSource addSu(Foam::fvMatrix<scalar>& Eqn)"); } void Foam::basicSource::setValue(Foam::fvMatrix<scalar>& Eqn) { - notImplemented - ( - "Foam::basicSource setValue(Foam::fvMatrix<scalar>& Eqn)" - ); + notImplemented("Foam::basicSource setValue(Foam::fvMatrix<scalar>& Eqn)"); } + // ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.H index 43c9e8bd47028e8fa32b05774ffdb1946d0be1f7..9862282479b3c0db10e3c4857353c9335f003540 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.H +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.H @@ -202,10 +202,7 @@ public: //- Return clone autoPtr<basicSource> clone() const { - notImplemented - ( - "autoPtr<basicSource> clone() const" - ); + notImplemented("autoPtr<basicSource> clone() const"); return autoPtr<basicSource>(NULL); } @@ -273,7 +270,7 @@ public: //- Return const access to the mesh database inline const fvMesh& mesh() const; - //- Return dictionay + //- Return dictionary inline const dictionary& dictCoeffs() const; //- Return const access to the source active flag @@ -340,7 +337,6 @@ public: //- Read source dictionary virtual bool read(const dictionary& dict) = 0; - }; diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceIO.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceIO.C index 8b425cf91c184669cdc76cc2791c2de33b21e25c..caec8b58cf78769260e7f08816fe352a2f77826c 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceIO.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceIO.C @@ -77,10 +77,9 @@ void Foam::basicSource::writeData(Ostream& os) const bool Foam::basicSource::read(const dictionary& dict) { - const dictionary& sourceDict = dict.subDict(name_); - active_ = readBool(sourceDict.lookup("active")); - timeStart_ = readScalar(sourceDict.lookup("timeStart")); - duration_ = readScalar(sourceDict.lookup("duration")); + active_ = readBool(dict.lookup("active")); + timeStart_ = readScalar(dict.lookup("timeStart")); + duration_ = readScalar(dict.lookup("duration")); return true; } diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C index e9827409c5aef81c4de1449e5ae03aa1cddbb58e..79af0251430efbafbb62641764fff90a7e7fdec3 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C @@ -51,12 +51,12 @@ Foam::basicSourceList::basicSourceList forAllConstIter(dictionary, dict, iter) { const word& name = iter().keyword(); - const dictionary& dict = iter().dict(); + const dictionary& sourceDict = iter().dict(); this->set ( i++, - basicSource::New(name, dict, mesh) + basicSource::New(name, sourceDict, mesh) ); } } @@ -108,7 +108,8 @@ bool Foam::basicSourceList::read(const dictionary& dict) bool allOk = true; forAll(*this, i) { - bool ok = this->operator[](i).read(dict); + basicSource& bs = this->operator[](i); + bool ok = bs.read(dict.subDict(bs.name())); allOk = (allOk && ok); } return allOk; @@ -117,12 +118,6 @@ bool Foam::basicSourceList::read(const dictionary& dict) bool Foam::basicSourceList::writeData(Ostream& os) const { - // Write size of list - os << nl << this->size(); - - // Write beginning of contents - os << nl << token::BEGIN_LIST; - // Write list contents forAll(*this, i) { @@ -130,9 +125,6 @@ bool Foam::basicSourceList::writeData(Ostream& os) const this->operator[](i).writeData(os); } - // Write end of contents - os << token::END_LIST << token::END_STATEMENT << nl; - // Check state of IOstream return os.good(); } diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.C index 2a14b793a8de3db75521adcff5d0a1ec05704b6b..898f2c41c3728ee47628edec240626365212f172 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.C @@ -72,11 +72,9 @@ void Foam::explicitSetValue::setFieldData(const dictionary& dict) } else { - FatalErrorIn - ( - "explicitSetValue::setFieldData" - ) << "header not OK " << io.name() - << exit(FatalError); + FatalErrorIn("explicitSetValue::setFieldData") + << "header not OK for field " << io.name() + << abort(FatalError); } } @@ -96,9 +94,9 @@ Foam::explicitSetValue::explicitSetValue ) : basicSource(name, modelType, dict, mesh), - dict_(dict.subDict(modelType + "Coeffs")) + coeffs_(dict.subDict(modelType + "Coeffs")) { - setFieldData(dict_.subDict("fieldData")); + setFieldData(coeffs_.subDict("fieldData")); } diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.H index 643e6b618243d40c4bb7401e13c7bdcc4e7a1352..df3d404f1ac373cdc1c584b6fcec73d08fbd2649 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.H +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.H @@ -89,8 +89,8 @@ protected: // Protected data - //- Sub dictionary for time activated explicit sources - const dictionary& dict_; + //- Coefficients dictionary + dictionary coeffs_; // Protected functions @@ -119,10 +119,7 @@ public: //- Return clone autoPtr<explicitSetValue> clone() const { - notImplemented - ( - "autoPtr<explicitSetValue> clone() const" - ); + notImplemented("autoPtr<explicitSetValue> clone() const"); return autoPtr<explicitSetValue>(NULL); } diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueIO.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueIO.C index dce0a53ddf69a4776aba28b21bc040ab76824773..24e9205606db82c5af9e1cebfe659628e7b2b6a6 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueIO.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueIO.C @@ -29,22 +29,8 @@ License void Foam::explicitSetValue::writeData(Ostream& os) const { - os << indent << name_ << nl - << indent << token::BEGIN_BLOCK << incrIndent << nl; - - if (scalarFields_.size() > 0) - { - os.writeKeyword("scalarFields") << scalarFields_ - << token::END_STATEMENT << nl; - } - - if (vectorFields_.size() > 0) - { - os.writeKeyword("vectorFields") << vectorFields_ - << token::END_STATEMENT << nl; - } - - os << decrIndent << indent << token::END_BLOCK << endl; + os << indent << name_ << endl; + dict_.write(os); } @@ -52,12 +38,8 @@ bool Foam::explicitSetValue::read(const dictionary& dict) { if (basicSource::read(dict)) { - const dictionary& sourceDict = dict.subDict(name()); - const dictionary& subDictCoeffs = sourceDict.subDict - ( - typeName + "Coeffs" - ); - setFieldData(subDictCoeffs.subDict("fieldData")); + coeffs_ = dict.subDict(typeName + "Coeffs"); + setFieldData(coeffs_.subDict("fieldData")); return true; } else diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.C index de01ebb5309c7e2bd4975a1ac4013ae3a6336350..af5555f3fb09af2f0c73a427a8930c8c4c10ff00 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.C @@ -112,8 +112,8 @@ Foam::explicitSource::explicitSource ) : basicSource(name, modelType, dict, mesh), - dict_(dict.subDict(modelType + "Coeffs")), - volumeMode_(volumeModeTypeNames_.read(dict_.lookup("volumeMode"))) + coeffs_(dict.subDict(modelType + "Coeffs")), + volumeMode_(volumeModeTypeNames_.read(coeffs_.lookup("volumeMode"))) { setFieldData(dict_.subDict("fieldData")); } diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.H index 9b88fac753e69c088d0b920d2d2928eb9681e772..05c3782254769ad49c32318e781dd132ab7c19b7 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.H +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.H @@ -110,8 +110,8 @@ protected: // Protected data - //- Sub dictionary for time activated explicit sources - const dictionary& dict_; + //- Coefficients dictionary + dictionary coeffs_; //- Volume mode volumeModeType volumeMode_; diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSourceIO.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSourceIO.C index 3b7b1b7c3521017e990f6fa591e59969cc27e507..f97ad329067091f0ba063410bc818d8e06e37e9c 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSourceIO.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSourceIO.C @@ -29,25 +29,8 @@ License void Foam::explicitSource::writeData(Ostream& os) const { - os << indent << name_ << nl - << indent << token::BEGIN_BLOCK << incrIndent << nl; - - os.writeKeyword("volumeMode") << volumeModeTypeNames_[volumeMode_] - << token::END_STATEMENT << nl; - - if (scalarFields_.size() > 0) - { - os.writeKeyword("scalarFields") << scalarFields_ - << token::END_STATEMENT << nl; - } - - if (vectorFields_.size() > 0) - { - os.writeKeyword("vectorFields") << vectorFields_ - << token::END_STATEMENT << nl; - } - - os << decrIndent << indent << token::END_BLOCK << endl; + os << indent << name_ << endl; + dict_.write(os); } @@ -55,12 +38,8 @@ bool Foam::explicitSource::read(const dictionary& dict) { if (basicSource::read(dict)) { - const dictionary& sourceDict = dict.subDict(name()); - const dictionary& subDictCoeffs = sourceDict.subDict - ( - typeName + "Coeffs" - ); - setFieldData(subDictCoeffs.subDict("fieldData")); + coeffs_ = dict.subDict(typeName + "Coeffs"); + setFieldData(coeffs_.subDict("fieldData")); return true; } else diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.C index e01a2a57361dbe299df301675f59953669972323..0444a9408b3dc3e81896634da6e641dbff45dfae 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.C @@ -53,10 +53,10 @@ Foam::radialActuationDiskSource::radialActuationDiskSource ) : actuationDiskSource(name, modelType, dict, mesh), - dict_(dict.subDict(modelType + "Coeffs")), + coeffsDict_(dict.subDict(modelType + "Coeffs")), coeffs_() { - dict_.lookup("coeffs") >> coeffs_; + coeffsDict_.lookup("coeffs") >> coeffs_; Info<< " - creating radial actuation disk zone: " << this->name() << endl; } @@ -114,14 +114,12 @@ bool Foam::radialActuationDiskSource::read(const dictionary& dict) { if (basicSource::read(dict)) { - const dictionary& sourceDict = dict.subDict(name()); - const dictionary& subDictCoeffs = - sourceDict.subDict(typeName + "Coeffs"); - subDictCoeffs.readIfPresent("diskDir", diskDir_); - subDictCoeffs.readIfPresent("Cp", Cp_); - subDictCoeffs.readIfPresent("Ct", Ct_); - subDictCoeffs.readIfPresent("diskArea", diskArea_); - subDictCoeffs.lookup("coeffs") >> coeffs_; + const dictionary& coeffsDict_ = dict.subDict(typeName + "Coeffs"); + coeffsDict_.readIfPresent("diskDir", diskDir_); + coeffsDict_.readIfPresent("Cp", Cp_); + coeffsDict_.readIfPresent("Ct", Ct_); + coeffsDict_.readIfPresent("diskArea", diskArea_); + coeffsDict_.lookup("coeffs") >> coeffs_; return true; } else diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.H index 701e752ff400ccb1650a70c5547cff7c05c75da0..245f002fa51446d0035f74c58987764ec3640783 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.H +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.H @@ -62,7 +62,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class radialActuationDiskSource Declaration + Class radialActuationDiskSource Declaration \*---------------------------------------------------------------------------*/ class radialActuationDiskSource @@ -71,8 +71,8 @@ class radialActuationDiskSource { // Private data - //- Sub dictionary with model information - const dictionary& dict_; + //- Coefficients dictionary + dictionary coeffsDict_; //- Coeffcients for the radial distribution FixedList<scalar, 3> coeffs_; @@ -123,7 +123,7 @@ public: // Public Functions - //-Source term to fvMatrix<vector> + //- Source term to fvMatrix<vector> virtual void addSu(fvMatrix<vector>& UEqn); diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C new file mode 100644 index 0000000000000000000000000000000000000000..a1f76d3e99344fb74acddb2f4e84d235b8449424 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C @@ -0,0 +1,195 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "bladeModel.H" +#include "unitConversion.H" +#include "Tuple2.H" +#include "vector.H" +#include "IFstream.H" + + +// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // + +bool Foam::bladeModel::readFromFile() const +{ + return fName_ != fileName::null; +} + + +void Foam::bladeModel::interpolateWeights +( + const scalar& xIn, + const List<scalar>& values, + label& i1, + label& i2, + scalar& ddx +) const +{ + scalar x = -GREAT; + label nElem = values.size(); + + i2 = 0; + while ((x < xIn) && (i2 < nElem)) + { + x = values[i2]; + i2++; + } + + if (i2 == 0) + { + i1 = i2; + ddx = 0.0; + return; + } + else if (i2 == values.size()) + { + i2 = values.size() - 1; + i1 = i2; + ddx = 0.0; + return; + } + else + { + i1 = i2 - 1; + ddx = (xIn - values[i1])/(values[i2] - values[i1]); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::bladeModel::bladeModel(const dictionary& dict) +: + profileName_(), + profileID_(), + radius_(), + twist_(), + chord_(), + fName_(fileName::null) +{ + List<Tuple2<word, vector> > data; + if (readFromFile()) + { + IFstream is(fName_); + is >> data; + } + else + { + dict.lookup("data") >> data; + } + + + if (data.size() > 0) + { + profileName_.setSize(data.size()); + profileID_.setSize(data.size()); + radius_.setSize(data.size()); + twist_.setSize(data.size()); + chord_.setSize(data.size()); + + forAll(data, i) + { + profileName_[i] = data[i].first(); + profileID_[i] = -1; + radius_[i] = data[i].second()[0]; + twist_[i] = degToRad(data[i].second()[1]); + chord_[i] = data[i].second()[2]; + } + } + else + { + FatalErrorIn + ( + "Foam::bladeModel::bladeModel" + "(" + "const dictionary&, " + "const word&" + ")" + ) << "No blade data specified" << exit(FatalError); + } +} + +// * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * // + +Foam::bladeModel::~bladeModel() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +const Foam::List<Foam::word>& Foam::bladeModel::profileName() const +{ + return profileName_; +} + + +const Foam::List<Foam::label>& Foam::bladeModel::profileID() const +{ + return profileID_; +} + + +const Foam::List<Foam::scalar>& Foam::bladeModel::radius() const +{ + return radius_; +} + + +const Foam::List<Foam::scalar>& Foam::bladeModel::twist() const +{ + return twist_; +} + + +const Foam::List<Foam::scalar>& Foam::bladeModel::chord() const +{ + return chord_; +} + + +Foam::List<Foam::label>& Foam::bladeModel::profileID() +{ + return profileID_; +} + + +void Foam::bladeModel::interpolate +( + const scalar radius, + scalar& twist, + scalar& chord, + label& i1, + label& i2, + scalar& invDr +) const +{ + interpolateWeights(radius, radius_, i1, i2, invDr); + + twist = invDr*(twist_[i2] - twist_[i1]) + twist_[i1]; + chord = invDr*(chord_[i2] - chord_[i1]) + chord_[i1]; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.H new file mode 100644 index 0000000000000000000000000000000000000000..df26e7bd401a5b1bdf617ec0620dc25ccac1ba4c --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.H @@ -0,0 +1,151 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Class + Foam::bladeModel + +Description + Blade model class + +SourceFiles + bladeModel.C + +\*---------------------------------------------------------------------------*/ + +#ifndef bladeModel_H +#define bladeModel_H + +#include "List.H" +#include "dictionary.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class bladeModel Declaration +\*---------------------------------------------------------------------------*/ + +class bladeModel +{ + +protected: + + // Protected data + + //- Corresponding profile name per section + List<word> profileName_; + + //- Corresponding profile ID per section + List<label> profileID_; + + //- Radius [m] + List<scalar> radius_; + + //- Twist [deg] on input, converted to [rad] + List<scalar> twist_; + + //- Chord [m] + List<scalar> chord_; + + //- File name (optional) + fileName fName_; + + + // Protected Member Functions + + //- Return ture if file name is set + bool readFromFile() const; + + //- Return the interpolation indices and gradient + void interpolateWeights + ( + const scalar& xIn, + const List<scalar>& values, + label& i1, + label& i2, + scalar& ddx + ) const; + + +public: + + //- Constructor + bladeModel(const dictionary& dict); + + + //- Destructor + virtual ~bladeModel(); + + + // Member functions + + // Access + + //- Return const access to the profile name list + const List<word>& profileName() const; + + //- Return const access to the profile ID list + const List<label>& profileID() const; + + //- Return const access to the radius list + const List<scalar>& radius() const; + + //- Return const access to the twist list + const List<scalar>& twist() const; + + //- Return const access to the chord list + const List<scalar>& chord() const; + + + // Edit + + //- Return non-const access to the profile ID list + List<label>& profileID(); + + + // Evaluation + + //- Return the twist and chord for a given radius + virtual void interpolate + ( + const scalar radius, + scalar& twist, + scalar& chord, + label& i1, + label& i2, + scalar& invDr + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C new file mode 100644 index 0000000000000000000000000000000000000000..ae80d03995d4e11bd5bd17043650240e70f31876 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C @@ -0,0 +1,148 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "lookupProfile.H" +#include "addToRunTimeSelectionTable.H" +#include "vector.H" +#include "unitConversion.H" +#include "IFstream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(lookupProfile, 0); + addToRunTimeSelectionTable(profileModel, lookupProfile, dictionary); +} + + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +void Foam::lookupProfile::interpolateWeights +( + const scalar& xIn, + const List<scalar>& values, + label& i1, + label& i2, + scalar& ddx +) const +{ + scalar x = -GREAT; + label nElem = values.size(); + + i2 = 0; + while ((x < xIn) && (i2 < nElem)) + { + x = values[i2]; + i2++; + } + + if (i2 == 0) + { + i1 = i2; + ddx = 0.0; + return; + } + else if (i2 == values.size()) + { + i2 = values.size() - 1; + i1 = i2; + ddx = 0.0; + return; + } + else + { + i1 = i2 - 1; + ddx = (xIn - values[i1])/(values[i2] - values[i1]); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::lookupProfile::lookupProfile +( + const dictionary& dict, + const word& modelName +) +: + profileModel(dict, modelName), + AOA_(), + Cd_(), + Cl_() +{ + List<vector> data; + if (readFromFile()) + { + IFstream is(fName_); + is >> data; + } + else + { + dict.lookup("data") >> data; + } + + if (data.size() > 0) + { + AOA_.setSize(data.size()); + Cd_.setSize(data.size()); + Cl_.setSize(data.size()); + + forAll(data, i) + { + AOA_[i] = degToRad(data[i][0]); + Cd_[i] = data[i][1]; + Cl_[i] = data[i][2]; + } + } + else + { + FatalErrorIn + ( + "Foam::lookupProfile::lookupProfile" + "(" + "const dictionary&, " + "const word&" + ")" + ) << "No profile data specified" << exit(FatalError); + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::lookupProfile::Cdl(const scalar alpha, scalar& Cd, scalar& Cl) const +{ + label i1 = -1; + label i2 = -1; + scalar invAlpha = -1.0; + interpolateWeights(alpha, AOA_, i1, i2, invAlpha); + + Cd = invAlpha*(Cd_[i2] - Cd_[i1]) + Cd_[i1]; + Cl = invAlpha*(Cl_[i2] - Cl_[i1]) + Cl_[i1]; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.H new file mode 100644 index 0000000000000000000000000000000000000000..deb5c1f725ab8096c9a0471e09d8784fba7e6061 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.H @@ -0,0 +1,124 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Class + Foam::lookupProfile + +Description + Look-up based profile data - drag and lift coefficients are lineraly + interpolated based on the supplied angle of attack + + Input in list format: + + data + ( + (AOA1 Cd1 Cl2) + (AOA2 Cd2 Cl2) + ... + (AOAN CdN CdN) + ); + + where: + AOA = angle of attack [deg] converted to [rad] internally + Cd = drag coefficient + Cl = lift coefficient + +SourceFiles + lookupProfile.C + +\*---------------------------------------------------------------------------*/ + +#ifndef lookupProfile_H +#define lookupProfile_H + +#include "profileModel.H" +#include "List.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class lookupProfile Declaration +\*---------------------------------------------------------------------------*/ + +class lookupProfile +: + public profileModel +{ + +protected: + + // Protected data + + //- List of angle-of-attack values [deg] on input, converted to [rad] + List<scalar> AOA_; + + //- List of drag coefficient values + List<scalar> Cd_; + + //- List of lift coefficient values + List<scalar> Cl_; + + + // Protected Member Functions + + //- Return the interpolation indices and gradient + void interpolateWeights + ( + const scalar& xIn, + const List<scalar>& values, + label& i1, + label& i2, + scalar& ddx + ) const; + + +public: + + //- Runtime type information + TypeName("lookup"); + + //- Constructor + lookupProfile(const dictionary& dict, const word& modelName); + + + // Member functions + + // Evaluation + + //- Return the Cd and Cl for a given angle-of-attack + virtual void Cdl(const scalar alpha, scalar& Cd, scalar& Cl) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModel.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModel.C new file mode 100644 index 0000000000000000000000000000000000000000..c9008d276c644533d6a446f453eb20b3d709fafe --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModel.C @@ -0,0 +1,99 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "profileModel.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(profileModel, 0); + defineRunTimeSelectionTable(profileModel, dictionary); +} + + +// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // + +bool Foam::profileModel::readFromFile() const +{ + return fName_ != fileName::null; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::profileModel::profileModel(const dictionary& dict, const word& name) +: + dict_(dict), + name_(name), + fName_(fileName::null) +{ + dict.readIfPresent("fileName", fName_); +} + +// * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * // + +Foam::profileModel::~profileModel() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +const Foam::word& Foam::profileModel::name() const +{ + return name_; +} + + +Foam::autoPtr<Foam::profileModel> Foam::profileModel::New +( + const dictionary& dict +) +{ + const word& modelName(dict.dictName()); + + const word modelType(dict.lookup("type")); + + Info<< " - creating " << modelType << " profile " << modelName << endl; + + dictionaryConstructorTable::iterator cstrIter = + dictionaryConstructorTablePtr_->find(modelType); + + if (cstrIter == dictionaryConstructorTablePtr_->end()) + { + FatalErrorIn("profileModel::New(const dictionary&)") + << "Unknown profile model type " << modelType + << nl << nl + << "Valid model types are :" << nl + << dictionaryConstructorTablePtr_->sortedToc() + << exit(FatalError); + } + + return autoPtr<profileModel>(cstrIter()(dict, modelName)); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModel.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModel.H new file mode 100644 index 0000000000000000000000000000000000000000..5c850ec1bb30a214b71b35586e0c09b8af57039f --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModel.H @@ -0,0 +1,136 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Class + Foam::profileModel + +Description + Base class for profile models + +SourceFiles + profileModel.C + +\*---------------------------------------------------------------------------*/ + +#ifndef profileModel_H +#define profileModel_H + +#include "autoPtr.H" +#include "runTimeSelectionTables.H" +#include "dictionary.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class profileModel Declaration +\*---------------------------------------------------------------------------*/ + +class profileModel +{ + +protected: + + // Protected data + + //- Coefficients dictionary + const dictionary dict_; + + //- Name of profile model + const word name_; + + //- File name (optional) + fileName fName_; + + + // Protected Member Functions + + //- Return ture if file name is set + bool readFromFile() const; + + +public: + + //- Runtime type information + TypeName("profileModel"); + + + // Declare run-time constructor selection table + declareRunTimeSelectionTable + ( + autoPtr, + profileModel, + dictionary, + ( + const dictionary& dict, + const word& modelName + ), + (dict, modelName) + ); + + + // Selectors + + //- Return a reference to the selected basicSource model + static autoPtr<profileModel> New(const dictionary& dict); + + + //- Constructor + profileModel(const dictionary& dict, const word& modelName); + + + //- Destructor + virtual ~profileModel(); + + + // Member functions + + // Access + + //- Return const access to the source name + const word& name() const; + + + // Evaluation + + //- Return the Cd and Cl for a given angle-of-attack + virtual void Cdl + ( + const scalar alpha, + scalar& Cd, + scalar& Cl + ) const = 0; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModelList.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModelList.C new file mode 100644 index 0000000000000000000000000000000000000000..904b01372b540aadf146c3e345744ba0d41298ad --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModelList.C @@ -0,0 +1,121 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "profileModelList.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::profileModelList::profileModelList +( + const dictionary& dict, + const bool readFields +) +: + PtrList<profileModel>(), + dict_(dict) +{ + if (readFields) + { + wordList modelNames(dict.toc()); + + Info<< " Constructing blade profiles:" << endl; + + if (modelNames.size() > 0) + { + this->setSize(modelNames.size()); + + forAll(modelNames, i) + { + const word& modelName = modelNames[i]; + + this->set + ( + i, + profileModel::New(dict.subDict(modelName)) + ); + } + } + else + { + Info<< " none" << endl; + } + } +} + + +// * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * // + +Foam::profileModelList::~profileModelList() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::profileModelList::connectBlades +( + const List<word>& names, + List<label>& addr +) const +{ + // construct the addressing between blade sections and profiles + forAll(names, bI) + { + label index = -1; + const word& profileName = names[bI]; + + forAll(*this, pI) + { + const profileModel& pm = this->operator[](pI); + + if (pm.name() == profileName) + { + index = pI; + break; + } + } + + if (index == -1) + { + List<word> profileNames(size()); + forAll(*this, i) + { + const profileModel& pm = this->operator[](i); + profileNames[i] = pm.name(); + } + + FatalErrorIn("void Foam::connectBlades(List<word>& names) const") + << "Profile " << profileName << " could not be found " + << "in profile list. Available profiles are" + << profileNames << exit(FatalError); + } + else + { + addr[bI] = index; + } + } +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModelList.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModelList.H new file mode 100644 index 0000000000000000000000000000000000000000..74326dced4b74ed907493363c2f066c831c212cc --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/profileModelList.H @@ -0,0 +1,91 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Class + Foam::profileModelList + +Description + Base class for profile models + +SourceFiles + profileModel.C + +\*---------------------------------------------------------------------------*/ + +#ifndef profileModelList_H +#define profileModelList_H + +#include "PtrList.H" +#include "profileModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class profileModelList Declaration +\*---------------------------------------------------------------------------*/ + +class profileModelList +: + public PtrList<profileModel> +{ + +protected: + + // Protected data + + //- Dictionary + const dictionary dict_; + + +public: + + //- Constructor + profileModelList(const dictionary& dict, const bool readFields = true); + + //- Destructor + ~profileModelList(); + + + // Member Functions + + //- Set blade->profile addressing + void connectBlades + ( + const List<word>& names, + List<label>& addr + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.C new file mode 100644 index 0000000000000000000000000000000000000000..865ed3ab62ec7459c76ad59b6e7bb1d417c233f1 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.C @@ -0,0 +1,116 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "seriesProfile.H" +#include "addToRunTimeSelectionTable.H" +#include "IFstream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(seriesProfile, 0); + addToRunTimeSelectionTable(profileModel, seriesProfile, dictionary); +} + + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +Foam::scalar Foam::seriesProfile::evaluate +( + const scalar& xIn, + const List<scalar>& values +) const +{ + scalar result = 0.0; + + forAll(values, i) + { + result += values[i]*cos((i+1)*xIn); + } + + return result; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::seriesProfile::seriesProfile +( + const dictionary& dict, + const word& modelName +) +: + profileModel(dict, modelName), + CdCoeffs_(), + ClCoeffs_() +{ + if (readFromFile()) + { + IFstream is(fName_); + is >> CdCoeffs_ >> ClCoeffs_; + } + else + { + dict.lookup("CdCoeffs") >> CdCoeffs_; + dict.lookup("ClCoeffs") >> ClCoeffs_; + } + + + if (CdCoeffs_.empty()) + { + FatalErrorIn + ( + "Foam::seriesProfile::seriesProfile" + "(" + "const dictionary&, " + "const word&" + ")" + ) << "CdCoeffs must be specified" << exit(FatalError); + } + if (ClCoeffs_.empty()) + { + FatalErrorIn + ( + "Foam::seriesProfile::seriesProfile" + "(" + "const dictionary&, " + "const word&" + ")" + ) << "ClCoeffs must be specified" << exit(FatalError); + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::seriesProfile::Cdl(const scalar alpha, scalar& Cd, scalar& Cl) const +{ + Cd = evaluate(alpha, CdCoeffs_); + Cl = evaluate(alpha, ClCoeffs_); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.H new file mode 100644 index 0000000000000000000000000000000000000000..73b474aae1f3bb22bbdbacc1b29ea35a914051fb --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.H @@ -0,0 +1,116 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Class + Foam::seriesProfile + +Description + Series-up based profile data - drag and lift coefficients computed as + sum of cosine series + + Cd = sum_i(CdCoeff)*cos(i*AOA) + Cl = sum_i(ClCoeff)*cos(i*AOA) + + where: + AOA = angle of attack [deg] converted to [rad] internally + Cd = drag coefficent + Cl = lift coefficent + + Input in two (arbitrary length) lists: + + CdCoeffs (coeff1 coeff2 ... coeffN); + ClCoeffs (coeff1 coeff2 ... coeffN); + +SourceFiles + seriesProfile.C + +\*---------------------------------------------------------------------------*/ + +#ifndef seriesProfile_H +#define seriesProfile_H + +#include "profileModel.H" +#include "List.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class seriesProfile Declaration +\*---------------------------------------------------------------------------*/ + +class seriesProfile +: + public profileModel +{ + +protected: + + // Protected data + + //- List of drag coefficient values + List<scalar> CdCoeffs_; + + //- List of lift coefficient values + List<scalar> ClCoeffs_; + + + // Protected Member Functions + + //- Evaluate + scalar evaluate + ( + const scalar& xIn, + const List<scalar>& values + ) const; + + +public: + + //- Runtime type information + TypeName("series"); + + //- Constructor + seriesProfile(const dictionary& dict, const word& modelName); + + + // Member functions + + // Evaluation + + //- Return the Cd and Cl for a given angle-of-attack + virtual void Cdl(const scalar alpha, scalar& Cd, scalar& Cl) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C new file mode 100644 index 0000000000000000000000000000000000000000..db4fd5d80336bb3c72f418c02e1b18bd0169c32c --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C @@ -0,0 +1,495 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "rotorDiskSource.H" +#include "addToRunTimeSelectionTable.H" +#include "mathematicalConstants.H" +#include "unitConversion.H" +#include "geometricOneField.H" + +using namespace Foam::constant; + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(rotorDiskSource, 0); + addToRunTimeSelectionTable(basicSource, rotorDiskSource, dictionary); + + template<> const char* NamedEnum<rotorDiskSource::geometryModeType, 2>:: + names[] = + { + "auto", + "specified" + }; + + const NamedEnum<rotorDiskSource::geometryModeType, 2> + rotorDiskSource::geometryModeTypeNames_; + + template<> const char* NamedEnum<rotorDiskSource::inletFlowType, 3>:: + names[] = + { + "fixed", + "surfaceNormal", + "local" + }; + + const NamedEnum<rotorDiskSource::inletFlowType, 3> + rotorDiskSource::inletFlowTypeNames_; +} + + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +void Foam::rotorDiskSource::checkData() +{ + switch (selectionMode()) + { + case smCellSet: + case smCellZone: + case smAll: + { + // set the profile ID for each blade section + profiles_.connectBlades(blade_.profileName(), blade_.profileID()); + switch (inletFlow_) + { + case ifFixed: + { + coeffs_.lookup("inletVelocity") >> inletVelocity_; + break; + } + case ifSurfaceNormal: + { + scalar UIn + ( + readScalar(coeffs_.lookup("inletNormalVelocity")) + ); + inletVelocity_ = -coordSys_.e3()*UIn; + break; + } + case ifLocal: + { + // do nothing + break; + } + default: + { + FatalErrorIn("void rotorDiskSource::checkData()") + << "Unknown inlet velocity type" << abort(FatalError); + } + } + + + break; + } + default: + { + FatalErrorIn("void rotorDiskSource::checkData()") + << "Source cannot be used with '" + << selectionModeTypeNames_[selectionMode()] + << "' mode. Please use one of: " << nl + << selectionModeTypeNames_[smCellSet] << nl + << selectionModeTypeNames_[smCellZone] << nl + << selectionModeTypeNames_[smAll] + << exit(FatalError); + } + } +} + + +void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct) +{ + // calculate rotor face areas + const vectorField& Sf = mesh_.Sf(); + const scalarField& magSf = mesh_.magSf(); + + boolList selectedCells(mesh_.nCells(), false); + boolList processedFaces(mesh_.nFaces(), false); + UIndirectList<bool>(selectedCells, cells_) = boolList(cells_.size(), true); + + vector n = vector::zero; + + label nFace = 0; + area_ = 0.0; + forAll(cells_, i) + { + const label cellI = cells_[i]; + const cell& cFaces = mesh_.cells()[cellI]; + forAll(cFaces, j) + { + const label faceI = cFaces[j]; + label own = mesh_.owner()[faceI]; + label nbr = mesh_.neighbour()[faceI]; + if + ( + !processedFaces[faceI] + && (selectedCells[own] != selectedCells[nbr]) + ) + { + if (selectedCells[own]) + { + if (((Sf[faceI]/magSf[faceI]) & axis) > 0.8) + { + area_[i] += magSf[faceI]; + n += Sf[faceI]; + nFace++; + } + } + else if (selectedCells[nbr]) + { + if (((-Sf[faceI]/magSf[faceI]) & axis) > 0.8) + { + area_[i] += magSf[faceI]; + n -= Sf[faceI]; + nFace++; + } + } + processedFaces[faceI] = true; + } + } + } + + if (correct && (nFace > 0)) + { + axis = n/mag(n); + } +} + + +void Foam::rotorDiskSource::createCoordinateSystem() +{ + // construct the local rotor co-prdinate system + vector origin(vector::zero); + vector axis(vector::zero); + vector refDir(vector::zero); + + geometryModeType gm = + geometryModeTypeNames_.read(coeffs_.lookup("geometryMode")); + + switch (gm) + { + case gmAuto: + { + // determine rotation origin + scalar sumV = 0.0; + const scalarField& V = mesh_.V(); + const vectorField& C = mesh_.C(); + forAll(cells_, i) + { + const label cellI = cells_[i]; + sumV += V[cellI]; + origin += V[cellI]*C[cellI]; + } + origin /= sumV; + + // determine first radial vector + vector dx1(vector::zero); + scalar magR = -GREAT; + forAll(cells_, i) + { + const label cellI = cells_[i]; + vector test = C[cellI] - origin; + if (mag(test) > magR) + { + dx1 = test; + magR = mag(test); + } + } + + // determine second radial vector and cross to determine axis + forAll(cells_, i) + { + const label cellI = cells_[i]; + vector dx2 = C[cellI] - origin; + if (mag(dx2) > 0.5*magR) + { + axis = dx1 ^ dx2; + if (mag(axis) > SMALL) + { + break; + } + } + } + axis /= mag(axis); + + // axis direction is somewhat arbitrary - check if user needs + // needs to reverse + bool reverse(readBool(coeffs_.lookup("reverseAxis"))); + if (reverse) + { + axis *= -1.0; + } + + coeffs_.lookup("refDirection") >> refDir; + + // set the face areas and apply correction to calculated axis + // e.g. if cellZone is more than a single layer in thickness + setFaceArea(axis, true); + + break; + } + case gmSpecified: + { + coeffs_.lookup("origin") >> origin; + coeffs_.lookup("axis") >> axis; + coeffs_.lookup("refDirection") >> refDir; + + setFaceArea(axis, false); + + break; + } + default: + { + FatalErrorIn + ( + "rotorDiskSource::createCoordinateSystem(const geometryMode&);" + ) << "Unknown geometryMode " << geometryModeTypeNames_[gm] + << ". Available geometry modes include " + << geometryModeTypeNames_ << exit(FatalError); + } + } + + coordSys_ = cylindricalCS("rotorCoordSys", origin, axis, refDir, false); + + const scalar sumArea = gSum(area_); + const scalar diameter = Foam::sqrt(4.0*sumArea/mathematical::pi); + Info<< " Rotor gometry:" << nl + << " - disk diameter = " << diameter << nl + << " - disk area = " << sumArea << nl + << " - origin = " << coordSys_.origin() << nl + << " - r-axis = " << coordSys_.e1() << nl + << " - psi-axis = " << coordSys_.e2() << nl + << " - z-axis = " << coordSys_.e3() << endl; +} + + +void Foam::rotorDiskSource::constructGeometry() +{ + const vectorField& C = mesh_.C(); + + const vector rDir = coordSys_.e1(); + const vector zDir = coordSys_.e3(); + + forAll(cells_, i) + { + const label cellI = cells_[i]; + + // position in rotor co-ordinate system + x_[i] = coordSys_.localPosition(C[cellI]); + + // cache max radius + rMax_ = max(rMax_, x_[i].x()); + + // determine swept angle relative to rDir axis + scalar psi = x_[i].y() - rDir.y(); + + // blade flap angle + scalar beta = flap_.beta0 - flap_.beta1*cos(psi) - flap_.beta2*sin(psi); + + // determine rotation tensor to convert into the rotor cone plane + scalar c = cos(-beta); + scalar s = sin(-beta); + R_[i] = tensor(1, 0, 0, 0, c, s, 0, -s, c); + + // geometric angle of attack - not including twist + alphag_[i] = trim_.alphaC - trim_.A*cos(psi) - trim_.B*sin(psi); + } +} + + +Foam::tmp<Foam::vectorField> Foam::rotorDiskSource::inflowVelocity +( + const volVectorField& U +) const +{ + switch (inletFlow_) + { + case ifFixed: + case ifSurfaceNormal: + { + return tmp<vectorField> + ( + new vectorField(mesh_.nCells(), inletVelocity_) + ); + + break; + } + case ifLocal: + { + return U.internalField(); + + break; + } + default: + { + FatalErrorIn + ( + "Foam::tmp<Foam::vectorField> " + "Foam::rotorDiskSource::inflowVelocity" + "(const volVectorField&) const" + ) << "Unknown inlet flow specification" << abort(FatalError); + } + } + + return tmp<vectorField>(new vectorField(mesh_.nCells(), vector::zero)); +} + + +// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * // + +Foam::rotorDiskSource::rotorDiskSource +( + const word& name, + const word& modelType, + const dictionary& dict, + const fvMesh& mesh + +) +: + basicSource(name, modelType, dict, mesh), + coeffs_(dict_.subDict(type() + "Coeffs")), + rhoName_("none"), + omega_(0.0), + nBlades_(0), + inletFlow_(ifLocal), + inletVelocity_(vector::zero), + tipEffect_(1.0), + flap_(), + trim_(), + blade_(coeffs_.subDict("blade")), + profiles_(coeffs_.subDict("profiles")), + x_(cells_.size(), vector::zero), + R_(cells_.size(), I), + alphag_(cells_.size(), 0.0), + area_(cells_.size(), 0.0), + coordSys_(false), + rMax_(0.0) +{ + read(dict); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::rotorDiskSource::~rotorDiskSource() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::rotorDiskSource::addSu(fvMatrix<vector>& UEqn) +{ + // add source to lhs of eqn + + const volVectorField& U = UEqn.psi(); + + if (UEqn.dimensions() == dimForce) + { + coeffs_.lookup("rhoName") >> rhoName_; + + const volScalarField& rho = + mesh_.lookupObject<volScalarField>(rhoName_); + + UEqn += calculateForces + ( + rho.internalField(), + inflowVelocity(U), + dimForce/dimVolume + ); + } + else + { + UEqn += calculateForces + ( + oneField(), + inflowVelocity(U), + dimForce/dimVolume/dimDensity + ); + } +} + + +void Foam::rotorDiskSource::addSu(fvMatrix<scalar>& UEqn) +{ + // do nothing +} + + +void Foam::rotorDiskSource::writeData(Ostream& os) const +{ + os << indent << name_ << endl; + dict_.write(os); +} + + +bool Foam::rotorDiskSource::read(const dictionary& dict) +{ + if (basicSource::read(dict)) + { + coeffs_ = dict.subDict(type() + "Coeffs"); + + scalar rpm(readScalar(coeffs_.lookup("rpm"))); + omega_ = rpm/60.0*mathematical::twoPi; + + coeffs_.lookup("nBlades") >> nBlades_; + + inletFlow_ = inletFlowTypeNames_.read(coeffs_.lookup("inletFlowType")); + + coeffs_.lookup("tipEffect") >> tipEffect_; + + const dictionary& flapCoeffs(coeffs_.subDict("flapCoeffs")); + flapCoeffs.lookup("beta0") >> flap_.beta0; + flapCoeffs.lookup("beta1") >> flap_.beta1; + flapCoeffs.lookup("beta2") >> flap_.beta2; + flap_.beta0 = degToRad(flap_.beta0); + + const dictionary& trimCoeffs(coeffs_.subDict("trimCoeffs")); + trimCoeffs.lookup("alphaC") >> trim_.alphaC; + trimCoeffs.lookup("A") >> trim_.A; + trimCoeffs.lookup("B") >> trim_.B; + trim_.alphaC = degToRad(trim_.alphaC); + + checkData(); + + createCoordinateSystem(); + + constructGeometry(); + + if (debug) + { + writeField("alphag", alphag_, true); + writeField("faceArea", area_, true); + } + + return true; + } + else + { + return false; + } +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H new file mode 100644 index 0000000000000000000000000000000000000000..0d646b98996b3174b4693da90e37346620adec46 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H @@ -0,0 +1,249 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Class + Foam::rotorDiskSource + +Description + Cell-zone based momemtum source + + Source approximates the mean effects of rotor forces on a cylindrical + region within the domain + +SourceFiles + rotorDiskSource.C + rotorDiskSourceTemplates.C + +\*---------------------------------------------------------------------------*/ + +#ifndef rotorDiskSource_H +#define rotorDiskSource_H + +#include "basicSource.H" +#include "cylindricalCS.H" +#include "NamedEnum.H" +#include "bladeModel.H" +#include "profileModelList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class rotorDiskSource Declaration +\*---------------------------------------------------------------------------*/ + +class rotorDiskSource +: + public basicSource +{ +public: + + enum geometryModeType + { + gmAuto, + gmSpecified + }; + static const NamedEnum<geometryModeType, 2> geometryModeTypeNames_; + + enum inletFlowType + { + ifFixed, + ifSurfaceNormal, + ifLocal + }; + static const NamedEnum<inletFlowType, 3> inletFlowTypeNames_; + + +protected: + + // Helper structures to encapsulate flap and trim data + + struct flapData + { + scalar beta0; // coning angle [deg] + scalar beta1; // lateral flapping coeff + scalar beta2; // longitudinal flapping coeff + }; + + struct trimData + { + scalar alphaC; // collective pitch angle [deg] + scalar A; // lateral cyclic coeff + scalar B; // longitudinal cyclic coeff + }; + + + // Protected data + + //- Coefficients dictionary + dictionary coeffs_; + + //- Name of density field + word rhoName_; + + //- Rotational speed [rad/s] + scalar omega_; + + //- Number of blades + label nBlades_; + + //- Inlet flow type + inletFlowType inletFlow_; + + //- Inlet velocity for specified iinflow + vector inletVelocity_; + + //- Tip effect [0-1] + // Ratio of blade radius beyond which lift=0 + scalar tipEffect_; + + //- Blade flap coefficients [rad/s] + flapData flap_; + + //- Blad trim coefficients + trimData trim_; + + //- Blade data + bladeModel blade_; + + //- Profile data + profileModelList profiles_; + + //- Cell centre positions in local rotor frame (Cartesian x, y, z) + List<point> x_; + + //- Rotation tensor for flap angle + List<tensor> R_; + + //- Geometric angle of attack [deg] + List<scalar> alphag_; + + //- Area [m2] + List<scalar> area_; + + //- Rotor co-ordinate system (r, theta, z) + cylindricalCS coordSys_; + + //- Maximum radius + scalar rMax_; + + + // Protected Member Functions + + //- Check data + void checkData(); + + //- Set the face areas per cell, and optionally correct the rotor axis + void setFaceArea(vector& axis, const bool correct); + + //- Create the co-ordinate system + void createCoordinateSystem(); + + //- Construct geometry + void constructGeometry(); + + //- Return the inlet flow field + tmp<vectorField> inflowVelocity(const volVectorField& U) const; + + //- Calculate forces + template<class RhoType> + tmp<volVectorField> calculateForces + ( + const RhoType& rho, + const vectorField& U, + const dimensionSet& dims + ); + + //- Helper function to write rotor values + template<class Type> + void writeField + ( + const word& name, + const List<Type>& values, + const bool writeNow = false + ) const; + + +public: + + //- Runtime type information + TypeName("rotorDisk"); + + + // Constructors + + + //- Construct from components + rotorDiskSource + ( + const word& name, + const word& modelType, + const dictionary& dict, + const fvMesh& mesh + ); + + + //- Destructor + virtual ~rotorDiskSource(); + + + // Member Functions + + + // Source term addition + + //- Source term to fvMatrix<vector> + virtual void addSu(fvMatrix<vector>& UEqn); + + //- Source term to fvMatrix<scalar> + virtual void addSu(fvMatrix<scalar>& UEqn); + + + // I-O + + //- Write the source properties + virtual void writeData(Ostream&) const; + + //- Read source dictionary + virtual bool read(const dictionary& dict); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "rotorDiskSourceTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // + diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceTemplates.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..5bc97c703f99f4ab4e4b316b0c1c9213a31c2fc5 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceTemplates.C @@ -0,0 +1,205 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "rotorDiskSource.H" +#include "addToRunTimeSelectionTable.H" +#include "mathematicalConstants.H" +#include "unitConversion.H" + +using namespace Foam::constant; + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +template<class Type> +void Foam::rotorDiskSource::writeField +( + const word& name, + const List<Type>& values, + const bool writeNow +) const +{ + typedef GeometricField<Type, fvPatchField, volMesh> fieldType; + + if (mesh_.time().outputTime() || writeNow) + { + tmp<fieldType> tfld + ( + new fieldType + ( + IOobject + ( + name, + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensioned<Type>("zero", dimless, pTraits<Type>::zero) + ) + ); + + Field<Type>& fld = tfld().internalField(); + + if (cells_.size() != values.size()) + { + FatalErrorIn("") << "cells_.size() != values_.size()" + << abort(FatalError); + } + + forAll(cells_, i) + { + const label cellI = cells_[i]; + fld[cellI] = values[i]; + } + + tfld().write(); + } +} + + +template<class RhoType> +Foam::tmp<Foam::volVectorField> Foam::rotorDiskSource::calculateForces +( + const RhoType& rho, + const vectorField& U, + const dimensionSet& dims +) +{ + tmp<volVectorField> tForce + ( + new volVectorField + ( + IOobject + ( + "rotorForce", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedVector("zero", dims, vector::zero) + ) + ); + + vectorField& force = tForce().internalField(); + const scalarField& V = mesh_.V(); + + + // logging info + scalar dragEff = 0.0; + scalar liftEff = 0.0; + scalar AOAmin = GREAT; + scalar AOAmax = -GREAT; + + forAll(cells_, i) + { + if (area_[i] > ROOTVSMALL) + { + const label cellI = cells_[i]; + + const scalar radius = x_[i].x(); + + // apply correction due to flap in cartesian frame + vector Uc = R_[i] & U[cellI]; + + // velocity in local reference frame + Uc = coordSys_.localVector(Uc); + + // set radial component of velocity to zero + Uc.x() = 0.0; + + // remove blade linear velocity from blade normal component + Uc.y() -= radius*omega_; + + // velocity magnitude + scalar magUc = mag(Uc); + + // determine blade data for this radius + // i1 = index of upper bound data point in blade list + scalar twist = 0.0; + scalar chord = 0.0; + label i1 = -1; + label i2 = -1; + scalar invDr = 0.0; + blade_.interpolate(radius, twist, chord, i1, i2, invDr); + + // effective angle of attack + scalar alphaEff = alphag_[i] + twist - atan(Uc.z()/Uc.y()); + AOAmin = min(AOAmin, alphaEff); + AOAmax = max(AOAmax, alphaEff); + + // determine profile data for this radius and angle of attack + const label profile1 = blade_.profileID()[i1]; + const label profile2 = blade_.profileID()[i2]; + + scalar Cd1 = 0.0; + scalar Cl1 = 0.0; + profiles_[profile1].Cdl(alphaEff, Cd1, Cl1); + + scalar Cd2 = 0.0; + scalar Cl2 = 0.0; + profiles_[profile2].Cdl(alphaEff, Cd2, Cl2); + + scalar Cd = invDr*(Cd2 - Cd1) + Cd1; + scalar Cl = invDr*(Cl2 - Cl1) + Cl1; + + // apply tip effect for blade lift + scalar tipFactor = 1.0; + if (radius/rMax_ > tipEffect_) + { + tipFactor = 0.0; + } + + // calculate forces + scalar pDyn = 0.5*rho[cellI]*sqr(magUc); + scalar f = pDyn*chord*nBlades_*area_[i]/(mathematical::twoPi); + vector localForce = vector(0.0, f*Cd, tipFactor*f*Cl); + + // accumulate forces + dragEff += localForce.y(); + liftEff += localForce.z(); + + // convert force to global cartesian co-ordinate system + force[cellI] = coordSys_.globalVector(localForce); + + force[cellI] /= V[cellI]; + } + } + + + Info<< type() << " output:" << nl + << " min/max(AOA) = " << radToDeg(AOAmin) << ", " + << radToDeg(AOAmax) << nl + << " Effective drag = " << dragEff << nl + << " Effective lift = " << liftEff << endl; + + + return tForce; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C index 09587c986d27a9e220bf18eae70367f9f0b142fa..d3d0116805e5d091d07f9a42d6a2819cce4a1df9 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C @@ -134,7 +134,7 @@ Foam::fv::gradScheme<Type>::grad regIOobject::store(tgGrad.ptr()); } - cachePrintMessage("Retreiving", name, vsf); + cachePrintMessage("Retrieving", name, vsf); GradFieldType& gGrad = const_cast<GradFieldType&> ( mesh().objectRegistry::template lookupObject<GradFieldType>(name) diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C index d080a624a2f68ca809fbcac43d5bdeaa1716ab1b..af83ca14295550c13f6540812abd0dd43daf5b78 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C @@ -239,7 +239,7 @@ void Foam::SurfaceFilmModel<CloudType>::cacheFilmFields diameterParcelPatch_ = filmModel.cloudDiameterTrans().boundaryField()[filmPatchI]; - filmModel.toPrimary(filmPatchI, diameterParcelPatch_); + filmModel.toPrimary(filmPatchI, diameterParcelPatch_, maxOp<scalar>()); UFilmPatch_ = filmModel.Us().boundaryField()[filmPatchI]; filmModel.toPrimary(filmPatchI, UFilmPatch_); diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C index 535583a9e13578de3a750ce3804e15122273fc17..2e6a4b4384da5fab4fd8602aa60c85a1295dd619 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C @@ -1333,11 +1333,12 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update template<class SourcePatch, class TargetPatch> -template<class Type> +template<class Type, class BinaryOp> Foam::tmp<Foam::Field<Type> > Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource ( - const Field<Type>& fld + const Field<Type>& fld, + const BinaryOp& bop ) const { if (fld.size() != tgtAddress_.size()) @@ -1377,7 +1378,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource forAll(faces, i) { - result[faceI] += work[faces[i]]*weights[i]; + result[faceI] = bop(result[faceI], work[faces[i]]*weights[i]); } } } @@ -1390,7 +1391,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource forAll(faces, i) { - result[faceI] += fld[faces[i]]*weights[i]; + result[faceI] = bop(result[faceI], fld[faces[i]]*weights[i]); } } } @@ -1400,24 +1401,26 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource template<class SourcePatch, class TargetPatch> -template<class Type> +template<class Type, class BinaryOp> Foam::tmp<Foam::Field<Type> > Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource ( - const tmp<Field<Type> >& tFld + const tmp<Field<Type> >& tFld, + const BinaryOp& bop ) const { - return interpolateToSource(tFld()); + return interpolateToSource(tFld(), bop); } template<class SourcePatch, class TargetPatch> -template<class Type> +template<class Type, class BinaryOp> Foam::tmp<Foam::Field<Type> > Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget ( - const Field<Type>& fld + const Field<Type>& fld, + const BinaryOp& bop ) const { if (fld.size() != srcAddress_.size()) @@ -1457,7 +1460,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget forAll(faces, i) { - result[faceI] += work[faces[i]]*weights[i]; + result[faceI] = bop(result[faceI], work[faces[i]]*weights[i]); } } } @@ -1470,7 +1473,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget forAll(faces, i) { - result[faceI] += fld[faces[i]]*weights[i]; + result[faceI] = bop(result[faceI], fld[faces[i]]*weights[i]); } } } @@ -1479,6 +1482,55 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget } +template<class SourcePatch, class TargetPatch> +template<class Type, class BinaryOp> +Foam::tmp<Foam::Field<Type> > +Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget +( + const tmp<Field<Type> >& tFld, + const BinaryOp& bop +) const +{ + return interpolateToTarget(tFld(), bop); +} + + +template<class SourcePatch, class TargetPatch> +template<class Type> +Foam::tmp<Foam::Field<Type> > +Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource +( + const Field<Type>& fld +) const +{ + return interpolateToSource(fld, sumOp<Type>()); +} + + +template<class SourcePatch, class TargetPatch> +template<class Type> +Foam::tmp<Foam::Field<Type> > +Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource +( + const tmp<Field<Type> >& tFld +) const +{ + return interpolateToSource(tFld, sumOp<Type>()); +} + + +template<class SourcePatch, class TargetPatch> +template<class Type> +Foam::tmp<Foam::Field<Type> > +Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget +( + const Field<Type>& fld +) const +{ + return interpolateToSource(fld, sumOp<Type>()); +} + + template<class SourcePatch, class TargetPatch> template<class Type> Foam::tmp<Foam::Field<Type> > @@ -1487,7 +1539,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget const tmp<Field<Type> >& tFld ) const { - return interpolateToTarget(tFld()); + return interpolateToSource(tFld, sumOp<Type>()); } diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H index b9d73f24ee3446bae541fea490636b7d0a56df8d..8df0d6fb4abaa105efca9f31ce858ecef62a1bcc 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H @@ -53,6 +53,7 @@ SourceFiles #include "treeBoundBox.H" #include "treeBoundBoxList.H" #include "globalIndex.H" +#include "ops.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -345,11 +346,46 @@ public: // Evaluation + //- Interpolate from target to source with supplied op + template<class Type, class BinaryOp> + tmp<Field<Type> > interpolateToSource + ( + const Field<Type>& fld, + const BinaryOp& bop + ) const; + + //- Interpolate from target tmp field to source with supplied op + template<class Type, class BinaryOp> + tmp<Field<Type> > interpolateToSource + ( + const tmp<Field<Type> >& tFld, + const BinaryOp& bop + ) const; + + //- Interpolate from source to target with supplied op + template<class Type, class BinaryOp> + tmp<Field<Type> > interpolateToTarget + ( + const Field<Type>& fld, + const BinaryOp& bop + ) const; + + //- Interpolate from source tmp field to target with supplied op + template<class Type, class BinaryOp> + tmp<Field<Type> > interpolateToTarget + ( + const tmp<Field<Type> >& tFld, + const BinaryOp& bop + ) const; + //- Interpolate from target to source template<class Type> - tmp<Field<Type> > interpolateToSource(const Field<Type>& fld) const; + tmp<Field<Type> > interpolateToSource + ( + const Field<Type>& fld + ) const; - //- Interpolate from target tmp field to source + //- Interpolate from target tmp field template<class Type> tmp<Field<Type> > interpolateToSource ( @@ -358,9 +394,12 @@ public: //- Interpolate from source to target template<class Type> - tmp<Field<Type> > interpolateToTarget(const Field<Type>& fld) const; + tmp<Field<Type> > interpolateToTarget + ( + const Field<Type>& fld + ) const; - //- Interpolate from source tmp field to target + //- Interpolate from source tmp field template<class Type> tmp<Field<Type> > interpolateToTarget ( diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index 51be5477869f8e60a1c123cde24fa37eaa8a5aa0..1c981f478877461e4a1af752630dbb48c6f9eb1b 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -30,8 +30,13 @@ meshSearch/meshSearch.C meshTools/meshTools.C -PointEdgeWave/PointEdgeWaveName.C -PointEdgeWave/pointEdgePoint.C +pWave = PointEdgeWave +$(pWave)/PointEdgeWaveName.C +$(pWave)/pointEdgePoint.C + +patchWave = PatchEdgeFaceWave +$(patchWave)/PatchEdgeFaceWaveName.C +$(patchWave)/patchEdgeFaceInfo.C regionSplit/regionSplit.C diff --git a/src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWave.C b/src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWave.C new file mode 100644 index 0000000000000000000000000000000000000000..0ac91b03c83f1b5b85931d08ecb20bfabc8cb527 --- /dev/null +++ b/src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWave.C @@ -0,0 +1,681 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "PatchEdgeFaceWave.H" +#include "polyMesh.H" +#include "globalMeshData.H" +#include "PatchTools.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +Foam::scalar Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +propagationTol_ = 0.01; + +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +Foam::label +Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +dummyTrackData_ = 12345; + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +// Update info for edgeI, at position pt, with information from +// neighbouring face. +// Updates: +// - changedEdge_, changedEdges_, +// - statistics: nEvals_, nUnvisitedEdges_ +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +bool Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +updateEdge +( + const label edgeI, + const label neighbourFaceI, + const Type& neighbourInfo, + Type& edgeInfo +) +{ + nEvals_++; + + bool wasValid = edgeInfo.valid(td_); + + bool propagate = + edgeInfo.updateEdge + ( + mesh_, + patch_, + edgeI, + neighbourFaceI, + neighbourInfo, + propagationTol_, + td_ + ); + + if (propagate) + { + if (!changedEdge_[edgeI]) + { + changedEdge_[edgeI] = true; + changedEdges_.append(edgeI); + } + } + + if (!wasValid && edgeInfo.valid(td_)) + { + --nUnvisitedEdges_; + } + + return propagate; +} + + +// Update info for faceI, at position pt, with information from +// neighbouring edge. +// Updates: +// - changedFace_, changedFaces_, +// - statistics: nEvals_, nUnvisitedFace_ +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +bool Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +updateFace +( + const label faceI, + const label neighbourEdgeI, + const Type& neighbourInfo, + Type& faceInfo +) +{ + nEvals_++; + + bool wasValid = faceInfo.valid(td_); + + bool propagate = + faceInfo.updateFace + ( + mesh_, + patch_, + faceI, + neighbourEdgeI, + neighbourInfo, + propagationTol_, + td_ + ); + + if (propagate) + { + if (!changedFace_[faceI]) + { + changedFace_[faceI] = true; + changedFaces_.append(faceI); + } + } + + if (!wasValid && faceInfo.valid(td_)) + { + --nUnvisitedFaces_; + } + + return propagate; +} + + +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +void Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +syncEdges() +{ + const globalMeshData& globalData = mesh_.globalData(); + const mapDistribute& map = globalData.globalEdgeSlavesMap(); + const PackedBoolList& cppOrientation = globalData.globalEdgeOrientation(); + + // Convert patch-edge data into cpp-edge data + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + //- Construct with all data in consistent orientation + List<Type> cppEdgeData(map.constructSize()); + + forAll(patchEdges_, i) + { + label patchEdgeI = patchEdges_[i]; + label coupledEdgeI = coupledEdges_[i]; + + if (changedEdge_[patchEdgeI]) + { + const Type& data = allEdgeInfo_[patchEdgeI]; + + // Patch-edge data needs to be converted into coupled-edge data + // (optionally flipped) and consistent in orientation with + // master of coupled edge (optionally flipped) + bool sameOrientation = + ( + sameEdgeOrientation_[i] + == cppOrientation[coupledEdgeI] + ); + + cppEdgeData[coupledEdgeI].updateEdge + ( + mesh_, + patch_, + data, + sameOrientation, + propagationTol_, + td_ + ); + } + } + + + // Synchronise + // ~~~~~~~~~~~ + + globalData.syncData + ( + cppEdgeData, + globalData.globalEdgeSlaves(), + globalData.globalEdgeTransformedSlaves(), + map, + globalData.globalTransforms(), + updateOp<PrimitivePatchType, Type, TrackingData> + ( + mesh_, + patch_, + propagationTol_, + td_ + ), + transformOp<PrimitivePatchType, Type, TrackingData> + ( + mesh_, + patch_, + propagationTol_, + td_ + ) + ); + + + // Back from cpp-edge to patch-edge data + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + forAll(patchEdges_, i) + { + label patchEdgeI = patchEdges_[i]; + label coupledEdgeI = coupledEdges_[i]; + + const Type& data = cppEdgeData[coupledEdgeI]; + + if (data.valid(td_)) + { + bool sameOrientation = + ( + sameEdgeOrientation_[i] + == cppOrientation[coupledEdgeI] + ); + + allEdgeInfo_[patchEdgeI].updateEdge + ( + mesh_, + patch_, + data, + sameOrientation, + propagationTol_, + td_ + ); + + if (!changedEdge_[patchEdgeI]) + { + changedEdges_.append(patchEdgeI); + changedEdge_[patchEdgeI] = true; + } + } + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Iterate, propagating changedEdgesInfo across patch, until no change (or +// maxIter reached). Initial edge values specified. +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +PatchEdgeFaceWave +( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const labelList& changedEdges, + const List<Type>& changedEdgesInfo, + + UList<Type>& allEdgeInfo, + UList<Type>& allFaceInfo, + + const label maxIter, + TrackingData& td +) +: + mesh_(mesh), + patch_(patch), + allEdgeInfo_(allEdgeInfo), + allFaceInfo_(allFaceInfo), + td_(td), + changedEdge_(patch_.nEdges()), + changedEdges_(patch_.size()), + changedFace_(patch_.size()), + changedFaces_(patch_.size()), + nEvals_(0), + nUnvisitedEdges_(patch_.nEdges()), + nUnvisitedFaces_(patch_.size()) +{ + // Calculate addressing between patch_ and mesh.globalData().coupledPatch() + // for ease of synchronisation + PatchTools::matchEdges + ( + patch_, + mesh_.globalData().coupledPatch(), + + patchEdges_, + coupledEdges_, + sameEdgeOrientation_ + ); + + + if (allEdgeInfo_.size() != patch_.nEdges()) + { + FatalErrorIn + ( + "PatchEdgeFaceWave<Type, TrackingData>::PatchEdgeFaceWave" + "(const polyMesh&, const labelList&, const List<Type>," + " List<Type>&, List<Type>&, const label maxIter)" + ) << "size of edgeInfo work array is not equal to the number" + << " of edges in the patch" << endl + << " edgeInfo :" << allEdgeInfo_.size() << endl + << " patch.nEdges:" << patch_.nEdges() + << exit(FatalError); + } + if (allFaceInfo_.size() != patch_.size()) + { + FatalErrorIn + ( + "PatchEdgeFaceWave<Type, TrackingData>::PatchEdgeFaceWave" + "(const polyMesh&, const labelList&, const List<Type>," + " List<Type>&, List<Type>&, const label maxIter)" + ) << "size of edgeInfo work array is not equal to the number" + << " of faces in the patch" << endl + << " faceInfo :" << allFaceInfo_.size() << endl + << " patch.size:" << patch_.size() + << exit(FatalError); + } + + + // Set from initial changed edges data + setEdgeInfo(changedEdges, changedEdgesInfo); + + if (debug) + { + Pout<< "Seed edges : " << changedEdges_.size() << endl; + } + + // Iterate until nothing changes + label iter = iterate(maxIter); + + if ((maxIter > 0) && (iter >= maxIter)) + { + FatalErrorIn + ( + "PatchEdgeFaceWave<Type, TrackingData>::PatchEdgeFaceWave" + "(const polyMesh&, const labelList&, const List<Type>," + " List<Type>&, List<Type>&, const label maxIter)" + ) << "Maximum number of iterations reached. Increase maxIter." << endl + << " maxIter:" << maxIter << endl + << " changedEdges:" << changedEdges_.size() << endl + << " changedFaces:" << changedFaces_.size() << endl + << exit(FatalError); + } +} + + +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +PatchEdgeFaceWave +( + const polyMesh& mesh, + const PrimitivePatchType& patch, + UList<Type>& allEdgeInfo, + UList<Type>& allFaceInfo, + TrackingData& td +) +: + mesh_(mesh), + patch_(patch), + allEdgeInfo_(allEdgeInfo), + allFaceInfo_(allFaceInfo), + td_(td), + changedEdge_(patch_.nEdges()), + changedEdges_(patch_.nEdges()), + changedFace_(patch_.size()), + changedFaces_(patch_.size()), + nEvals_(0), + nUnvisitedEdges_(patch_.nEdges()), + nUnvisitedFaces_(patch_.size()) +{ + // Calculate addressing between patch_ and mesh.globalData().coupledPatch() + // for ease of synchronisation + PatchTools::matchEdges + ( + patch_, + mesh_.globalData().coupledPatch(), + + patchEdges_, + coupledEdges_, + sameEdgeOrientation_ + ); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +getUnsetEdges() const +{ + return nUnvisitedEdges_; +} + + +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +getUnsetFaces() const +{ + return nUnvisitedFaces_; +} + + +// Copy edge information into member data +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +void Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +setEdgeInfo +( + const labelList& changedEdges, + const List<Type>& changedEdgesInfo +) +{ + forAll(changedEdges, changedEdgeI) + { + label edgeI = changedEdges[changedEdgeI]; + + bool wasValid = allEdgeInfo_[edgeI].valid(td_); + + // Copy info for edgeI + allEdgeInfo_[edgeI] = changedEdgesInfo[changedEdgeI]; + + // Maintain count of unset edges + if (!wasValid && allEdgeInfo_[edgeI].valid(td_)) + { + --nUnvisitedEdges_; + } + + // Mark edgeI as changed, both on list and on edge itself. + + if (!changedEdge_[edgeI]) + { + changedEdge_[edgeI] = true; + changedEdges_.append(edgeI); + } + } +} + + +// Propagate information from face to edge. Return number of edges changed. +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +faceToEdge() +{ + changedEdges_.clear(); + changedEdge_ = false; + + forAll(changedFaces_, changedFaceI) + { + label faceI = changedFaces_[changedFaceI]; + + if (!changedFace_[faceI]) + { + FatalErrorIn("PatchEdgeFaceWave<Type, TrackingData>::faceToEdge()") + << "face " << faceI + << " not marked as having been changed" << nl + << "This might be caused by multiple occurences of the same" + << " seed edge." << abort(FatalError); + } + + const Type& neighbourWallInfo = allFaceInfo_[faceI]; + + // Evaluate all connected edges + const labelList& fEdges = patch_.faceEdges()[faceI]; + + forAll(fEdges, fEdgeI) + { + label edgeI = fEdges[fEdgeI]; + + Type& currentWallInfo = allEdgeInfo_[edgeI]; + + if (!currentWallInfo.equal(neighbourWallInfo, td_)) + { + updateEdge + ( + edgeI, + faceI, + neighbourWallInfo, + currentWallInfo + ); + } + } + } + + + syncEdges(); + + + if (debug) + { + Pout<< "Changed edges : " << changedEdges_.size() << endl; + } + + return returnReduce(changedEdges_.size(), sumOp<label>()); +} + + +// Propagate information from edge to face. Return number of faces changed. +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +edgeToFace() +{ + changedFaces_.clear(); + changedFace_ = false; + + const labelListList& edgeFaces = patch_.edgeFaces(); + + forAll(changedEdges_, changedEdgeI) + { + label edgeI = changedEdges_[changedEdgeI]; + + if (!changedEdge_[edgeI]) + { + FatalErrorIn("PatchEdgeFaceWave<Type, TrackingData>::edgeToFace()") + << "edge " << edgeI + << " not marked as having been changed" << nl + << "This might be caused by multiple occurences of the same" + << " seed edge." << abort(FatalError); + } + + const Type& neighbourWallInfo = allEdgeInfo_[edgeI]; + + // Evaluate all connected faces + + const labelList& eFaces = edgeFaces[edgeI]; + forAll(eFaces, eFaceI) + { + label faceI = eFaces[eFaceI]; + + Type& currentWallInfo = allFaceInfo_[faceI]; + + if (!currentWallInfo.equal(neighbourWallInfo, td_)) + { + updateFace + ( + faceI, + edgeI, + neighbourWallInfo, + currentWallInfo + ); + } + } + } + + if (debug) + { + Pout<< "Changed faces : " << changedFaces_.size() << endl; + } + + return returnReduce(changedFaces_.size(), sumOp<label>()); +} + + +// Iterate +template +< + class PrimitivePatchType, + class Type, + class TrackingData +> +Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>:: +iterate +( + const label maxIter +) +{ + // Make sure coupled edges contain same info + syncEdges(); + + nEvals_ = 0; + + label iter = 0; + + while (iter < maxIter) + { + if (debug) + { + Pout<< "Iteration " << iter << endl; + } + + label nFaces = edgeToFace(); + + if (debug) + { + Pout<< "Total changed faces : " << nFaces << endl; + } + + if (nFaces == 0) + { + break; + } + + label nEdges = faceToEdge(); + + if (debug) + { + Pout<< "Total changed edges : " << nEdges << nl + << "Total evaluations : " << nEvals_ << nl + << "Remaining unvisited edges : " << nUnvisitedEdges_ << nl + << "Remaining unvisited faces : " << nUnvisitedFaces_ << nl + << endl; + } + + if (nEdges == 0) + { + break; + } + + iter++; + } + + return iter; +} + + +// ************************************************************************* // diff --git a/src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWave.H b/src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWave.H new file mode 100644 index 0000000000000000000000000000000000000000..ff99788f7e6bc0d39321bf462ada2e2795e27f96 --- /dev/null +++ b/src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWave.H @@ -0,0 +1,368 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Class + Foam::PatchEdgeFaceWave + +Description + Wave propagation of information along patch. Every iteration + information goes through one layer of faces. Templated on information + that is transferred. + +SourceFiles + PatchEdgeFaceWave.C + +\*---------------------------------------------------------------------------*/ + +#ifndef PatchEdgeFaceWave_H +#define PatchEdgeFaceWave_H + +#include "scalarField.H" +#include "PackedBoolList.H" +#include "PrimitivePatch.H" +#include "vectorTensorTransform.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class polyMesh; + +/*---------------------------------------------------------------------------*\ + Class PatchEdgeFaceWaveName Declaration +\*---------------------------------------------------------------------------*/ + +TemplateName(PatchEdgeFaceWave); + + +/*---------------------------------------------------------------------------*\ + Class PatchEdgeFaceWave Declaration +\*---------------------------------------------------------------------------*/ + +template +< + class PrimitivePatchType, + class Type, + class TrackingData = int +> +class PatchEdgeFaceWave +: + public PatchEdgeFaceWaveName +{ + // Private static data + + //- Relative tolerance. Stop propagation if relative changes + // less than this tolerance (responsability for checking this is + // up to Type implementation) + static scalar propagationTol_; + + //- Used as default trackdata value to satisfy default template + // argument. + static label dummyTrackData_; + + + // Private data + + //- Reference to mesh + const polyMesh& mesh_; + + //- Reference to patch + const PrimitivePatchType& patch_; + + //- Wall information for all edges + UList<Type>& allEdgeInfo_; + + //- Information on all patch faces + UList<Type>& allFaceInfo_; + + //- Additional data to be passed into container + TrackingData& td_; + + //- Has edge changed + PackedBoolList changedEdge_; + + //- List of changed edges + DynamicList<label> changedEdges_; + + //- Has face changed + PackedBoolList changedFace_; + + //- List of changed faces + DynamicList<label> changedFaces_; + + //- Number of evaluations + label nEvals_; + + //- Number of unvisited faces/edges + label nUnvisitedEdges_; + label nUnvisitedFaces_; + + + // Addressing between edges of patch_ and globalData.coupledPatch() + labelList patchEdges_; + labelList coupledEdges_; + PackedBoolList sameEdgeOrientation_; + + + // Private Member Functions + + //- Updates edgeInfo with information from neighbour. Updates all + // statistics. + bool updateEdge + ( + const label edgeI, + const label neighbourFaceI, + const Type& neighbourInfo, + Type& edgeInfo + ); + + //- Updates faceInfo with information from neighbour. Updates all + // statistics. + bool updateFace + ( + const label faceI, + const label neighbourEdgeI, + const Type& neighbourInfo, + Type& faceInfo + ); + + //- Update coupled edges + void syncEdges(); + + //- Disallow default bitwise copy construct + PatchEdgeFaceWave(const PatchEdgeFaceWave&); + + //- Disallow default bitwise assignment + void operator=(const PatchEdgeFaceWave&); + + +public: + + // Static Functions + + //- Access to tolerance + static scalar propagationTol() + { + return propagationTol_; + } + + //- Change tolerance + static void setPropagationTol(const scalar tol) + { + propagationTol_ = tol; + } + + + // Constructors + + //- Construct from patch, list of changed edges with the Type + // for these edges. Gets work arrays to operate on, one of size + // number of patch edges, the other number of patch faces. + // Iterates until nothing changes or maxIter reached. + // (maxIter can be 0) + PatchEdgeFaceWave + ( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const labelList& initialEdges, + const List<Type>& initialEdgesInfo, + UList<Type>& allEdgeInfo, + UList<Type>& allFaceInfo, + const label maxIter, + TrackingData& td = dummyTrackData_ + ); + + //- Construct from patch. Use setEdgeInfo and iterate() to do + // actual calculation + PatchEdgeFaceWave + ( + const polyMesh& mesh, + const PrimitivePatchType& patch, + UList<Type>& allEdgeInfo, + UList<Type>& allFaceInfo, + TrackingData& td = dummyTrackData_ + ); + + + // Member Functions + + //- Access allEdgeInfo + UList<Type>& allEdgeInfo() const + { + return allEdgeInfo_; + } + + //- Access allFaceInfo + UList<Type>& allFaceInfo() const + { + return allFaceInfo_; + } + + //- Additional data to be passed into container + const TrackingData& data() const + { + return td_; + } + + //- Get number of unvisited faces, i.e. faces that were not (yet) + // reached from walking across patch. This can happen from + // - not enough iterations done + // - a disconnected patch + // - a patch without walls in it + label getUnsetFaces() const; + + label getUnsetEdges() const; + + //- Copy initial data into allEdgeInfo_ + void setEdgeInfo + ( + const labelList& changedEdges, + const List<Type>& changedEdgesInfo + ); + + //- Propagate from edge to face. Returns total number of faces + // (over all processors) changed. + label edgeToFace(); + + //- Propagate from face to edge. Returns total number of edges + // (over all processors) changed. + label faceToEdge(); + + //- Iterate until no changes or maxIter reached. Returns actual + // number of iterations. + label iterate(const label maxIter); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//- Update operation +template +< + class PrimitivePatchType, + class Type, + class TrackingData = int +> +class updateOp +{ + //- Additional data to be passed into container + const polyMesh& mesh_; + const PrimitivePatchType& patch_; + const scalar tol_; + TrackingData& td_; + +public: + updateOp + ( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const scalar tol, + TrackingData& td + ) + : + mesh_(mesh), + patch_(patch), + tol_(tol), + td_(td) + {} + + void operator()(Type& x, const Type& y) const + { + if (y.valid(td_)) + { + x.updateEdge(mesh_, patch_, y, true, tol_, td_); + } + } +}; + + +//- Transform operation +template +< + class PrimitivePatchType, + class Type, + class TrackingData = int +> +class transformOp +{ + //- Additional data to be passed into container + const polyMesh& mesh_; + const PrimitivePatchType& patch_; + const scalar tol_; + TrackingData& td_; + +public: + transformOp + ( + const polyMesh& mesh, + const PrimitivePatchType& patch, + const scalar tol, + TrackingData& td + ) + : + mesh_(mesh), + patch_(patch), + tol_(tol), + td_(td) + {} + + void operator() + ( + const vectorTensorTransform& vt, + const bool forward, + List<Type>& fld + ) const + { + if (forward) + { + forAll(fld, i) + { + fld[i].transform(mesh_, patch_, vt.R(), tol_, td_); + } + } + else + { + forAll(fld, i) + { + fld[i].transform(mesh_, patch_, vt.R().T(), tol_, td_); + } + } + } +}; + +} // End namespace Foam + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "PatchEdgeFaceWave.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWaveName.C b/src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWaveName.C new file mode 100644 index 0000000000000000000000000000000000000000..5d3f3e1d581a7219bb0f1b9ce2f8f39b86753fc5 --- /dev/null +++ b/src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWaveName.C @@ -0,0 +1,32 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "PatchEdgeFaceWave.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(Foam::PatchEdgeFaceWaveName, 0); + +// ************************************************************************* // diff --git a/src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfo.C b/src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfo.C new file mode 100644 index 0000000000000000000000000000000000000000..f6e72fed99820e7cfd1c6a3f88d40e0611913592 --- /dev/null +++ b/src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfo.C @@ -0,0 +1,50 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "patchEdgeFaceInfo.H" + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + +Foam::Ostream& Foam::operator<< +( + Foam::Ostream& os, + const Foam::patchEdgeFaceInfo& wDist +) +{ + return os << wDist.origin() << wDist.distSqr(); +} + + +Foam::Istream& Foam::operator>> +( + Foam::Istream& is, + Foam::patchEdgeFaceInfo& wDist +) +{ + return is >> wDist.origin_ >> wDist.distSqr_; +} + + +// ************************************************************************* // diff --git a/src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfo.H b/src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfo.H new file mode 100644 index 0000000000000000000000000000000000000000..f02120d7d9a5655da97acfbf706956d9f44e8afd --- /dev/null +++ b/src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfo.H @@ -0,0 +1,212 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Class + Foam::patchEdgeFaceInfo + +Description + +SourceFiles + patchEdgeFaceInfoI.H + patchEdgeFaceInfo.C + +\*---------------------------------------------------------------------------*/ + +#ifndef patchEdgeFaceInfo_H +#define patchEdgeFaceInfo_H + +#include "point.H" +#include "label.H" +#include "scalar.H" +#include "tensor.H" +#include "pTraits.H" +#include "primitivePatch.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class polyPatch; +class polyMesh; + +/*---------------------------------------------------------------------------*\ + Class patchEdgeFaceInfo Declaration +\*---------------------------------------------------------------------------*/ + +class patchEdgeFaceInfo +{ + // Private data + + //- position of nearest wall center + point origin_; + + //- normal distance (squared) from point to origin + scalar distSqr_; + + + // Private Member Functions + + //- Evaluate distance to point. Update distSqr, origin from whomever + // is nearer pt. Return true if w2 is closer to point, + // false otherwise. + template<class TrackingData> + inline bool update + ( + const point&, + const patchEdgeFaceInfo& w2, + const scalar tol, + TrackingData& td + ); + + //- Combine current with w2. Update distSqr, origin if w2 has smaller + // quantities and returns true. + template<class TrackingData> + inline bool update + ( + const patchEdgeFaceInfo& w2, + const scalar tol, + TrackingData& td + ); + + +public: + + // Constructors + + //- Construct null + inline patchEdgeFaceInfo(); + + //- Construct from origin, distance + inline patchEdgeFaceInfo(const point&, const scalar); + + //- Construct as copy + inline patchEdgeFaceInfo(const patchEdgeFaceInfo&); + + + // Member Functions + + // Access + + inline const point& origin() const; + + inline scalar distSqr() const; + + + // Needed by meshWave + + //- Check whether origin has been changed at all or + // still contains original (invalid) value. + template<class TrackingData> + inline bool valid(TrackingData& td) const; + + //- Apply rotation matrix + template<class TrackingData> + inline void transform + ( + const polyMesh& mesh, + const primitivePatch& patch, + const tensor& rotTensor, + const scalar tol, + TrackingData& td + ); + + //- Influence of face on edge + template<class TrackingData> + inline bool updateEdge + ( + const polyMesh& mesh, + const primitivePatch& patch, + const label edgeI, + const label faceI, + const patchEdgeFaceInfo& faceInfo, + const scalar tol, + TrackingData& td + ); + + //- New information for edge (from e.g. coupled edge) + template<class TrackingData> + inline bool updateEdge + ( + const polyMesh& mesh, + const primitivePatch& patch, + const patchEdgeFaceInfo& edgeInfo, + const bool sameOrientation, + const scalar tol, + TrackingData& td + ); + + //- Influence of edge on face. + template<class TrackingData> + inline bool updateFace + ( + const polyMesh& mesh, + const primitivePatch& patch, + const label faceI, + const label edgeI, + const patchEdgeFaceInfo& edgeInfo, + const scalar tol, + TrackingData& td + ); + + //- Same (like operator==) + template<class TrackingData> + inline bool equal(const patchEdgeFaceInfo&, TrackingData& td) const; + + + // Member Operators + + // Needed for List IO + inline bool operator==(const patchEdgeFaceInfo&) const; + inline bool operator!=(const patchEdgeFaceInfo&) const; + + + // IOstream Operators + + friend Ostream& operator<<(Ostream&, const patchEdgeFaceInfo&); + friend Istream& operator>>(Istream&, patchEdgeFaceInfo&); +}; + + +//- Data associated with patchEdgeFaceInfo type are contiguous +template<> +inline bool contiguous<patchEdgeFaceInfo>() +{ + return true; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "patchEdgeFaceInfoI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfoI.H b/src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfoI.H new file mode 100644 index 0000000000000000000000000000000000000000..881e586626a1f69760486c948a9a66ed0da341dd --- /dev/null +++ b/src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfoI.H @@ -0,0 +1,268 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "polyMesh.H" +#include "transform.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +// Update this with w2 if w2 nearer to pt. +template<class TrackingData> +inline bool Foam::patchEdgeFaceInfo::update +( + const point& pt, + const patchEdgeFaceInfo& w2, + const scalar tol, + TrackingData& td +) +{ + scalar dist2 = magSqr(pt - w2.origin()); + + if (!valid(td)) + { + // current not yet set so use any value + distSqr_ = dist2; + origin_ = w2.origin(); + + return true; + } + + scalar diff = distSqr_ - dist2; + + if (diff < 0) + { + // already nearer to pt + return false; + } + + if ((diff < SMALL) || ((distSqr_ > SMALL) && (diff/distSqr_ < tol))) + { + // don't propagate small changes + return false; + } + else + { + // update with new values + distSqr_ = dist2; + origin_ = w2.origin(); + + return true; + } +} + + +// Update this with w2 (information on same edge) +template<class TrackingData> +inline bool Foam::patchEdgeFaceInfo::update +( + const patchEdgeFaceInfo& w2, + const scalar tol, + TrackingData& td +) +{ + if (!valid(td)) + { + // current not yet set so use any value + distSqr_ = w2.distSqr(); + origin_ = w2.origin(); + + return true; + } + + scalar diff = distSqr_ - w2.distSqr(); + + if (diff < 0) + { + // already nearer to pt + return false; + } + + if ((diff < SMALL) || ((distSqr_ > SMALL) && (diff/distSqr_ < tol))) + { + // don't propagate small changes + return false; + } + else + { + // update with new values + distSqr_ = w2.distSqr(); + origin_ = w2.origin(); + + return true; + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Null constructor +inline Foam::patchEdgeFaceInfo::patchEdgeFaceInfo() +: + origin_(point::max), + distSqr_(sqr(GREAT)) +{} + + +// Construct from origin, distance +inline Foam::patchEdgeFaceInfo::patchEdgeFaceInfo +( + const point& origin, + const scalar distSqr +) +: + origin_(origin), + distSqr_(distSqr) +{} + + +// Construct as copy +inline Foam::patchEdgeFaceInfo::patchEdgeFaceInfo(const patchEdgeFaceInfo& wpt) +: + origin_(wpt.origin()), + distSqr_(wpt.distSqr()) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline const Foam::point& Foam::patchEdgeFaceInfo::origin() const +{ + return origin_; +} + + +inline Foam::scalar Foam::patchEdgeFaceInfo::distSqr() const +{ + return distSqr_; +} + + +template<class TrackingData> +inline bool Foam::patchEdgeFaceInfo::valid(TrackingData& td) const +{ + return origin_ != point::max; +} + + +template<class TrackingData> +inline void Foam::patchEdgeFaceInfo::transform +( + const polyMesh& mesh, + const primitivePatch& patch, + const tensor& rotTensor, + const scalar tol, + TrackingData& td +) +{ + origin_ = Foam::transform(rotTensor, origin_); +} + + +template<class TrackingData> +inline bool Foam::patchEdgeFaceInfo::updateEdge +( + const polyMesh& mesh, + const primitivePatch& patch, + const label edgeI, + const label faceI, + const patchEdgeFaceInfo& faceInfo, + const scalar tol, + TrackingData& td +) +{ + const edge& e = patch.edges()[edgeI]; + point eMid = + 0.5 + * ( + patch.points()[patch.meshPoints()[e[0]]] + + patch.points()[patch.meshPoints()[e[1]]] + ); + return update(eMid, faceInfo, tol, td); +} + + +template<class TrackingData> +inline bool Foam::patchEdgeFaceInfo::updateEdge +( + const polyMesh& mesh, + const primitivePatch& patch, + const patchEdgeFaceInfo& edgeInfo, + const bool sameOrientation, + const scalar tol, + TrackingData& td +) +{ + return update(edgeInfo, tol, td); +} + + +template<class TrackingData> +inline bool Foam::patchEdgeFaceInfo::updateFace +( + const polyMesh& mesh, + const primitivePatch& patch, + const label faceI, + const label edgeI, + const patchEdgeFaceInfo& edgeInfo, + const scalar tol, + TrackingData& td +) +{ + return update(patch.faceCentres()[faceI], edgeInfo, tol, td); +} + + +template <class TrackingData> +inline bool Foam::patchEdgeFaceInfo::equal +( + const patchEdgeFaceInfo& rhs, + TrackingData& td +) const +{ + return operator==(rhs); +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +inline bool Foam::patchEdgeFaceInfo::operator== +( + const Foam::patchEdgeFaceInfo& rhs +) const +{ + return origin() == rhs.origin(); +} + + +inline bool Foam::patchEdgeFaceInfo::operator!= +( + const Foam::patchEdgeFaceInfo& rhs +) const +{ + return !(*this == rhs); +} + + +// ************************************************************************* // diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H index 68c432e97aadc3be08e4cdfea6b34cabe048b8ff..b2f0f4eb92e1764521eacfb57897a184411d9bdd 100644 --- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H +++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H @@ -339,6 +339,39 @@ public: } + //- Wrapper around map/interpolate data distribution with supplied op + template<class Type, class BinaryOp> + void distribute(List<Type>& lst, const BinaryOp& bop) const + { + switch (mode_) + { + case NEARESTPATCHFACEAMI: + { + lst = AMI().interpolateToSource + ( + Field<Type>(lst.xfer()), + bop + ); + break; + } + default: + { + map().distribute + ( + Pstream::defaultCommsType, + map().schedule(), + map().constructSize(), + map().subMap(), + map().constructMap(), + lst, + bop, + pTraits<Type>::zero + ); + } + } + } + + //- Wrapper around map/interpolate data distribution template<class Type> void reverseDistribute(List<Type>& lst) const @@ -359,6 +392,40 @@ public: } + //- Wrapper around map/interpolate data distribution with supplied op + template<class Type, class BinaryOp> + void reverseDistribute(List<Type>& lst, const BinaryOp& bop) const + { + switch (mode_) + { + case NEARESTPATCHFACEAMI: + { + lst = AMI().interpolateToTarget + ( + Field<Type>(lst.xfer()), + bop + ); + break; + } + default: + { + label cSize = patch_.size(); + map().distribute + ( + Pstream::defaultCommsType, + map().schedule(), + cSize, + map().constructMap(), + map().subMap(), + lst, + bop, + pTraits<Type>::zero + ); + } + } + } + + //- Return reference to the parallel distribution map const mapDistribute& map() const { diff --git a/src/postProcessing/functionObjects/field/Make/files b/src/postProcessing/functionObjects/field/Make/files index c1497c5e5484b0fa8ac62e29f7afdbdd0fbe4e05..9d465d2663337a86c4472089e6cd2ad239373e38 100644 --- a/src/postProcessing/functionObjects/field/Make/files +++ b/src/postProcessing/functionObjects/field/Make/files @@ -3,6 +3,9 @@ fieldAverage/fieldAverageItem/fieldAverageItem.C fieldAverage/fieldAverageItem/fieldAverageItemIO.C fieldAverage/fieldAverageFunctionObject/fieldAverageFunctionObject.C +fieldCoordinateSystemTransform/fieldCoordinateSystemTransform.C +fieldCoordinateSystemTransform/fieldCoordinateSystemTransformFunctionObject.C + fieldMinMax/fieldMinMax.C fieldMinMax/fieldMinMaxFunctionObject.C diff --git a/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/IOfieldCoordinateSystemTransform.H b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/IOfieldCoordinateSystemTransform.H new file mode 100644 index 0000000000000000000000000000000000000000..36216f399bf18258482528045045a7a11f340afb --- /dev/null +++ b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/IOfieldCoordinateSystemTransform.H @@ -0,0 +1,50 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Typedef + Foam::IOfieldCoordinateSystemTransform + +Description + Instance of the generic IOOutputFilter for fieldCoordinateSystemTransform. + +\*---------------------------------------------------------------------------*/ + +#ifndef IOfieldCoordinateSystemTransform_H +#define IOfieldCoordinateSystemTransform_H + +#include "fieldCoordinateSystemTransform.H" +#include "IOOutputFilter.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + typedef IOOutputFilter<fieldCoordinateSystemTransform> + IOfieldCoordinateSystemTransform; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransform.C b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransform.C new file mode 100644 index 0000000000000000000000000000000000000000..afb592bc78dab27c912d39f5fdab4a760b072e87 --- /dev/null +++ b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransform.C @@ -0,0 +1,118 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "fieldCoordinateSystemTransform.H" +#include "dictionary.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(Foam::fieldCoordinateSystemTransform, 0); + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::fieldCoordinateSystemTransform::fieldCoordinateSystemTransform +( + const word& name, + const objectRegistry& obr, + const dictionary& dict, + const bool loadFromFiles +) +: + name_(name), + obr_(obr), + active_(true), + fieldSet_(), + coordSys_(dict, obr) +{ + // Check if the available mesh is an fvMesh otherise deactivate + if (!isA<fvMesh>(obr_)) + { + active_ = false; + WarningIn + ( + "fieldCoordinateSystemTransform::fieldCoordinateSystemTransform" + "(" + "const word&, " + "const objectRegistry&, " + "const dictionary&, " + "const bool" + ")" + ) << "No fvMesh available, deactivating." + << endl; + } + + read(dict); + + Info<< type() << ":" << nl + << " Applying transformation from global Cartesian to local " + << coordSys_ << nl << endl; +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::fieldCoordinateSystemTransform::~fieldCoordinateSystemTransform() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::fieldCoordinateSystemTransform::read(const dictionary& dict) +{ + if (active_) + { + dict.lookup("fields") >> fieldSet_; + } +} + + +void Foam::fieldCoordinateSystemTransform::execute() +{ + // Do nothing +} + + +void Foam::fieldCoordinateSystemTransform::end() +{ + // Do nothing +} + + +void Foam::fieldCoordinateSystemTransform::write() +{ + forAll(fieldSet_, fieldI) + { + // If necessary load field + transform<scalar>(fieldSet_[fieldI]); + transform<vector>(fieldSet_[fieldI]); + transform<sphericalTensor>(fieldSet_[fieldI]); + transform<symmTensor>(fieldSet_[fieldI]); + transform<tensor>(fieldSet_[fieldI]); + } +} + + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransform.H b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransform.H new file mode 100644 index 0000000000000000000000000000000000000000..138d2d09d5fbe65cc6b424de9aa5337b441ff6df --- /dev/null +++ b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransform.H @@ -0,0 +1,163 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Class + Foam::fieldCoordinateSystemTransform + +Description + Transforms fields from global cartesian co-ordinates to local co-ordinate + system + +SourceFiles + fieldCoordinateSystemTransform.C + IOfieldCoordinateSystemTransform.H + +\*---------------------------------------------------------------------------*/ + +#ifndef fieldCoordinateSystemTransform_H +#define fieldCoordinateSystemTransform_H + +#include "OFstream.H" +#include "pointFieldFwd.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "coordinateSystem.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class objectRegistry; +class dictionary; +class mapPolyMesh; + +/*---------------------------------------------------------------------------*\ + Class fieldCoordinateSystemTransform Declaration +\*---------------------------------------------------------------------------*/ + +class fieldCoordinateSystemTransform +{ +protected: + + // Protected data + + //- Name + word name_; + + const objectRegistry& obr_; + + //- on/off switch + bool active_; + + //- Fields to transform + wordList fieldSet_; + + //- Co-ordinate system to transform to + coordinateSystem coordSys_; + + + // Protected Member Functions + + //- Disallow default bitwise copy construct + fieldCoordinateSystemTransform(const fieldCoordinateSystemTransform&); + + //- Disallow default bitwise assignment + void operator=(const fieldCoordinateSystemTransform&); + + template<class Type> + void transform(const word& fieldName) const; + + template<class Type> + void transformField(const Type& field) const; + + +public: + + //- Runtime type information + TypeName("fieldCoordinateSystemTransform"); + + + // Constructors + + //- Construct for given objectRegistry and dictionary. + // Allow the possibility to load fields from files + fieldCoordinateSystemTransform + ( + const word& name, + const objectRegistry&, + const dictionary&, + const bool loadFromFiles = false + ); + + + //- Destructor + virtual ~fieldCoordinateSystemTransform(); + + + // Member Functions + + //- Return name of the fieldCoordinateSystemTransform object + virtual const word& name() const + { + return name_; + } + + //- Read the input data + virtual void read(const dictionary&); + + //- Execute, currently does nothing + virtual void execute(); + + //- Execute at the final time-loop, currently does nothing + virtual void end(); + + //- Write + virtual void write(); + + //- Update for changes of mesh + virtual void updateMesh(const mapPolyMesh&) + {} + + //- Update for changes of mesh + virtual void movePoints(const pointField&) + {} +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "fieldCoordinateSystemTransformTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformFunctionObject.C similarity index 55% rename from applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C rename to src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformFunctionObject.C index ea8f509a292d1761a0e9665f17153aa0171c87fb..be3f5624b474da841e9ddb374c0aed611ecf4a16 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C +++ b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformFunctionObject.C @@ -23,45 +23,23 @@ License \*---------------------------------------------------------------------------*/ -#include "patchPointEdgeCirculator.H" +#include "fieldCoordinateSystemTransformFunctionObject.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -const Foam::patchPointEdgeCirculator -Foam::patchPointEdgeCirculator::endConstIter -( - *reinterpret_cast<primitiveFacePatch*>(0), // primitiveFacePatch - *reinterpret_cast<PackedBoolList*>(0), // PackedBoolList - -1, // edgeID - -1, // index - -1 // pointID -); - - -Foam::Ostream& Foam::operator<< -( - Ostream& os, - const InfoProxy<patchPointEdgeCirculator>& ip -) +namespace Foam { - const patchPointEdgeCirculator& c = ip.t_; - - const pointField& pts = c.patch_.localPoints(); - const edge& e = c.patch_.edges()[c.edgeID_]; - label faceI = c.faceID(); - - os << "around point:" << c.pointID_ - << " coord:" << pts[c.pointID_] - << " at edge:" << c.edgeID_ - << " coord:" << pts[e.otherVertex(c.pointID_)] - << " in direction of face:" << faceI; - - if (faceI != -1) - { - os << " fc:" << c.patch_.localFaces()[faceI].centre(pts); - } - return os; + defineNamedTemplateTypeNameAndDebug + ( + fieldCoordinateSystemTransformFunctionObject, 0 + ); + + addToRunTimeSelectionTable + ( + functionObject, + fieldCoordinateSystemTransformFunctionObject, + dictionary + ); } - // ************************************************************************* // diff --git a/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformFunctionObject.H b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformFunctionObject.H new file mode 100644 index 0000000000000000000000000000000000000000..20cf0fccbf9b73b590cfe41b5d61af920162c6ab --- /dev/null +++ b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformFunctionObject.H @@ -0,0 +1,54 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +Typedef + Foam::fieldCoordinateSystemTransformFunctionObject + +Description + FunctionObject wrapper around fieldCoordinateSystemTransform to allow + them to be created via the functions entry within controlDict. + +SourceFiles + fieldCoordinateSystemTransformFunctionObject.C + +\*---------------------------------------------------------------------------*/ + +#ifndef fieldCoordinateSystemTransformFunctionObject_H +#define fieldCoordinateSystemTransformFunctionObject_H + +#include "fieldCoordinateSystemTransform.H" +#include "OutputFilterFunctionObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + typedef OutputFilterFunctionObject<fieldCoordinateSystemTransform> + fieldCoordinateSystemTransformFunctionObject; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformTemplates.C b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..c1ecc72707bce985a0705929006be29edaf1a993 --- /dev/null +++ b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformTemplates.C @@ -0,0 +1,158 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "fieldCoordinateSystemTransform.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "Time.H" +#include "transformGeometricField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<class Type> +void Foam::fieldCoordinateSystemTransform::transformField +( + const Type& field +) const +{ + const word& fieldName = field.name() + "Transformed"; + + dimensionedTensor R("R", field.dimensions(), coordSys_.R()); + + if (obr_.foundObject<Type>(fieldName)) + { + Type& transField = + const_cast<Type&>(obr_.lookupObject<Type>(fieldName)); + + transField == field; + + forAll(field, i) + { + Foam::transform(transField, R, transField); + } + + transField.write(); + } + else + { + Type& transField = obr_.store + ( + new Type + ( + IOobject + ( + fieldName, + obr_.time().timeName(), + obr_, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + field + ) + ); + + forAll(field, i) + { + Foam::transform(transField, R, transField); + } + + transField.write(); + } +} + + +template<class Type> +void Foam::fieldCoordinateSystemTransform::transform +( + const word& fieldName +) const +{ + typedef GeometricField<Type, fvPatchField, volMesh> vfType; + typedef GeometricField<Type, fvsPatchField, surfaceMesh> sfType; + + if (obr_.foundObject<vfType>(fieldName)) + { + if (debug) + { + Info<< type() << ": Field " << fieldName << " already in database" + << endl; + } + + transformField<vfType>(obr_.lookupObject<vfType>(fieldName)); + } + else if (obr_.foundObject<sfType>(fieldName)) + { + if (debug) + { + Info<< type() << ": Field " << fieldName << " already in database" + << endl; + } + + transformField<sfType>(obr_.lookupObject<sfType>(fieldName)); + } + else + { + IOobject fieldHeader + ( + fieldName, + obr_.time().timeName(), + obr_, + IOobject::MUST_READ, + IOobject::NO_WRITE + ); + + if + ( + fieldHeader.headerOk() + && fieldHeader.headerClassName() == vfType::typeName + ) + { + if (debug) + { + Info<< type() << ": Field " << fieldName << " read from file" + << endl; + } + + transformField<vfType>(obr_.lookupObject<vfType>(fieldName)); + } + else if + ( + fieldHeader.headerOk() + && fieldHeader.headerClassName() == sfType::typeName + ) + { + if (debug) + { + Info<< type() << ": Field " << fieldName << " read from file" + << endl; + } + + transformField<sfType>(obr_.lookupObject<sfType>(fieldName)); + } + } +} + + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/postProcessingDict b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/postProcessingDict new file mode 100644 index 0000000000000000000000000000000000000000..e0d17f6271ab7d9e767ffe454713e2c917eaae7e --- /dev/null +++ b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/postProcessingDict @@ -0,0 +1,50 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object postProcessingDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +functions +{ + fieldCoordinateSystemTransform1 + { + // Type of functionObject + type fieldCoordinateSystemTransform; + + // Where to load it from (if not already in solver) + functionObjectLibs ("libfieldCoordinateSystemTransform.so"); + + // Function object enabled flag + enabled true; + + // When to output the average fields + outputControl outputTime; + + // Fields to be transformed - runTime modifiable + fields + ( + U + UMean + UPrime2Mean + ); + + coordinateSystem + { + origin (0.001 0 0); + e1 (1 0.15 0); + e3 (0 0 -1); + } + } +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/field/streamLine/controlDict b/src/postProcessing/functionObjects/field/streamLine/controlDict index 048b7af4f3376fa0e69e9c740107315d06d7a58f..4d51dd82e1f758f5af36f68e927088bb87b14d90 100644 --- a/src/postProcessing/functionObjects/field/streamLine/controlDict +++ b/src/postProcessing/functionObjects/field/streamLine/controlDict @@ -60,7 +60,7 @@ functions setFormat vtk; //gnuplot; //xmgr; //raw; //jplot; // Velocity field to use for tracking. - U U; + UName U; // Interpolation method. Default is cellPoint. See sampleDict. //interpolationScheme pointMVC; diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLine.C b/src/postProcessing/functionObjects/field/streamLine/streamLine.C index ae902f9c54f1ab1ba50713fc15cc795835b95c48..03a5b8b539c4b803a9b050dda20a0f9cb3e43ef6 100644 --- a/src/postProcessing/functionObjects/field/streamLine/streamLine.C +++ b/src/postProcessing/functionObjects/field/streamLine/streamLine.C @@ -318,7 +318,22 @@ void Foam::streamLine::read(const dictionary& dict) { //dict_ = dict; dict.lookup("fields") >> fields_; - UName_ = dict.lookupOrDefault<word>("U", "U"); + if (dict.found("UName")) + { + dict.lookup("UName") >> UName_; + } + else + { + UName_ = "U"; + if (dict.found("U")) + { + IOWarningIn("streamLine::read(const dictionary&)", dict) + << "Using deprecated entry \"U\"." + << " Please use \"UName\" instead." + << endl; + dict.lookup("U") >> UName_; + } + } dict.lookup("trackForward") >> trackForward_; dict.lookup("lifeTime") >> lifeTime_; if (lifeTime_ < 1) diff --git a/src/regionModels/regionModel/regionModel/regionModel.H b/src/regionModels/regionModel/regionModel/regionModel.H index 2af98421d108216d76b5606f0ff7d3bf4abc28cc..49f708bcd81e38562234da286387e11e3e5f0dfe 100644 --- a/src/regionModels/regionModel/regionModel/regionModel.H +++ b/src/regionModels/regionModel/regionModel/regionModel.H @@ -231,6 +231,24 @@ public: List<Type>& primaryFieldField ) const; + //- Convert a local region field to the primary region with op + template<class Type, class BinaryOp> + void toPrimary + ( + const label regionPatchI, + List<Type>& regionField, + const BinaryOp& bop + ) const; + + //- Convert a primary region field to the local region with op + template<class Type, class BinaryOp> + void toRegion + ( + const label regionPatchI, + List<Type>& primaryFieldField, + const BinaryOp& bop + ) const; + // Evolution diff --git a/src/regionModels/regionModel/regionModel/regionModelTemplates.C b/src/regionModels/regionModel/regionModel/regionModelTemplates.C index 977da7dd8167df5a3b5acf703cc09f434602c5b6..ab013ad8e08f09cff6df2b8142805c7605813b5f 100644 --- a/src/regionModels/regionModel/regionModel/regionModelTemplates.C +++ b/src/regionModels/regionModel/regionModel/regionModelTemplates.C @@ -77,4 +77,69 @@ void Foam::regionModels::regionModel::toRegion } +template<class Type, class BinaryOp> +void Foam::regionModels::regionModel::toPrimary +( + const label regionPatchI, + List<Type>& regionField, + const BinaryOp& bop +) const +{ + forAll(intCoupledPatchIDs_, i) + { + if (intCoupledPatchIDs_[i] == regionPatchI) + { + const mappedPatchBase& mpb = + refCast<const mappedPatchBase> + ( + regionMesh().boundaryMesh()[regionPatchI] + ); + mpb.reverseDistribute(regionField, bop); + return; + } + } + + FatalErrorIn + ( + "const void toPrimary" + "(" + "const label, " + "List<Type>&, " + "const BinaryOp&" + ") const" + ) << "Region patch ID " << regionPatchI << " not found in region mesh" + << abort(FatalError); +} + + +template<class Type, class BinaryOp> +void Foam::regionModels::regionModel::toRegion +( + const label regionPatchI, + List<Type>& primaryField, + const BinaryOp& bop +) const +{ + forAll(intCoupledPatchIDs_, i) + { + if (intCoupledPatchIDs_[i] == regionPatchI) + { + const mappedPatchBase& mpb = + refCast<const mappedPatchBase> + ( + regionMesh().boundaryMesh()[regionPatchI] + ); + mpb.distribute(primaryField, bop); + return; + } + } + + FatalErrorIn + ( + "const void toRegion(const label, List<Type>&, const BinaryOp&) const" + ) << "Region patch ID " << regionPatchI << " not found in region mesh" + << abort(FatalError); +} + + // ************************************************************************* // diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C index 4fede11db6d5d5ea981b29856360fda57d7f15be..2a396388c816392b2a30f7f06b01bd9787bcf0b5 100644 --- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C +++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C @@ -54,10 +54,6 @@ externalWallHeatFluxTemperatureFvPatchScalarField::operationModeNames; } // End namespace Foam - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::externalWallHeatFluxTemperatureFvPatchScalarField:: @@ -118,7 +114,7 @@ externalWallHeatFluxTemperatureFvPatchScalarField oldMode_ = fixedHeatFlux; q_ = scalarField("q", dict, p.size()); } - else if(dict.found("h") && dict.found("Ta") && !dict.found("q")) + else if (dict.found("h") && dict.found("Ta") && !dict.found("q")) { oldMode_ = fixedHeatTransferCoeff; h_ = scalarField("h", dict, p.size()); @@ -209,7 +205,7 @@ void Foam::externalWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs() { q = q_; } - else if(oldMode_ == fixedHeatTransferCoeff) + else if (oldMode_ == fixedHeatTransferCoeff) { q = (Ta_ - *this)*h_; } diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/streamLines b/tutorials/incompressible/simpleFoam/motorBike/system/streamLines index 1f9488d8c2fce94156d99abf5403eb02936aa99b..82c3f1d68c2c92f7c1c55f84dbf15627305ffc88 100644 --- a/tutorials/incompressible/simpleFoam/motorBike/system/streamLines +++ b/tutorials/incompressible/simpleFoam/motorBike/system/streamLines @@ -17,7 +17,7 @@ streamLines setFormat vtk; //gnuplot; //xmgr; //raw; //jplot; // Velocity field to use for tracking. - U U; + UName U; // Tracked forwards (+U) or backwards (-U) trackForward true; diff --git a/tutorials/incompressible/simpleFoam/pitzDaily/system/controlDict b/tutorials/incompressible/simpleFoam/pitzDaily/system/controlDict index a3baf4a417fd5498ff9d42a15a57f0fbd1ea5b38..aaedd23dab9b1b592cc799596df2f301c0917407 100644 --- a/tutorials/incompressible/simpleFoam/pitzDaily/system/controlDict +++ b/tutorials/incompressible/simpleFoam/pitzDaily/system/controlDict @@ -61,7 +61,7 @@ functions setFormat vtk; //gnuplot; //xmgr; //raw; //jplot; // Velocity field to use for tracking. - U U; + UName U; // Tracked forwards (+U) or backwards (-U) trackForward true; diff --git a/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/system/controlDict b/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/system/controlDict index 89dd2ff19ab190e6f0dafddc3c2b615d0bb5e2f0..c2b5d2e493f326b7c385fc98ee0957a1d3ea4427 100644 --- a/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/system/controlDict +++ b/tutorials/incompressible/simpleFoam/pitzDailyExptInlet/system/controlDict @@ -61,7 +61,7 @@ functions setFormat vtk; //gnuplot; //xmgr; //raw; //jplot; // Velocity field to use for tracking. - U U; + UName U; // Tracked forwards (+U) or backwards (-U) trackForward true;