diff --git a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C index 7950d6a0068e0b31dee4137e9b09d0154eaa9fa8..557c8aa1fb79d3f7c50e5a20dd9307dc9800ab06 100644 --- a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C +++ b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C @@ -66,6 +66,7 @@ Description #include "Switch.H" #include "bound.H" #include "dynamicRefineFvMesh.H" +#include "pimpleControl.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -83,6 +84,8 @@ int main(int argc, char *argv[]) #include "compressibleCourantNo.H" #include "setInitialDeltaT.H" + pimpleControl pimple(mesh); + scalar StCoNum = 0.0; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -94,7 +97,6 @@ int main(int argc, char *argv[]) while (runTime.run()) { #include "readTimeControls.H" - #include "readPISOControls.H" #include "compressibleCourantNo.H" #include "setDeltaT.H" @@ -163,26 +165,35 @@ int main(int argc, char *argv[]) #include "rhoEqn.H" - #include "UEqn.H" - // --- PISO loop - for (int corr=1; corr<=nCorr; corr++) + // --- Pressure-velocity PIMPLE corrector loop + for (pimple.start(); pimple.loop(); pimple++) { - #include "bEqn.H" - #include "ftEqn.H" - #include "huEqn.H" - #include "hEqn.H" + #include "UEqn.H" + - if (!ign.ignited()) + // --- PISO loop + for (int corr=1; corr<=pimple.nCorr(); corr++) { - hu == h; + #include "bEqn.H" + #include "ftEqn.H" + #include "huEqn.H" + #include "hEqn.H" + + if (!ign.ignited()) + { + hu == h; + } + + #include "pEqn.H" } - #include "pEqn.H" + if (pimple.turbCorr()) + { + turbulence->correct(); + } } - turbulence->correct(); - runTime.write(); Info<< "\nExecutionTime = " 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/etc/config/settings.csh b/etc/config/settings.csh index 93f3e506abdff716446fbfd173f00d76a5e42fb6..2b0a132ebca31976263356d13850deae79525b33 100644 --- a/etc/config/settings.csh +++ b/etc/config/settings.csh @@ -203,10 +203,6 @@ case ThirdParty: switch ("$WM_COMPILER") case Gcc: case Gcc++0x: - set gcc_version=gcc-4.4.3 - set gmp_version=gmp-5.0.1 - set mpfr_version=mpfr-2.4.2 - breaksw case Gcc46: case Gcc46++0x: set gcc_version=gcc-4.6.1 @@ -487,9 +483,31 @@ case QSMPI: breaksw case SGIMPI: - setenv FOAM_MPI ${MPI_ROOT##*/} + if ( ! $?MPI_ROOT) setenv MPI_ROOT /dummy + + if ( ! -d "$MPI_ROOT" ) then + echo "Warning in $WM_PROJECT_DIR/etc/config/settings.csh:" + echo " MPI_ROOT not a valid mpt installation directory." + echo " Please set MPI_ROOT to the mpt installation directory." + echo " (usually done by loading the mpt module)" + echo " MPI_ROOT currently set to '$MPI_ROOT'" + endif + + if ( "${MPI_ROOT:h}/" == $MPI_ROOT ) then + setenv MPI_ROOT ${MPI_ROOT:h} + endif + + setenv FOAM_MPI ${MPI_ROOT:t} setenv MPI_ARCH_PATH $MPI_ROOT + + if ($?FOAM_VERBOSE && $?prompt) then + echo "Using SGI MPT:" + echo " MPI_ROOT : $MPI_ROOT" + echo " FOAM_MPI : $FOAM_MPI" + endif + + _foamAddPath $MPI_ARCH_PATH/bin _foamAddLib $MPI_ARCH_PATH/lib breaksw diff --git a/etc/config/settings.sh b/etc/config/settings.sh index 5111ef405a3152cf02d526e1322127e6a346e0a3..e56afffa179418e1b97ff2d6205c6eade342526e 100644 --- a/etc/config/settings.sh +++ b/etc/config/settings.sh @@ -222,12 +222,7 @@ fi case "${foamCompiler}" in OpenFOAM | ThirdParty) case "$WM_COMPILER" in - Gcc | Gcc++0x) - gcc_version=gcc-4.4.3 - gmp_version=gmp-5.0.1 - mpfr_version=mpfr-2.4.2 - ;; - Gcc46 | Gcc46++0x) + Gcc | Gcc++0x | Gcc46 | Gcc46++0x) gcc_version=gcc-4.6.1 gmp_version=gmp-5.0.2 mpfr_version=mpfr-3.0.1 @@ -511,9 +506,30 @@ QSMPI) ;; SGIMPI) + lastCharID=$(( ${#MPI_ROOT} - 1 )) + if [ "${MPI_ROOT:$lastCharID:1}" == '/' ] + then + MPI_ROOT=${MPI_ROOT:0:$lastCharID} + fi + export FOAM_MPI=${MPI_ROOT##*/} export MPI_ARCH_PATH=$MPI_ROOT + if [ ! -d "$MPI_ROOT" -o -z "$MPI_ARCH_PATH" ] + then + echo "Warning in $WM_PROJECT_DIR/etc/config/settings.sh:" 1>&2 + echo " MPI_ROOT not a valid mpt installation directory or ending in a '/'." 1>&2 + echo " Please set MPI_ROOT to the mpt installation directory." 1>&2 + echo " MPI_ROOT currently set to '$MPI_ROOT'" 1>&2 + fi + + if [ "$FOAM_VERBOSE" -a "$PS1" ] + then + echo "Using SGI MPT:" + echo " MPI_ROOT : $MPI_ROOT" + echo " FOAM_MPI : $FOAM_MPI" + fi + _foamAddPath $MPI_ARCH_PATH/bin _foamAddLib $MPI_ARCH_PATH/lib ;; diff --git a/etc/controlDict b/etc/controlDict index c367f5d85329c3751b5297a6ae88008d9d531b66..7a14e9036c0c1e10af492423710bf0cf488fdbf0 100644 --- a/etc/controlDict +++ b/etc/controlDict @@ -58,6 +58,11 @@ OptimisationSwitches commsType nonBlocking; //scheduled; //blocking; floatTransfer 0; nProcsSimpleSum 0; + + // Force dumping (at next timestep) upon signal (-1 to disable) + writeNowSignal -1; //10; + // Force dumping (at next timestep) upon signal (-1 to disable) and exit + stopAtWriteNowSignal -1; } diff --git a/src/OSspecific/POSIX/Make/files b/src/OSspecific/POSIX/Make/files index c7396b63452db5569a2c52116371d2fa3e50aaba..90dc5bc92ec2924ed922296b620d8be742091f91 100644 --- a/src/OSspecific/POSIX/Make/files +++ b/src/OSspecific/POSIX/Make/files @@ -2,6 +2,8 @@ signals/sigFpe.C signals/sigSegv.C signals/sigInt.C signals/sigQuit.C +signals/sigStopAtWriteNow.C +signals/sigWriteNow.C regExp.C timer.C fileStat.C diff --git a/src/OSspecific/POSIX/POSIX.C b/src/OSspecific/POSIX/POSIX.C index 7ab46939cdbe5b08aaacd6f64eca01aaf5392254..e214e7eecabcf09054838adbc52bfe67ab84a420 100644 --- a/src/OSspecific/POSIX/POSIX.C +++ b/src/OSspecific/POSIX/POSIX.C @@ -123,7 +123,7 @@ bool Foam::setEnv } -Foam::word Foam::hostName(bool full) +Foam::string Foam::hostName(bool full) { char buf[128]; ::gethostname(buf, sizeof(buf)); @@ -142,7 +142,7 @@ Foam::word Foam::hostName(bool full) } -Foam::word Foam::domainName() +Foam::string Foam::domainName() { char buf[128]; ::gethostname(buf, sizeof(buf)); @@ -159,11 +159,11 @@ Foam::word Foam::domainName() } } - return word::null; + return string::null; } -Foam::word Foam::userName() +Foam::string Foam::userName() { struct passwd* pw = ::getpwuid(::getuid()); @@ -173,7 +173,7 @@ Foam::word Foam::userName() } else { - return word::null; + return string::null; } } @@ -209,7 +209,7 @@ Foam::fileName Foam::home() } -Foam::fileName Foam::home(const word& userName) +Foam::fileName Foam::home(const string& userName) { struct passwd* pw; @@ -1069,7 +1069,7 @@ void Foam::fdClose(const int fd) bool Foam::ping ( - const word& destName, + const string& destName, const label destPort, const label timeOut ) @@ -1083,7 +1083,7 @@ bool Foam::ping { FatalErrorIn ( - "Foam::ping(const word&, ...)" + "Foam::ping(const string&, ...)" ) << "gethostbyname error " << h_errno << " for host " << destName << abort(FatalError); } @@ -1097,7 +1097,7 @@ bool Foam::ping { FatalErrorIn ( - "Foam::ping(const word&, const label)" + "Foam::ping(const string&, const label)" ) << "socket error" << abort(FatalError); } @@ -1149,7 +1149,7 @@ bool Foam::ping } -bool Foam::ping(const word& hostname, const label timeOut) +bool Foam::ping(const string& hostname, const label timeOut) { return ping(hostname, 222, timeOut) || ping(hostname, 22, timeOut); } diff --git a/src/OSspecific/POSIX/signals/sigStopAtWriteNow.C b/src/OSspecific/POSIX/signals/sigStopAtWriteNow.C new file mode 100644 index 0000000000000000000000000000000000000000..e973bc1bd8fc90039bb43542cf5b3a7dce8000b3 --- /dev/null +++ b/src/OSspecific/POSIX/signals/sigStopAtWriteNow.C @@ -0,0 +1,155 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "sigStopAtWriteNow.H" +#include "error.H" +#include "JobInfo.H" +#include "IOstreams.H" +#include "Time.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +// Signal number to catch +int Foam::sigStopAtWriteNow::signal_ +( + debug::optimisationSwitch("stopAtWriteNowSignal", -1) +); + +static Foam::Time const* runTimePtr_ = NULL; + + +struct sigaction Foam::sigStopAtWriteNow::oldAction_; + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::sigStopAtWriteNow::sigHandler(int) +{ + // Reset old handling + if (sigaction(signal_, &oldAction_, NULL) < 0) + { + FatalErrorIn + ( + "Foam::sigStopAtWriteNow::sigHandler(int)" + ) << "Cannot reset " << signal_ << " trapping" + << abort(FatalError); + } + + // Update jobInfo file + jobInfo.signalEnd(); + + Info<< "sigStopAtWriteNow :" + << " setting up write and stop at end of the next iteration" + << nl << endl; + runTimePtr_->stopAt(Time::saWriteNow); + + //// Throw signal (to old handler) + //raise(signal_); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sigStopAtWriteNow::sigStopAtWriteNow(){} + + +Foam::sigStopAtWriteNow::sigStopAtWriteNow +( + const bool verbose, + const Time& runTime +) +{ + 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; + + struct sigaction newAction; + newAction.sa_handler = sigHandler; + newAction.sa_flags = SA_NODEFER; + sigemptyset(&newAction.sa_mask); + if (sigaction(signal_, &newAction, &oldAction_) < 0) + { + FatalErrorIn + ( + "Foam::sigStopAtWriteNow::sigStopAtWriteNow" + "(const bool, const Time&)" + ) << "Cannot set " << signal_ << " trapping" + << abort(FatalError); + } + + if (verbose) + { + Info<< "sigStopAtWriteNow :" + << " Enabling writing and stopping upon signal " << signal_ + << endl; + } + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sigStopAtWriteNow::~sigStopAtWriteNow() +{ + // Reset old handling + if (signal_ > 0) + { + if (sigaction(signal_, &oldAction_, NULL) < 0) + { + FatalErrorIn + ( + "Foam::sigStopAtWriteNow::~sigStopAtWriteNow()" + ) << "Cannot reset " << signal_ << " trapping" + << abort(FatalError); + } + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::sigStopAtWriteNow::active() const +{ + return signal_ > 0; +} + + +// ************************************************************************* // diff --git a/src/OSspecific/POSIX/signals/sigStopAtWriteNow.H b/src/OSspecific/POSIX/signals/sigStopAtWriteNow.H new file mode 100644 index 0000000000000000000000000000000000000000..4c07248eb21be8c29325914f1dcbc99bb72c86d0 --- /dev/null +++ b/src/OSspecific/POSIX/signals/sigStopAtWriteNow.H @@ -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/>. + +Class + Foam::sigStopAtWriteNow + +Description + Signal handler for interupt defined by + OptimisationSwitches::stopAtWriteNowSignal + + Write and stop the job. + +SourceFiles + sigStopAtWriteNow.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sigStopAtWriteNow_H +#define sigStopAtWriteNow_H + +#include <signal.h> + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class Time; + +/*---------------------------------------------------------------------------*\ + Class sigStopAtWriteNow Declaration +\*---------------------------------------------------------------------------*/ + +class sigStopAtWriteNow +{ + // Private data + + //- number of signal to use + static int signal_; + + //- Saved old signal trapping setting + static struct sigaction oldAction_; + + // Private Member Functions + + static void sigHandler(int); + + +public: + + // Constructors + + //- Construct null + sigStopAtWriteNow(); + + //- Construct from components + sigStopAtWriteNow(const bool verbose, const Time& runTime); + + + //- Destructor + ~sigStopAtWriteNow(); + + + // Member functions + + //- Is active? + bool active() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OSspecific/POSIX/signals/sigWriteNow.C b/src/OSspecific/POSIX/signals/sigWriteNow.C new file mode 100644 index 0000000000000000000000000000000000000000..6ad98a6190d88bac2a87ad29edf3416da049d917 --- /dev/null +++ b/src/OSspecific/POSIX/signals/sigWriteNow.C @@ -0,0 +1,122 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "sigWriteNow.H" +#include "error.H" +#include "JobInfo.H" +#include "IOstreams.H" +#include "Time.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +// Signal number to catch +int Foam::sigWriteNow::signal_ +( + debug::optimisationSwitch("writeNowSignal", -1) +); + +static Foam::Time* runTimePtr_ = NULL; + + +struct sigaction Foam::sigWriteNow::oldAction_; + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::sigWriteNow::sigHandler(int) +{ + Info<< "sigWriteNow :" + << " setting up write at end of the next iteration" << nl << endl; + runTimePtr_->writeOnce(); + + //// Throw signal (to old handler) + //raise(signal_); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sigWriteNow::sigWriteNow() +{} + + +Foam::sigWriteNow::sigWriteNow(const bool verbose, Time& runTime) +{ + if (signal_ >= 0) + { + // Store runTime + runTimePtr_ = &runTime; + + struct sigaction newAction; + newAction.sa_handler = sigHandler; + newAction.sa_flags = SA_NODEFER; + sigemptyset(&newAction.sa_mask); + if (sigaction(signal_, &newAction, &oldAction_) < 0) + { + FatalErrorIn + ( + "Foam::sigWriteNow::sigWriteNow(const bool, const Time&)" + ) << "Cannot set " << signal_ << " trapping" + << abort(FatalError); + } + + if (verbose) + { + Info<< "sigWriteNow :" + << " Enabling writing upon signal " << signal_ + << endl; + } + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sigWriteNow::~sigWriteNow() +{ + // Reset old handling + if (signal_ > 0) + { + if (sigaction(signal_, &oldAction_, NULL) < 0) + { + FatalErrorIn + ( + "Foam::sigWriteNow::~sigWriteNow()" + ) << "Cannot reset " << signal_ << " trapping" + << abort(FatalError); + } + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::sigWriteNow::active() const +{ + return signal_ > 0; +} + + +// ************************************************************************* // diff --git a/src/OSspecific/POSIX/signals/sigWriteNow.H b/src/OSspecific/POSIX/signals/sigWriteNow.H new file mode 100644 index 0000000000000000000000000000000000000000..477bae825f4b2a48ffd06d3b9203b7af57f40896 --- /dev/null +++ b/src/OSspecific/POSIX/signals/sigWriteNow.H @@ -0,0 +1,101 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::sigWriteNow + +Description + Signal handler for interupt defined by OptimisationSwitches::writeNowSignal + + Write once and continue. + +SourceFiles + sigWriteNow.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sigWriteNow_H +#define sigWriteNow_H + +#include <signal.h> + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class Time; + +/*---------------------------------------------------------------------------*\ + Class sigWriteNow Declaration +\*---------------------------------------------------------------------------*/ + +class sigWriteNow +{ + // Private data + + //- number of signal to use + static int signal_; + + //- Saved old signal trapping setting + static struct sigaction oldAction_; + + // Private Member Functions + + static void sigHandler(int); + + +public: + + friend class sigStopAtWriteNow; + + // Constructors + + //- Construct null + sigWriteNow(); + + //- Construct from components + sigWriteNow(const bool verbose, Time& runTime); + + + //- Destructor + ~sigWriteNow(); + + + // Member functions + + //- Is active? + bool active() const; + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/db/Time/Time.C b/src/OpenFOAM/db/Time/Time.C index 727f21c98ceb1217d2eb3751e55c076049651d04..0051d23153104f301817f973b09f813bc0de8932 100644 --- a/src/OpenFOAM/db/Time/Time.C +++ b/src/OpenFOAM/db/Time/Time.C @@ -81,10 +81,17 @@ void Foam::Time::adjustDeltaT() { if (writeControl_ == wcAdjustableRunTime) { + scalar interval = writeInterval_; + if (secondaryWriteControl_ == wcAdjustableRunTime) + { + interval = min(interval, secondaryWriteInterval_); + } + + scalar timeToNextWrite = max ( 0.0, - (outputTimeIndex_ + 1)*writeInterval_ - (value() - startTime_) + (outputTimeIndex_ + 1)*interval - (value() - startTime_) ); scalar nSteps = timeToNextWrite/deltaT_ - SMALL; @@ -252,8 +259,13 @@ Foam::Time::Time stopAt_(saEndTime), writeControl_(wcTimeStep), writeInterval_(GREAT), + secondaryWriteControl_(wcTimeStep), + secondaryWriteInterval_(labelMax), purgeWrite_(0), + writeOnce_(false), subCycling_(false), + sigWriteNow_(true, *this), + sigStopAtWriteNow_(true, *this), writeFormat_(IOstream::ASCII), writeVersion_(IOstream::currentVersion), @@ -339,8 +351,13 @@ Foam::Time::Time stopAt_(saEndTime), writeControl_(wcTimeStep), writeInterval_(GREAT), + secondaryWriteControl_(wcTimeStep), + secondaryWriteInterval_(labelMax), purgeWrite_(0), + writeOnce_(false), subCycling_(false), + sigWriteNow_(true, *this), + sigStopAtWriteNow_(true, *this), writeFormat_(IOstream::ASCII), writeVersion_(IOstream::currentVersion), @@ -429,8 +446,13 @@ Foam::Time::Time stopAt_(saEndTime), writeControl_(wcTimeStep), writeInterval_(GREAT), + secondaryWriteControl_(wcTimeStep), + secondaryWriteInterval_(labelMax), purgeWrite_(0), + writeOnce_(false), subCycling_(false), + sigWriteNow_(true, *this), + sigStopAtWriteNow_(true, *this), writeFormat_(IOstream::ASCII), writeVersion_(IOstream::currentVersion), @@ -521,7 +543,10 @@ Foam::Time::Time stopAt_(saEndTime), writeControl_(wcTimeStep), writeInterval_(GREAT), + secondaryWriteControl_(wcTimeStep), + secondaryWriteInterval_(labelMax), purgeWrite_(0), + writeOnce_(false), subCycling_(false), writeFormat_(IOstream::ASCII), @@ -944,6 +969,35 @@ Foam::Time& Foam::Time::operator++() if (!subCycling_) { + if (sigStopAtWriteNow_.active() || sigWriteNow_.active()) + { + // A signal might have been sent on one processor only + // Reduce so all decide the same. + + label flag = 0; + if (sigStopAtWriteNow_.active() && stopAt_ == saWriteNow) + { + flag += 1; + } + if (sigWriteNow_.active() && writeOnce_) + { + flag += 2; + } + reduce(flag, maxOp<label>()); + + if (flag & 1) + { + stopAt_ = saWriteNow; + } + if (flag & 2) + { + writeOnce_ = true; + } + } + + + outputTime_ = false; + switch (writeControl_) { case wcTimeStep: @@ -964,10 +1018,6 @@ Foam::Time& Foam::Time::operator++() outputTime_ = true; outputTimeIndex_ = outputIndex; } - else - { - outputTime_ = false; - } } break; @@ -983,10 +1033,6 @@ Foam::Time& Foam::Time::operator++() outputTime_ = true; outputTimeIndex_ = outputIndex; } - else - { - outputTime_ = false; - } } break; @@ -1002,14 +1048,69 @@ Foam::Time& Foam::Time::operator++() outputTime_ = true; outputTimeIndex_ = outputIndex; } - else + } + break; + } + + + // Adapt for secondaryWrite controls + switch (secondaryWriteControl_) + { + case wcTimeStep: + outputTime_ = + outputTime_ + || !(timeIndex_ % label(secondaryWriteInterval_)); + break; + + case wcRunTime: + case wcAdjustableRunTime: + { + label outputIndex = label + ( + ((value() - startTime_) + 0.5*deltaT_) + / secondaryWriteInterval_ + ); + + if (outputIndex > outputTimeIndex_) + { + outputTime_ = true; + outputTimeIndex_ = outputIndex; + } + } + break; + + case wcCpuTime: + { + label outputIndex = label + ( + returnReduce(elapsedCpuTime(), maxOp<double>()) + / secondaryWriteInterval_ + ); + if (outputIndex > outputTimeIndex_) { - outputTime_ = false; + outputTime_ = true; + outputTimeIndex_ = outputIndex; + } + } + break; + + case wcClockTime: + { + label outputIndex = label + ( + returnReduce(label(elapsedClockTime()), maxOp<label>()) + / secondaryWriteInterval_ + ); + if (outputIndex > outputTimeIndex_) + { + outputTime_ = true; + outputTimeIndex_ = outputIndex; } } break; } + // see if endTime needs adjustment to stop at the next run()/end() check if (!end()) { @@ -1027,6 +1128,14 @@ Foam::Time& Foam::Time::operator++() endTime_ = value(); } } + + // Override outputTime if one-shot writing + if (writeOnce_) + { + outputTime_ = true; + writeOnce_ = false; + } + } return *this; diff --git a/src/OpenFOAM/db/Time/Time.H b/src/OpenFOAM/db/Time/Time.H index a84b558a38ff04ef59e141fae2571ffb8155ffaa..0cfe4c51fa9888c5bbebfecf8d271d939207fee8 100644 --- a/src/OpenFOAM/db/Time/Time.H +++ b/src/OpenFOAM/db/Time/Time.H @@ -52,6 +52,8 @@ SourceFiles #include "dlLibraryTable.H" #include "functionObjectList.H" #include "fileMonitor.H" +#include "sigWriteNow.H" +#include "sigStopAtWriteNow.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -130,15 +132,35 @@ protected: scalar writeInterval_; + // Additional writing + + writeControls secondaryWriteControl_; + + scalar secondaryWriteInterval_; + + label purgeWrite_; mutable FIFOStack<word> previousOutputTimes_; + // One-shot writing + bool writeOnce_; + //- Is the time currently being sub-cycled? bool subCycling_; //- If time is being sub-cycled this is the previous TimeState autoPtr<TimeState> prevTimeState_; + + // Signal handlers for secondary writing + + //- Enable one-shot writing upon signal + sigWriteNow sigWriteNow_; + + //- Enable write and clean exit upon signal + sigStopAtWriteNow sigStopAtWriteNow_; + + //- Time directory name format static fmtflags format_; @@ -357,12 +379,16 @@ public: IOstream::compressionType ) const; - //- Write the objects now and continue the run + //- Write the objects now (not at end of iteration) and continue + // the run bool writeNow(); - //- Write the objects now and end the run + //- Write the objects now (not at end of iteration) and end the run bool writeAndEnd(); + //- Write the objects once (one shot) and continue the run + void writeOnce(); + // Access diff --git a/src/OpenFOAM/db/Time/TimeIO.C b/src/OpenFOAM/db/Time/TimeIO.C index 905e508a6aa4d989d5e6764ba698e4ff43d13113..9ed963127d8f7eba9b3dbf92189cd64b2b7075a9 100644 --- a/src/OpenFOAM/db/Time/TimeIO.C +++ b/src/OpenFOAM/db/Time/TimeIO.C @@ -58,6 +58,45 @@ void Foam::Time::readDict() controlDict_.lookup("writeFrequency") >> writeInterval_; } + + // Additional writing + if (controlDict_.found("secondaryWriteControl")) + { + secondaryWriteControl_ = writeControlNames_.read + ( + controlDict_.lookup("secondaryWriteControl") + ); + + if + ( + controlDict_.readIfPresent + ( + "secondaryWriteInterval", + secondaryWriteInterval_ + ) + ) + { + if + ( + secondaryWriteControl_ + == wcTimeStep && label(secondaryWriteInterval_) < 1 + ) + { + FatalIOErrorIn("Time::readDict()", controlDict_) + << "secondaryWriteInterval < 1" + << " for secondaryWriteControl timeStep" + << exit(FatalIOError); + } + } + else + { + controlDict_.lookup("secondaryWriteFrequency") + >> secondaryWriteInterval_; + } + } + + + if (oldWriteInterval != writeInterval_) { switch (writeControl_) @@ -310,4 +349,10 @@ bool Foam::Time::writeAndEnd() } +void Foam::Time::writeOnce() +{ + writeOnce_ = true; +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/db/objectRegistry/objectRegistry.C b/src/OpenFOAM/db/objectRegistry/objectRegistry.C index f140045d06d995684b431637aac3de6cab3a5d9a..6fcaed2e2eab3df27a80a878edccd0d3e469ed0b 100644 --- a/src/OpenFOAM/db/objectRegistry/objectRegistry.C +++ b/src/OpenFOAM/db/objectRegistry/objectRegistry.C @@ -164,10 +164,13 @@ Foam::label Foam::objectRegistry::getEvent() const if (event_ == labelMax) { - WarningIn("objectRegistry::getEvent() const") - << "Event counter has overflowed. " - << "Resetting counter on all dependent objects." << nl - << "This might cause extra evaluations." << endl; + if (objectRegistry::debug) + { + WarningIn("objectRegistry::getEvent() const") + << "Event counter has overflowed. " + << "Resetting counter on all dependent objects." << nl + << "This might cause extra evaluations." << endl; + } // Reset event counter curEvent = 1; diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C index 75bc1cb4dd82c24e77fef88fc625f3d908af6863..b2ba573dde9ca0fb0e5851e483136b98c7c156b9 100644 --- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C +++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C @@ -124,11 +124,7 @@ DimensionedField<Type, GeoMesh>::DimensionedField const DimensionedField<Type, GeoMesh>& df ) : -# ifdef ConstructFromTmp regIOobject(df), -# else - regIOobject(df, true), -# endif Field<Type>(df), mesh_(df.mesh_), dimensions_(df.dimensions_) @@ -142,7 +138,7 @@ DimensionedField<Type, GeoMesh>::DimensionedField bool reUse ) : - regIOobject(df, true), + regIOobject(df, reUse), Field<Type>(df, reUse), mesh_(df.mesh_), dimensions_(df.dimensions_) @@ -169,7 +165,7 @@ DimensionedField<Type, GeoMesh>::DimensionedField const tmp<DimensionedField<Type, GeoMesh> >& tdf ) : - regIOobject(tdf(), true), + regIOobject(tdf(), tdf.isTmp()), Field<Type> ( const_cast<DimensionedField<Type, GeoMesh>&>(tdf()), diff --git a/src/OpenFOAM/global/argList/argList.C b/src/OpenFOAM/global/argList/argList.C index 3c42659bf20572a86372c9c7ec0f9162bf9f4129..5d8758e2d2978145f83e5e15935c7b254ab7e582 100644 --- a/src/OpenFOAM/global/argList/argList.C +++ b/src/OpenFOAM/global/argList/argList.C @@ -726,7 +726,7 @@ Foam::argList::argList } - wordList slaveProcs; + stringList slaveProcs; // collect slave machine/pid if (parRunControl_.parRun()) @@ -734,7 +734,7 @@ Foam::argList::argList if (Pstream::master()) { slaveProcs.setSize(Pstream::nProcs() - 1); - word slaveMachine; + string slaveMachine; label slavePid; label procI = 0; diff --git a/src/OpenFOAM/include/OSspecific.H b/src/OpenFOAM/include/OSspecific.H index dc12ef753f8e1b94ae24f429c5893b4136a90a73..d39a1c54cbb2e4b931cdf768215525a23ef14e3a 100644 --- a/src/OpenFOAM/include/OSspecific.H +++ b/src/OpenFOAM/include/OSspecific.H @@ -69,13 +69,13 @@ bool setEnv(const word& name, const std::string& value, const bool overwrite); //- Return the system's host name, as per hostname(1) // Optionally with the full name (as per the '-f' option) -word hostName(const bool full=false); +string hostName(const bool full=false); //- Return the system's domain name, as per hostname(1) with the '-d' option -word domainName(); +string domainName(); //- Return the user's login name -word userName(); +string userName(); //- Is user administrator bool isAdministrator(); @@ -84,7 +84,7 @@ bool isAdministrator(); fileName home(); //- Return home directory path name for a particular user -fileName home(const word& userName); +fileName home(const string& userName); //- Return current working directory path name fileName cwd(); @@ -189,10 +189,10 @@ unsigned int sleep(const unsigned int); void fdClose(const int); //- Check if machine is up by pinging given port -bool ping(const word&, const label port, const label timeOut); +bool ping(const string&, const label port, const label timeOut); //- Check if machine is up by pinging port 22 (ssh) and 222 (rsh) -bool ping(const word&, const label timeOut=10); +bool ping(const string&, const label timeOut=10); //- Execute the specified command int system(const std::string& command); diff --git a/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C b/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C index 69f07a6372a42d631591380e45316baa0b3a1977..5e8909549815b851539bb97315976c77649158ea 100644 --- a/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C +++ b/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C @@ -88,7 +88,7 @@ Foam::interpolationTable<Type>::interpolationTable(const fileName& fName) List<Tuple2<scalar, Type> >(), boundsHandling_(interpolationTable::WARN), fileName_(fName), - reader_(new openFoamTableReader<Type>()) + reader_(new openFoamTableReader<Type>(dictionary())) { readTable(); } 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/OpenFOAM/primitives/strings/stringOps/stringOps.C b/src/OpenFOAM/primitives/strings/stringOps/stringOps.C index 32e88b16f80b03d8bb9b85796679d617efc7a226..e7493b78e954694fa8b4f5e956b51a86f512133b 100644 --- a/src/OpenFOAM/primitives/strings/stringOps/stringOps.C +++ b/src/OpenFOAM/primitives/strings/stringOps/stringOps.C @@ -548,7 +548,7 @@ Foam::string& Foam::stringOps::inplaceExpand // ~OpenFOAM => site/user OpenFOAM configuration directory // ~user => home directory for specified user - word user; + string user; fileName file; if ((begVar = s.find('/')) != string::npos) diff --git a/src/Pstream/dummy/UIPread.C b/src/Pstream/dummy/UIPread.C index 1dc91fe038d3ccdb6104b706bf983b426d714092..fc3b4466e2a32bad3e4e4c1c4b74e39f78bbf488 100644 --- a/src/Pstream/dummy/UIPread.C +++ b/src/Pstream/dummy/UIPread.C @@ -53,15 +53,42 @@ Foam::UIPstream::UIPstream { notImplemented ( - "UIPstream::UIPstream" - "(" - "const commsTypes," - "const int fromProcNo," - "DynamicList<char>&," - "label&," - "const int tag," - "const bool," - "streamFormat, versionNumber" + "UIPstream::UIPstream\n" + "(\n" + "const commsTypes,\n" + "const int,\n" + "DynamicList<char>&,\n" + "label&,\n" + "const int,\n" + "const bool,\n" + "streamFormat,\n" + "versionNumber\n" + ")" + ); +} + + +Foam::UIPstream::UIPstream +( + const int fromProcNo, + PstreamBuffers& buffers +) +: + UPstream(buffers.commsType_), + Istream(buffers.format_, buffers.version_), + fromProcNo_(fromProcNo), + externalBuf_(buffers.recvBuf_[fromProcNo]), + externalBufPosition_(buffers.recvBufPos_[fromProcNo]), + tag_(buffers.tag_), + clearAtEnd_(true), + messageSize_(0) +{ + notImplemented + ( + "UIPstream::UIPstream\n" + "(\n" + "const int,\n" + "PstreamBuffers&\n" ")" ); } diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 971740c7dab76d78af367655372dfa0d47096db5..3bcd7b80d55dfec16641962a6dbf6ff0e12c12fd 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -334,6 +334,7 @@ $(snGradSchemes)/snGradScheme/snGradSchemes.C $(snGradSchemes)/correctedSnGrad/correctedSnGrads.C $(snGradSchemes)/limitedSnGrad/limitedSnGrads.C $(snGradSchemes)/uncorrectedSnGrad/uncorrectedSnGrads.C +$(snGradSchemes)/orthogonalSnGrad/orthogonalSnGrads.C /* $(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGradData.C $(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C @@ -388,6 +389,7 @@ $(basicSource)/basicSource/IObasicSourceList.C $(basicSource)/actuationDiskSource/actuationDiskSource.C $(basicSource)/radialActuationDiskSource/radialActuationDiskSource.C $(basicSource)/explicitSource/explicitSource.C +$(basicSource)/explicitSetValue/explicitSetValue.C LIB = $(FOAM_LIBBIN)/libfiniteVolume diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C index a5c8208d7b21d0fab774d96fbd378bcbda62f335..1b1436767e9bc85a130a485a8997064a0142ba8c 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C @@ -48,14 +48,14 @@ void Foam::actuationDiskSource::addActuationDiskAxialInertialResistance E.xx() = uniDiskDir.x(); E.yy() = uniDiskDir.y(); E.zz() = uniDiskDir.z(); - const vectorField U1((1.0 - a)*U); + forAll(cells, i) { - T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U1[cells[i]])*a/(1.0 - a); + T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U[cells[i]])*a/(1.0 - a); } forAll(cells, i) { - Usource[cells[i]] += ((Vcells[cells[i]]/V())*T[i]*E) & U1[cells[i]]; + Usource[cells[i]] += ((Vcells[cells[i]]/V())*T[i]*E) & U[cells[i]]; } } diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C index 6bca9b3e3095cf3c8175e8e8514e0f7cf883f116..114e6e58cbe692549e742f4d515bf2cbbe962c74 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C @@ -116,10 +116,15 @@ void Foam::basicSource::setCellSet() label globalCellI = returnReduce(cellI, maxOp<label>()); if (globalCellI < 0) { - WarningIn("TimeActivatedExplicitSource<Type>::setCellIds()") - << "Unable to find owner cell for point " << points_[i] - << endl; + WarningIn + ( + "TimeActivatedExplicitSource<Type>::setCellIds()" + ) + << "Unable to find owner cell for point " << points_[i] + << endl; + } + } cells_ = selectedCells.toc(); @@ -270,4 +275,30 @@ bool Foam::basicSource::isActive() } +void 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)" + ); +} + + +void 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 0ab5002fade4328ae373c145a169efd60c90c501..43c9e8bd47028e8fa32b05774ffdb1946d0be1f7 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.H +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.H @@ -56,13 +56,13 @@ Description selectionMode points; cellSet c0; + points // list of points when selectionMode = points + ( + (-0.088 0.007 -0.02) + (-0.028 0.007 -0.02) + ); explicitSourceCoeffs { - points // list of points when selectionMode = points - ( - (-0.088 0.007 -0.02) - (-0.028 0.007 -0.02) - ); volumeMode specific; //absolute fieldData //field data { @@ -323,10 +323,14 @@ public: // Evaluation //- Add source term to vector fvMatrix - virtual void addSu(fvMatrix<vector>& Eqn) = 0; + virtual void addSu(fvMatrix<vector>& Eqn); //- Add source term to scalar fvMatrix - virtual void addSu(fvMatrix<scalar>& Eqn) = 0; + virtual void addSu(fvMatrix<scalar>& Eqn); + + //- Set constant value on field + virtual void setValue(fvMatrix<scalar>& Eq); + // I-O diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C index baefe473f8eac1abdd065814e6ed12c1916d1224..e9827409c5aef81c4de1449e5ae03aa1cddbb58e 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C @@ -90,6 +90,19 @@ void Foam::basicSourceList::addSu(fvMatrix<vector>& Eqn) } +void Foam::basicSourceList::setValue(fvMatrix<scalar>& Eqn) +{ + + forAll(*this, i) + { + if (this->operator[](i).isActive()) + { + this->operator[](i).setValue(Eqn); + } + } +} + + bool Foam::basicSourceList::read(const dictionary& dict) { bool allOk = true; diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.H index 2acf4b8e96b1066fc90763f4ba995947adf3ad7e..30e6501de129f4724e3fe993f3b96231ef6b26a1 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.H +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.H @@ -92,6 +92,9 @@ public: //- Add source terms to vector fvMatrix void addSu(fvMatrix<vector>& Eq); + //- Set constant value on field + void setValue(fvMatrix<scalar>& Eq); + // I-O diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.C new file mode 100644 index 0000000000000000000000000000000000000000..2a14b793a8de3db75521adcff5d0a1ec05704b6b --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.C @@ -0,0 +1,117 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "explicitSetValue.H" +#include "fvMesh.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" +#include "HashSet.H" + +// * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(explicitSetValue, 0); + addToRunTimeSelectionTable + ( + basicSource, + explicitSetValue, + dictionary + ); +} + + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +void Foam::explicitSetValue::setFieldData(const dictionary& dict) +{ + scalarFields_.clear(); + vectorFields_.clear(); + + wordList fieldTypes(dict.toc().size()); + wordList fieldNames(dict.toc().size()); + + forAll(dict.toc(), i) + { + const word& fieldName = dict.toc()[i]; + IOobject io + ( + fieldName, + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ); + if (io.headerOk()) + { + fieldTypes[i] = io.headerClassName(); + fieldNames[i] = dict.toc()[i]; + } + else + { + FatalErrorIn + ( + "explicitSetValue::setFieldData" + ) << "header not OK " << io.name() + << exit(FatalError); + } + } + + addField(scalarFields_, fieldTypes, fieldNames, dict); + addField(vectorFields_, fieldTypes, fieldNames, dict); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::explicitSetValue::explicitSetValue +( + const word& name, + const word& modelType, + const dictionary& dict, + const fvMesh& mesh +) +: + basicSource(name, modelType, dict, mesh), + dict_(dict.subDict(modelType + "Coeffs")) +{ + setFieldData(dict_.subDict("fieldData")); +} + + +void Foam::explicitSetValue::setValue(fvMatrix<scalar>& Eqn) +{ + setFieldValue(Eqn, scalarFields_[Eqn.psi().name()]); +} + + +void Foam::explicitSetValue::setValue(fvMatrix<vector>& Eqn) +{ + setFieldValue(Eqn, vectorFields_[Eqn.psi().name()]); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.H new file mode 100644 index 0000000000000000000000000000000000000000..643e6b618243d40c4bb7401e13c7bdcc4e7a1352 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.H @@ -0,0 +1,184 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::explicitSetValue + +Description + Explicit set values on fields. + + Sources described by: + + explicitSetValueCoeffs + { + fieldData // field data - usage for multiple fields + { + k 30.7; + epsilon 1.5; + } + } + +SourceFiles + explicitSetValue.C + +\*---------------------------------------------------------------------------*/ + +#ifndef explicitSetValue_H +#define explicitSetValue_H + +#include "cellSet.H" +#include "volFieldsFwd.H" +#include "DimensionedField.H" +#include "basicSource.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class explicitSetValue Declaration +\*---------------------------------------------------------------------------*/ + +class explicitSetValue +: + public basicSource +{ + // Private data + + //- List of field types + HashTable<scalar> scalarFields_; + HashTable<vector> vectorFields_; + + //- Set value to field + template<class Type> + void setFieldValue(fvMatrix<Type>&, const Type&) const; + + //- Add field names and values to field table for types. + template<class Type> + void addField + ( + HashTable<Type>& fields, + const wordList& fieldTypes, + const wordList& fieldNames, + const dictionary& dict + ); + + +protected: + + // Protected data + + //- Sub dictionary for time activated explicit sources + const dictionary& dict_; + + + // Protected functions + + //- Set the local field data + void setFieldData(const dictionary& dict); + + +public: + + //- Runtime type information + TypeName("explicitSetValue"); + + + // Constructors + + //- Construct from components + explicitSetValue + ( + const word& name, + const word& modelType, + const dictionary& dict, + const fvMesh& mesh + ); + + //- Return clone + autoPtr<explicitSetValue> clone() const + { + notImplemented + ( + "autoPtr<explicitSetValue> clone() const" + ); + return autoPtr<explicitSetValue>(NULL); + } + + + // Member Functions + + + // Edit + + //- Return points + inline const List<point>& points() const; + + + // Evaluation + + //- Set value on vector field + virtual void setValue(fvMatrix<vector>& UEqn); + + //- Set value on scalar field + virtual void setValue(fvMatrix<scalar>& UEqn); + + + // I-O + + //- Write the source properties + virtual void writeData(Ostream&) const; + + //- Read fieldData in sub-dictionary + virtual bool read(const dictionary& dict); + + //- Ostream operator + friend Ostream& operator<< + ( + Ostream& os, + const explicitSetValue& source + ); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "explicitSetValueIO.C" +#include "explicitSetValueI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "explicitSetValueTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueI.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueI.H new file mode 100644 index 0000000000000000000000000000000000000000..35aa075140ce370c0dcd761729793081d8189d21 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueI.H @@ -0,0 +1,37 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "explicitSetValue.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline const Foam::List<Foam::point>& +Foam::explicitSetValue::points() const +{ + return points_; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueIO.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueIO.C new file mode 100644 index 0000000000000000000000000000000000000000..dce0a53ddf69a4776aba28b21bc040ab76824773 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueIO.C @@ -0,0 +1,79 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "explicitSetValue.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +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; +} + + +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")); + return true; + } + else + { + return false; + } +} + + +// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // + +Foam::Ostream& Foam::operator<<(Ostream& os, const explicitSetValue& source) +{ + source.writeData(os); + return os; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueTemplates.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..f67c3183436edcb72c1c5171f3fe0fb86a93800a --- /dev/null +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueTemplates.C @@ -0,0 +1,105 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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/>. + +\*---------------------------------------------------------------------------*/ + +template <class Type> +void Foam::explicitSetValue::setFieldValue +( + fvMatrix<Type>& Eqn, + const Type& value +) const +{ + Type data = value; + + DimensionedField<Type, volMesh> rhs + ( + IOobject + ( + "rhs", + Eqn.psi().mesh().time().timeName(), + Eqn.psi().mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + Eqn.psi().mesh(), + dimensioned<Type> + ( + "zero", + dimless, + pTraits<Type>::zero + ) + ); + + List<Type> values(this->cells().size()); + + forAll (values, i) + { + values[i] = data; + } + + Eqn.setValues(this->cells(), values); +} + + +template <class Type> +void Foam::explicitSetValue::addField +( + HashTable<Type>& fields, + const wordList& fieldTypes, + const wordList& fieldNames, + const dictionary& fieldDataDict +) +{ + typedef GeometricField<Type, fvPatchField, volMesh> geometricField; + + forAll (fieldTypes, fieldI) + { + word fieldName = fieldNames[fieldI]; + word fieldType = fieldTypes[fieldI]; + + if + ( + ( + fieldType + == GeometricField<Type, fvPatchField, volMesh>::typeName + ) && + ( + this->mesh().foundObject<geometricField>(fieldName) + ) + ) + { + Type fieldValue = fieldDataDict.lookupOrDefault<Type> + ( + fieldName, + pTraits<Type>::zero + ); + + fields.insert(fieldName, fieldValue); + } + } +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.C index 3299f1c5566b2b672c6a2fa3a1c46ae8fbd23238..de01ebb5309c7e2bd4975a1ac4013ae3a6336350 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.C @@ -75,7 +75,7 @@ void Foam::explicitSource::setFieldData(const dictionary& dict) IOobject io ( fieldName, - this->mesh().time().timeName(0), + this->mesh().time().timeName(), this->mesh(), IOobject::NO_READ, IOobject::NO_WRITE, diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C index 6c8469fda8958d4276e14a586cb83974c02512b1..47ac4d0f4083862825de8baf0832ac88e8620ba5 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C @@ -53,7 +53,6 @@ addRadialActuationDiskAxialInertialResistance E.xx() = uniDiskDir.x(); E.yy() = uniDiskDir.y(); E.zz() = uniDiskDir.z(); - const vectorField U1((1.0 - a)*U); const Field<vector> zoneCellCentres(mesh().cellCentres(), cells); const Field<scalar> zoneCellVolumes(mesh().cellVolumes(), cells); @@ -68,7 +67,7 @@ addRadialActuationDiskAxialInertialResistance forAll(cells, i) { - T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U1[cells[i]])*a/(1.0 - a); + T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U[cells[i]])*a/(1.0 - a); scalar r = mag(mesh().cellCentres()[cells[i]] - avgCentre); @@ -79,7 +78,7 @@ addRadialActuationDiskAxialInertialResistance forAll(cells, i) { - Usource[cells[i]] += ((Vcells[cells[i]]/V())*Tr[i]*E) & U1[cells[i]]; + Usource[cells[i]] += ((Vcells[cells[i]]/V())*Tr[i]*E) & U[cells[i]]; } if (debug) diff --git a/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C index eab4f831462022ce305cbad2586be693ae76864f..4bf1a4dd6b50be8dcc9b24fd2910f1ea749c0aca 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C @@ -110,38 +110,34 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const const labelUList& owner = mesh_.owner(); const labelUList& neighbour = mesh_.neighbour(); - // Build the d-vectors - surfaceVectorField d - ( - mesh_.Sf()/(mesh_.magSf()*mesh_.deltaCoeffs()) - ); - - if (!mesh_.orthogonal()) - { - d -= mesh_.correctionVectors()/mesh_.deltaCoeffs(); - } - + const volVectorField& C = mesh.C(); // Set up temporary storage for the dd tensor (before inversion) symmTensorField dd(mesh_.nCells(), symmTensor::zero); - forAll(owner, faceI) + forAll(owner, facei) { - const symmTensor wdd(1.0/magSqr(d[faceI])*sqr(d[faceI])); + label own = owner[facei]; + label nei = neighbour[facei]; + + vector d = C[nei] - C[own]; - dd[owner[faceI]] += wdd; - dd[neighbour[faceI]] += wdd; + const symmTensor wdd(1.0/magSqr(d[facei])*sqr(d[facei])); + + dd[own] += wdd; + dd[nei] += wdd; } // Visit the boundaries. Coupled boundaries are taken into account // in the construction of d vectors. - forAll(d.boundaryField(), patchI) + forAll(lsP.boundaryField(), patchi) { - const fvsPatchVectorField& pd = d.boundaryField()[patchI]; - - const fvPatch& p = pd.patch(); + const fvPatch& p = lsP.boundaryField()[patchi].patch(); const labelUList& faceCells = p.faceCells(); + // Build the d-vectors + vectorField pd = p.delta(); + forAll(pd, patchFaceI) { dd[faceCells[patchFaceI]] += @@ -232,24 +228,30 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const // Revisit all faces and calculate the lsP and lsN vectors - forAll(owner, faceI) + forAll(owner, facei) { - lsP[faceI] = - (1.0/magSqr(d[faceI]))*(invDd[owner[faceI]] & d[faceI]); + label own = owner[facei]; + label nei = neighbour[facei]; + + vector d = C[nei] - C[own]; - lsN[faceI] = - ((-1.0)/magSqr(d[faceI]))*(invDd[neighbour[faceI]] & d[faceI]); + lsP[facei] = + (1.0/magSqr(d[facei]))*(invDd[owner[facei]] & d); + + lsN[facei] = + ((-1.0)/magSqr(d[facei]))*(invDd[neighbour[facei]] & d); } forAll(lsP.boundaryField(), patchI) { - const fvsPatchVectorField& pd = d.boundaryField()[patchI]; - fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchI]; const fvPatch& p = patchLsP.patch(); const labelUList& faceCells = p.faceCells(); + // Build the d-vectors + vectorField pd = p.delta(); + forAll(p, patchFaceI) { patchLsP[patchFaceI] = diff --git a/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C index b8a3b428ddb63f8fc9be3a11c75465170f0ce320..ead1ed0c2b69b72ff892b4269fe6a4cffa7fb8dd 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C @@ -121,24 +121,12 @@ Foam::fv::fourthGrad<Type>::calcGrad const scalarField& lambdap = lambda.boundaryField()[patchi]; - // Build the d-vectors - vectorField pd - ( - mesh.Sf().boundaryField()[patchi] - / ( - mesh.magSf().boundaryField()[patchi] - * mesh.deltaCoeffs().boundaryField()[patchi] - ) - ); + const fvPatch& p = fGrad.boundaryField()[patchi].patch(); - if (!mesh.orthogonal()) - { - pd -= mesh.correctionVectors().boundaryField()[patchi] - /mesh.deltaCoeffs().boundaryField()[patchi]; - } + const labelUList& faceCells = p.faceCells(); - const labelUList& faceCells = - fGrad.boundaryField()[patchi].patch().faceCells(); + // Build the d-vectors + vectorField pd = p.delta(); const Field<GradType> neighbourSecondfGrad ( 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/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C index 293c0e1d64ac556ff2d07d573ab0028989701c13..6002a07e1c1a2f95a39bed80ec4921a3da8a5052 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C @@ -130,20 +130,7 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const const labelUList& faceCells = p.patch().faceCells(); // Build the d-vectors - vectorField pd - ( - mesh.Sf().boundaryField()[patchi] - / ( - mesh.magSf().boundaryField()[patchi] - * mesh.deltaCoeffs().boundaryField()[patchi] - ) - ); - - if (!mesh.orthogonal()) - { - pd -= mesh.correctionVectors().boundaryField()[patchi] - /mesh.deltaCoeffs().boundaryField()[patchi]; - } + vectorField pd = p.delta(); if (p.coupled()) { @@ -196,21 +183,7 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const const labelUList& faceCells = p.faceCells(); // Build the d-vectors - vectorField pd - ( - mesh.Sf().boundaryField()[patchi] - /( - mesh.magSf().boundaryField()[patchi] - *mesh.deltaCoeffs().boundaryField()[patchi] - ) - ); - - if (!mesh.orthogonal()) - { - pd -= mesh.correctionVectors().boundaryField()[patchi] - /mesh.deltaCoeffs().boundaryField()[patchi]; - } - + vectorField pd = p.delta(); if (p.coupled()) { diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.H b/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.H index ca1b8dac432661701359a9a6c35003116bf431e9..529e55bdbe9f49276a20dbc4a5e3f09d73aa029f 100644 --- a/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.H +++ b/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.H @@ -211,7 +211,7 @@ public: // Add the patch constructor functions to the hash tables -#define makeFvLaplacianTypeScheme(SS, Type, GType) \ +#define makeFvLaplacianTypeScheme(SS, GType, Type) \ \ typedef SS<Type, GType> SS##Type##GType; \ defineNamedTemplateTypeNameAndDebug(SS##Type##GType, 0); \ @@ -224,13 +224,20 @@ public: #define makeFvLaplacianScheme(SS) \ \ makeFvLaplacianTypeScheme(SS, scalar, scalar) \ -makeFvLaplacianTypeScheme(SS, scalar, symmTensor) \ -makeFvLaplacianTypeScheme(SS, scalar, tensor) \ -makeFvLaplacianTypeScheme(SS, vector, scalar) \ -makeFvLaplacianTypeScheme(SS, sphericalTensor, scalar) \ makeFvLaplacianTypeScheme(SS, symmTensor, scalar) \ makeFvLaplacianTypeScheme(SS, tensor, scalar) \ - +makeFvLaplacianTypeScheme(SS, scalar, vector) \ +makeFvLaplacianTypeScheme(SS, symmTensor, vector) \ +makeFvLaplacianTypeScheme(SS, tensor, vector) \ +makeFvLaplacianTypeScheme(SS, scalar, sphericalTensor) \ +makeFvLaplacianTypeScheme(SS, symmTensor, sphericalTensor) \ +makeFvLaplacianTypeScheme(SS, tensor, sphericalTensor) \ +makeFvLaplacianTypeScheme(SS, scalar, symmTensor) \ +makeFvLaplacianTypeScheme(SS, symmTensor, symmTensor) \ +makeFvLaplacianTypeScheme(SS, tensor, symmTensor) \ +makeFvLaplacianTypeScheme(SS, scalar, tensor) \ +makeFvLaplacianTypeScheme(SS, symmTensor, tensor) \ +makeFvLaplacianTypeScheme(SS, tensor, tensor) // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C index 83cddf531effbbbe5f4736b85d436fdb2db0faf0..ba7ddcf6fa46e62244908d007cdafa5a01a02a8e 100644 --- a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C +++ b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C @@ -50,7 +50,7 @@ Foam::fv::correctedSnGrad<Type>::fullGradCorrection // construct GeometricField<Type, fvsPatchField, surfaceMesh> tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tssf = - mesh.correctionVectors() + mesh.nonOrthCorrectionVectors() & linear<typename outerProduct<vector, Type>::type>(mesh).interpolate ( gradScheme<Type>::New @@ -88,7 +88,7 @@ Foam::fv::correctedSnGrad<Type>::correction IOobject::NO_WRITE ), mesh, - vf.dimensions()*mesh.deltaCoeffs().dimensions() + vf.dimensions()*mesh.nonOrthDeltaCoeffs().dimensions() ) ); GeometricField<Type, fvsPatchField, surfaceMesh>& ssf = tssf(); @@ -98,7 +98,7 @@ Foam::fv::correctedSnGrad<Type>::correction ssf.replace ( cmpt, - mesh.correctionVectors() + mesh.nonOrthCorrectionVectors() & linear < typename diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.H b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.H index a962855181e08521b5f02205cd14fee0b830c1b4..1361ed9de27fdce7db5c39374d1ea3ee4e821039 100644 --- a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.H +++ b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.H @@ -96,13 +96,13 @@ public: const GeometricField<Type, fvPatchField, volMesh>& ) const { - return this->mesh().deltaCoeffs(); + return this->mesh().nonOrthDeltaCoeffs(); } //- Return true if this scheme uses an explicit correction virtual bool corrected() const { - return !this->mesh().orthogonal(); + return true; } //- Return the explicit correction to the correctedSnGrad diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/limitedSnGrad/limitedSnGrad.H b/src/finiteVolume/finiteVolume/snGradSchemes/limitedSnGrad/limitedSnGrad.H index 0f14823dc28d889602e308d0fc70aee5e3e17853..9dbe0f01143ec668280c78d189c7d3ed133c512e 100644 --- a/src/finiteVolume/finiteVolume/snGradSchemes/limitedSnGrad/limitedSnGrad.H +++ b/src/finiteVolume/finiteVolume/snGradSchemes/limitedSnGrad/limitedSnGrad.H @@ -120,13 +120,13 @@ public: const GeometricField<Type, fvPatchField, volMesh>& ) const { - return this->mesh().deltaCoeffs(); + return this->mesh().nonOrthDeltaCoeffs(); } //- Return true if this scheme uses an explicit correction virtual bool corrected() const { - return !this->mesh().orthogonal(); + return true; } //- Return the explicit correction to the limitedSnGrad diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrad.C new file mode 100644 index 0000000000000000000000000000000000000000..39866b78ab5f1f9d205a78c0530d555c4bff4519 --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrad.C @@ -0,0 +1,76 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + Simple central-difference snGrad scheme without non-orthogonal correction. + +\*---------------------------------------------------------------------------*/ + +#include "orthogonalSnGrad.H" +#include "volFields.H" +#include "surfaceFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template<class Type> +orthogonalSnGrad<Type>::~orthogonalSnGrad() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class Type> +tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > +orthogonalSnGrad<Type>::correction +( + const GeometricField<Type, fvPatchField, volMesh>& +) const +{ + notImplemented + ( + "orthogonalSnGrad<Type>::correction" + "(const GeometricField<Type, fvPatchField, volMesh>&)" + ); + return tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >(NULL); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrad.H b/src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrad.H new file mode 100644 index 0000000000000000000000000000000000000000..5bd76fe5161ab022e929401fb7b36a6ba0a3a51e --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrad.H @@ -0,0 +1,133 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::fv::orthogonalSnGrad + +Description + Simple central-difference snGrad scheme without non-orthogonal correction. + +SourceFiles + orthogonalSnGrad.C + +\*---------------------------------------------------------------------------*/ + +#ifndef orthogonalSnGrad_H +#define orthogonalSnGrad_H + +#include "snGradScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class orthogonalSnGrad Declaration +\*---------------------------------------------------------------------------*/ + +template<class Type> +class orthogonalSnGrad +: + public snGradScheme<Type> +{ + // Private Member Functions + + //- Disallow default bitwise assignment + void operator=(const orthogonalSnGrad&); + + +public: + + //- Runtime type information + TypeName("uncorrected"); + + + // Constructors + + //- Construct from mesh + orthogonalSnGrad(const fvMesh& mesh) + : + snGradScheme<Type>(mesh) + {} + + + //- Construct from mesh and data stream + orthogonalSnGrad(const fvMesh& mesh, Istream&) + : + snGradScheme<Type>(mesh) + {} + + + //- Destructor + virtual ~orthogonalSnGrad(); + + + // Member Functions + + //- Return the interpolation weighting factors for the given field + virtual tmp<surfaceScalarField> deltaCoeffs + ( + const GeometricField<Type, fvPatchField, volMesh>& + ) const + { + return this->mesh().deltaCoeffs(); + } + + //- Return true if this scheme uses an explicit correction + virtual bool corrected() const + { + return false; + } + + //- Return the explicit correction to the orthogonalSnGrad + // for the given field + virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > + correction(const GeometricField<Type, fvPatchField, volMesh>&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "orthogonalSnGrad.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrads.C b/src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrads.C new file mode 100644 index 0000000000000000000000000000000000000000..4597bb43ffeaf88235a495e1e2ce36f363f95b11 --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrads.C @@ -0,0 +1,42 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + Simple central-difference snGrad scheme without non-orthogonal correction. + +\*---------------------------------------------------------------------------*/ + +#include "orthogonalSnGrad.H" +#include "fvMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace fv +{ + makeSnGradScheme(orthogonalSnGrad) +} +} + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H index e5059d638269020d985b74c315051789572b2d17..b4892095771727fdaaf39ea14e718501a9b4f932 100644 --- a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H +++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H @@ -109,7 +109,7 @@ public: const GeometricField<Type, fvPatchField, volMesh>& ) const { - return this->mesh().deltaCoeffs(); + return this->mesh().nonOrthDeltaCoeffs(); } //- Return true if this scheme uses an explicit correction diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C index b0462cc047976d92d6c09a0ef743bc22d25751fb..fe0858e1c6a6e9d4557e2b5ffdc85ad3336615d9 100644 --- a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C +++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C @@ -247,7 +247,7 @@ Foam::label Foam::quadraticFitSnGradData::calcFit scalarList singVals(minSize_); label nSVDzeros = 0; - const scalar deltaCoeff = mesh().deltaCoeffs()[faci]; + const scalar deltaCoeff = deltaCoeffs()[faci]; bool goodFit = false; for (int iIt = 0; iIt < 10 && !goodFit; iIt++) diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/uncorrectedSnGrad/uncorrectedSnGrad.H b/src/finiteVolume/finiteVolume/snGradSchemes/uncorrectedSnGrad/uncorrectedSnGrad.H index e52db735f16b9f9a3eda66f6540439b113c6872c..082eeabbd0af37c9f5295878037503314f419b66 100644 --- a/src/finiteVolume/finiteVolume/snGradSchemes/uncorrectedSnGrad/uncorrectedSnGrad.H +++ b/src/finiteVolume/finiteVolume/snGradSchemes/uncorrectedSnGrad/uncorrectedSnGrad.H @@ -96,7 +96,7 @@ public: const GeometricField<Type, fvPatchField, volMesh>& ) const { - return this->mesh().deltaCoeffs(); + return this->mesh().nonOrthDeltaCoeffs(); } //- Return true if this scheme uses an explicit correction diff --git a/src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.H index 255b20043c5e035efefa23d32287c71b196a7fd7..7fc59ffb9502b3e2785d9281f2361db190a7e4bc 100644 --- a/src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.H @@ -66,9 +66,6 @@ protected: //- Make patch weighting factors virtual void makeWeights(scalarField&) const = 0; - //- Make patch face - neighbour cell distances - virtual void makeDeltaCoeffs(scalarField&) const = 0; - public: diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C index c453c91910fc20946044580716dae8afbcb0b253..d127cf3a022750d308ebe59512c5045d20ecfbdb 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C @@ -57,25 +57,6 @@ void Foam::cyclicFvPatch::makeWeights(scalarField& w) const } -// Make patch face - neighbour cell distances -void Foam::cyclicFvPatch::makeDeltaCoeffs(scalarField& dc) const -{ - //const cyclicPolyPatch& nbrPatch = cyclicPolyPatch_.neighbPatch(); - const cyclicFvPatch& nbrPatch = neighbFvPatch(); - - const scalarField deltas(nf() & fvPatch::delta()); - const scalarField nbrDeltas(nbrPatch.nf() & nbrPatch.fvPatch::delta()); - - forAll(deltas, facei) - { - scalar di = deltas[facei]; - scalar dni = nbrDeltas[facei]; - - dc[facei] = 1.0/(di + dni); - } -} - - // Return delta (P to N) vectors across coupled patch Foam::tmp<Foam::vectorField> Foam::cyclicFvPatch::delta() const { diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.H index 107a085927dbe51569f32436708e48865da77924..fe0decf711a262d60b795f0ad03835fd44b7aa54 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.H @@ -66,9 +66,6 @@ protected: //- Make patch weighting factors void makeWeights(scalarField&) const; - //- Make patch face - neighbour cell distances - void makeDeltaCoeffs(scalarField&) const; - public: diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C index a7f9f6592e12b4e49dd9889749509873108c0a65..9717c96e6c10094984da271398f4535c7ce67226 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C @@ -60,27 +60,6 @@ void Foam::cyclicAMIFvPatch::makeWeights(scalarField& w) const } -void Foam::cyclicAMIFvPatch::makeDeltaCoeffs(scalarField& dc) const -{ - const cyclicAMIFvPatch& nbrPatch = neighbFvPatch(); - - const scalarField deltas(nf() & fvPatch::delta()); - - const scalarField nbrDeltas - ( - interpolate(nbrPatch.nf() & nbrPatch.fvPatch::delta()) - ); - - forAll(deltas, faceI) - { - scalar di = deltas[faceI]; - scalar dni = nbrDeltas[faceI]; - - dc[faceI] = 1.0/(di + dni); - } -} - - Foam::tmp<Foam::vectorField> Foam::cyclicAMIFvPatch::delta() const { const vectorField patchD(fvPatch::delta()); diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H index 3e29a2be42aacbd49754ea8542e91f08e7397ab4..c53518fb23a7b6266ded7ae2325d75607e5a797f 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H @@ -66,9 +66,6 @@ protected: //- Make patch weighting factors void makeWeights(scalarField&) const; - //- Make patch face - neighbour cell distances - void makeDeltaCoeffs(scalarField&) const; - public: diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.C index cbd3575ff1681eb048f1ac2b053a792a4c19ac44..9a2c1568a971c177c422301559c65385db661cea 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.C @@ -65,19 +65,6 @@ void processorFvPatch::makeWeights(scalarField& w) const } -void processorFvPatch::makeDeltaCoeffs(scalarField& dc) const -{ - if (Pstream::parRun()) - { - dc = (1.0 - weights())/(nf() & fvPatch::delta()); - } - else - { - dc = 1.0/(nf() & fvPatch::delta()); - } -} - - tmp<vectorField> processorFvPatch::delta() const { if (Pstream::parRun()) diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.H index 0058299e77204b2ceb87135fe56ce1661afa1f1f..d0f3b04c7ccf20cd8f88e437d8a6521d63e24db6 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.H @@ -65,9 +65,6 @@ protected: //- Make patch weighting factors void makeWeights(scalarField&) const; - //- Make patch face - neighbour cell distances - void makeDeltaCoeffs(scalarField&) const; - public: diff --git a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C index f9d1b8aca70042ad542626da8ce84168231816d1..1c9a6ca20c9371d440dfe980b6010b26f2ef6b28 100644 --- a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C @@ -150,12 +150,6 @@ void Foam::fvPatch::makeWeights(scalarField& w) const } -void Foam::fvPatch::makeDeltaCoeffs(scalarField& dc) const -{ - dc = 1.0/(nf() & delta()); -} - - void Foam::fvPatch::initMovePoints() {} diff --git a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H index d94363f05f62db8b65f51c36375bec94f751dcf1..be7ee1bda475e0e2210c11ae142b13a65b45b785 100644 --- a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H @@ -86,9 +86,6 @@ protected: //- Make patch weighting factors virtual void makeWeights(scalarField&) const; - //- Make patch face - neighbour cell distances - virtual void makeDeltaCoeffs(scalarField&) const; - //- Initialise the patches for moving points virtual void initMovePoints(); diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C index b8f108e73c29387b43cc3d6d811b92aa5dd70740..a7ea3a3374fa5e8cc125a2ee605b77455df238b3 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C @@ -122,20 +122,7 @@ Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter ); // Build the d-vectors - vectorField pd - ( - mesh.Sf().boundaryField()[patchi] - / ( - mesh.magSf().boundaryField()[patchi] - * mesh.deltaCoeffs().boundaryField()[patchi] - ) - ); - - if (!mesh.orthogonal()) - { - pd -= mesh.correctionVectors().boundaryField()[patchi] - /mesh.deltaCoeffs().boundaryField()[patchi]; - } + vectorField pd = CDweights.boundaryField()[patchi].patch().delta(); forAll(pLim, face) { diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C index 183a7684bd604c409c145da70ab4393e35e7bbc3..b513c50e4a060efdc55f61e1e8ddd27699f972bd 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.C @@ -111,20 +111,7 @@ Foam::linearUpwind<Type>::correction ); // Build the d-vectors - vectorField pd - ( - mesh.Sf().boundaryField()[patchi] - / ( - mesh.magSf().boundaryField()[patchi] - * mesh.deltaCoeffs().boundaryField()[patchi] - ) - ); - - if (!mesh.orthogonal()) - { - pd -= mesh.correctionVectors().boundaryField()[patchi] - /mesh.deltaCoeffs().boundaryField()[patchi]; - } + vectorField pd = Cf.boundaryField()[patchi].patch().delta(); forAll(pOwner, facei) { diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C index 747cd83f7f196a5d06fc4ebd5d9263ac957d2cc7..9feff9adaa1409eba853ce0b6d58b85ae7cc3f79 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C @@ -82,21 +82,17 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const const surfaceVectorField& Sf = mesh_.Sf(); const labelUList& owner = mesh_.owner(); + const labelUList& neighbour = mesh_.neighbour(); - // Build the d-vectors - surfaceVectorField d(Sf/(mesh_.magSf()*mesh_.deltaCoeffs())); - - if (!mesh_.orthogonal()) + forAll(owner, facei) { - d -= mesh_.correctionVectors()/mesh_.deltaCoeffs(); - } + label own = owner[facei]; + label nei = neighbour[facei]; - forAll(owner, faceI) - { - vector Cpf = Cf[faceI] - C[owner[faceI]]; + vector d = C[nei] - C[own]; + vector Cpf = Cf[facei] - C[own]; - SkewCorrVecs[faceI] = - Cpf - ((Sf[faceI] & Cpf)/(Sf[faceI] & d[faceI]))*d[faceI]; + SkewCorrVecs[facei] = Cpf - ((Sf[facei] & Cpf)/(Sf[facei] & d))*d; } @@ -115,7 +111,7 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const const labelUList& faceCells = p.faceCells(); const vectorField& patchFaceCentres = Cf.boundaryField()[patchI]; const vectorField& patchSf = Sf.boundaryField()[patchI]; - const vectorField& patchD = d.boundaryField()[patchI]; + const vectorField patchD = p.delta(); forAll(p, patchFaceI) { @@ -136,7 +132,7 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const if (Sf.internalField().size()) { - skewCoeff = max(mag(SkewCorrVecs)/mag(d)).value(); + skewCoeff = max(mag(SkewCorrVecs)*mesh_.deltaCoeffs()).value(); } if (debug) @@ -182,7 +178,7 @@ const Foam::surfaceVectorField& Foam::skewCorrectionVectors::operator()() const if (!skew()) { FatalErrorIn("skewCorrectionVectors::operator()()") - << "Cannot return correctionVectors; mesh is not skewed" + << "Cannot return skewCorrectionVectors; mesh is not skewed" << abort(FatalError); } diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C index a18ed2367cdebf497e01d8f8fe20af452d17eac1..fe0643c1951bf15056fdf64f6f46039ba560a6f8 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C @@ -31,43 +31,38 @@ Description #include "surfaceFields.H" #include "demandDrivenData.H" #include "coupledFvPatch.H" -#include "unitConversion.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -defineTypeNameAndDebug(surfaceInterpolation, 0); +defineTypeNameAndDebug(Foam::surfaceInterpolation, 0); // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // -void surfaceInterpolation::clearOut() +void Foam::surfaceInterpolation::clearOut() { - deleteDemandDrivenData(weightingFactors_); - deleteDemandDrivenData(differenceFactors_); - deleteDemandDrivenData(correctionVectors_); + deleteDemandDrivenData(weights_); + deleteDemandDrivenData(deltaCoeffs_); + deleteDemandDrivenData(nonOrthDeltaCoeffs_); + deleteDemandDrivenData(nonOrthCorrectionVectors_); } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -surfaceInterpolation::surfaceInterpolation(const fvMesh& fvm) +Foam::surfaceInterpolation::surfaceInterpolation(const fvMesh& fvm) : mesh_(fvm), - weightingFactors_(NULL), - differenceFactors_(NULL), - orthogonal_(false), - correctionVectors_(NULL) + weights_(NULL), + deltaCoeffs_(NULL), + nonOrthDeltaCoeffs_(NULL), + nonOrthCorrectionVectors_(NULL) {} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // -surfaceInterpolation::~surfaceInterpolation() +Foam::surfaceInterpolation::~surfaceInterpolation() { clearOut(); } @@ -75,66 +70,67 @@ surfaceInterpolation::~surfaceInterpolation() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const surfaceScalarField& surfaceInterpolation::weights() const +const Foam::surfaceScalarField& +Foam::surfaceInterpolation::weights() const { - if (!weightingFactors_) + if (!weights_) { makeWeights(); } - return (*weightingFactors_); + return (*weights_); } -const surfaceScalarField& surfaceInterpolation::deltaCoeffs() const +const Foam::surfaceScalarField& +Foam::surfaceInterpolation::deltaCoeffs() const { - if (!differenceFactors_) + if (!deltaCoeffs_) { makeDeltaCoeffs(); } - return (*differenceFactors_); + return (*deltaCoeffs_); } -bool surfaceInterpolation::orthogonal() const +const Foam::surfaceScalarField& +Foam::surfaceInterpolation::nonOrthDeltaCoeffs() const { - if (orthogonal_ == false && !correctionVectors_) + if (!nonOrthDeltaCoeffs_) { - makeCorrectionVectors(); + makeNonOrthDeltaCoeffs(); } - return orthogonal_; + return (*nonOrthDeltaCoeffs_); } -const surfaceVectorField& surfaceInterpolation::correctionVectors() const +const Foam::surfaceVectorField& +Foam::surfaceInterpolation::nonOrthCorrectionVectors() const { - if (orthogonal()) + if (!nonOrthCorrectionVectors_) { - FatalErrorIn("surfaceInterpolation::correctionVectors()") - << "cannot return correctionVectors; mesh is orthogonal" - << abort(FatalError); + makeNonOrthCorrectionVectors(); } - return (*correctionVectors_); + return (*nonOrthCorrectionVectors_); } // Do what is neccessary if the mesh has moved -bool surfaceInterpolation::movePoints() +bool Foam::surfaceInterpolation::movePoints() { - deleteDemandDrivenData(weightingFactors_); - deleteDemandDrivenData(differenceFactors_); - - orthogonal_ = false; - deleteDemandDrivenData(correctionVectors_); + deleteDemandDrivenData(weights_); + deleteDemandDrivenData(deltaCoeffs_); + deleteDemandDrivenData(nonOrthDeltaCoeffs_); + deleteDemandDrivenData(nonOrthCorrectionVectors_); return true; } -void surfaceInterpolation::makeWeights() const +void Foam::surfaceInterpolation::makeWeights() const { if (debug) { @@ -143,20 +139,18 @@ void surfaceInterpolation::makeWeights() const << endl; } - - weightingFactors_ = new surfaceScalarField + weights_ = new surfaceScalarField ( IOobject ( - "weightingFactors", + "weights", mesh_.pointsInstance(), mesh_ ), mesh_, dimless ); - surfaceScalarField& weightingFactors = *weightingFactors_; - + surfaceScalarField& weights = *weights_; // Set local references to mesh data // (note that we should not use fvMesh sliced fields at this point yet @@ -170,7 +164,7 @@ void surfaceInterpolation::makeWeights() const const vectorField& Sf = mesh_.faceAreas(); // ... and reference to the internal field of the weighting factors - scalarField& w = weightingFactors.internalField(); + scalarField& w = weights.internalField(); forAll(owner, facei) { @@ -188,11 +182,10 @@ void surfaceInterpolation::makeWeights() const { mesh_.boundary()[patchi].makeWeights ( - weightingFactors.boundaryField()[patchi] + weights.boundaryField()[patchi] ); } - if (debug) { Info<< "surfaceInterpolation::makeWeights() : " @@ -202,7 +195,7 @@ void surfaceInterpolation::makeWeights() const } -void surfaceInterpolation::makeDeltaCoeffs() const +void Foam::surfaceInterpolation::makeDeltaCoeffs() const { if (debug) { @@ -215,18 +208,63 @@ void surfaceInterpolation::makeDeltaCoeffs() const // needed to make sure deltaCoeffs are calculated for parallel runs. weights(); - differenceFactors_ = new surfaceScalarField + deltaCoeffs_ = new surfaceScalarField ( IOobject ( - "differenceFactors_", + "deltaCoeffs", mesh_.pointsInstance(), mesh_ ), mesh_, dimless/dimLength ); - surfaceScalarField& DeltaCoeffs = *differenceFactors_; + surfaceScalarField& DeltaCoeffs = *deltaCoeffs_; + + + // Set local references to mesh data + const volVectorField& C = mesh_.C(); + const labelUList& owner = mesh_.owner(); + const labelUList& neighbour = mesh_.neighbour(); + + forAll(owner, facei) + { + DeltaCoeffs[facei] = 1.0/mag(C[neighbour[facei]] - C[owner[facei]]); + } + + forAll(DeltaCoeffs.boundaryField(), patchi) + { + DeltaCoeffs.boundaryField()[patchi] = + 1.0/mag(mesh_.boundary()[patchi].delta()); + } +} + + +void Foam::surfaceInterpolation::makeNonOrthDeltaCoeffs() const +{ + if (debug) + { + Info<< "surfaceInterpolation::makeNonOrthDeltaCoeffs() : " + << "Constructing differencing factors array for face gradient" + << endl; + } + + // Force the construction of the weighting factors + // needed to make sure deltaCoeffs are calculated for parallel runs. + weights(); + + nonOrthDeltaCoeffs_ = new surfaceScalarField + ( + IOobject + ( + "nonOrthDeltaCoeffs", + mesh_.pointsInstance(), + mesh_ + ), + mesh_, + dimless/dimLength + ); + surfaceScalarField& nonOrthDeltaCoeffs = *nonOrthDeltaCoeffs_; // Set local references to mesh data @@ -242,49 +280,49 @@ void surfaceInterpolation::makeDeltaCoeffs() const vector unitArea = Sf[facei]/magSf[facei]; // Standard cell-centre distance form - //DeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta); + //NonOrthDeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta); // Slightly under-relaxed form - //DeltaCoeffs[facei] = 1.0/mag(delta); + //NonOrthDeltaCoeffs[facei] = 1.0/mag(delta); // More under-relaxed form - //DeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL); + //NonOrthDeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL); // Stabilised form for bad meshes - DeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta)); + nonOrthDeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta)); } - forAll(DeltaCoeffs.boundaryField(), patchi) + forAll(nonOrthDeltaCoeffs.boundaryField(), patchi) { - mesh_.boundary()[patchi].makeDeltaCoeffs - ( - DeltaCoeffs.boundaryField()[patchi] - ); + vectorField delta = mesh_.boundary()[patchi].delta(); + + nonOrthDeltaCoeffs.boundaryField()[patchi] = + 1.0/max(mesh_.boundary()[patchi].nf() & delta, 0.05*mag(delta)); } } -void surfaceInterpolation::makeCorrectionVectors() const +void Foam::surfaceInterpolation::makeNonOrthCorrectionVectors() const { if (debug) { - Info<< "surfaceInterpolation::makeCorrectionVectors() : " + Info<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : " << "Constructing non-orthogonal correction vectors" << endl; } - correctionVectors_ = new surfaceVectorField + nonOrthCorrectionVectors_ = new surfaceVectorField ( IOobject ( - "correctionVectors", + "nonOrthCorrectionVectors", mesh_.pointsInstance(), mesh_ ), mesh_, dimless ); - surfaceVectorField& corrVecs = *correctionVectors_; + surfaceVectorField& corrVecs = *nonOrthCorrectionVectors_; // Set local references to mesh data const volVectorField& C = mesh_.C(); @@ -292,14 +330,14 @@ void surfaceInterpolation::makeCorrectionVectors() const const labelUList& neighbour = mesh_.neighbour(); const surfaceVectorField& Sf = mesh_.Sf(); const surfaceScalarField& magSf = mesh_.magSf(); - const surfaceScalarField& DeltaCoeffs = deltaCoeffs(); + const surfaceScalarField& NonOrthDeltaCoeffs = nonOrthDeltaCoeffs(); forAll(owner, facei) { vector unitArea = Sf[facei]/magSf[facei]; vector delta = C[neighbour[facei]] - C[owner[facei]]; - corrVecs[facei] = unitArea - delta*DeltaCoeffs[facei]; + corrVecs[facei] = unitArea - delta*NonOrthDeltaCoeffs[facei]; } // Boundary correction vectors set to zero for boundary patches @@ -308,18 +346,18 @@ void surfaceInterpolation::makeCorrectionVectors() const forAll(corrVecs.boundaryField(), patchi) { - fvsPatchVectorField& patchcorrVecs = corrVecs.boundaryField()[patchi]; + fvsPatchVectorField& patchCorrVecs = corrVecs.boundaryField()[patchi]; - if (!patchcorrVecs.coupled()) + if (!patchCorrVecs.coupled()) { - patchcorrVecs = vector::zero; + patchCorrVecs = vector::zero; } else { - const fvsPatchScalarField& patchDeltaCoeffs - = DeltaCoeffs.boundaryField()[patchi]; + const fvsPatchScalarField& patchNonOrthDeltaCoeffs + = NonOrthDeltaCoeffs.boundaryField()[patchi]; - const fvPatch& p = patchcorrVecs.patch(); + const fvPatch& p = patchCorrVecs.patch(); const vectorField patchDeltas(mesh_.boundary()[patchi].delta()); @@ -331,60 +369,19 @@ void surfaceInterpolation::makeCorrectionVectors() const const vector& delta = patchDeltas[patchFacei]; - patchcorrVecs[patchFacei] = - unitArea - delta*patchDeltaCoeffs[patchFacei]; + patchCorrVecs[patchFacei] = + unitArea - delta*patchNonOrthDeltaCoeffs[patchFacei]; } } } - scalar NonOrthogCoeff = 0.0; - - // Calculate the non-orthogonality for meshes with 1 face or more - if (returnReduce(magSf.size(), sumOp<label>()) > 0) - { - NonOrthogCoeff = radToDeg - ( - asin - ( - min - ( - (sum(magSf*mag(corrVecs))/sum(magSf)).value(), - 1.0 - ) - ) - ); - } - - if (debug) - { - Info<< "surfaceInterpolation::makeCorrectionVectors() : " - << "non-orthogonality coefficient = " << NonOrthogCoeff << " deg." - << endl; - } - - //NonOrthogCoeff = 0.0; - - if (NonOrthogCoeff < 0.1) - { - orthogonal_ = true; - deleteDemandDrivenData(correctionVectors_); - } - else - { - orthogonal_ = false; - } - if (debug) { - Info<< "surfaceInterpolation::makeCorrectionVectors() : " + Info<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : " << "Finished constructing non-orthogonal correction vectors" << endl; } } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - // ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.H b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.H index 17333be4abe1ab54768d88e611a719a7229e1e12..4a8601cb32a4d4ead857db07cf3cf2314ebb0f6b 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.H @@ -59,17 +59,17 @@ class surfaceInterpolation // Demand-driven data - //- Central-differencing weighting factors - mutable surfaceScalarField* weightingFactors_; + //- Linear difference weighting factors + mutable surfaceScalarField* weights_; - //- Face-gradient difference factors - mutable surfaceScalarField* differenceFactors_; + //- Cell-centre difference coefficients + mutable surfaceScalarField* deltaCoeffs_; - //- Is mesh orthogonal - mutable bool orthogonal_; + //- Non-orthogonal cell-centre difference coefficients + mutable surfaceScalarField* nonOrthDeltaCoeffs_; //- Non-orthogonality correction vectors - mutable surfaceVectorField* correctionVectors_; + mutable surfaceVectorField* nonOrthCorrectionVectors_; // Private Member Functions @@ -80,8 +80,11 @@ class surfaceInterpolation //- Construct face-gradient difference factors void makeDeltaCoeffs() const; + //- Construct face-gradient difference factors + void makeNonOrthDeltaCoeffs() const; + //- Construct non-orthogonality correction vectors - void makeCorrectionVectors() const; + void makeNonOrthCorrectionVectors() const; protected: @@ -112,17 +115,18 @@ public: // Member functions - //- Return reference to weighting factors array + //- Return reference to linear difference weighting factors const surfaceScalarField& weights() const; - //- Return reference to difference factors array + //- Return reference to cell-centre difference coefficients const surfaceScalarField& deltaCoeffs() const; - //- Return whether mesh is orthogonal or not - bool orthogonal() const; + //- Return reference to non-orthogonal cell-centre difference + // coefficients + const surfaceScalarField& nonOrthDeltaCoeffs() const; - //- Return reference to non-orthogonality correction vectors array - const surfaceVectorField& correctionVectors() const; + //- Return reference to non-orthogonality correction vectors + const surfaceVectorField& nonOrthCorrectionVectors() const; //- Do what is neccessary if the mesh has moved bool movePoints(); diff --git a/src/fvMotionSolver/Make/files b/src/fvMotionSolver/Make/files index 87e8e39c255b2b319371f1f10d409132224e50f8..bea2f1065699047d4470169b4a9b92cbb1e1e0c7 100644 --- a/src/fvMotionSolver/Make/files +++ b/src/fvMotionSolver/Make/files @@ -33,5 +33,6 @@ pointPatchFields/derived/oscillatingDisplacement/oscillatingDisplacementPointPat pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C +pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.C LIB = $(FOAM_LIBBIN)/libfvMotionSolvers diff --git a/src/fvMotionSolver/pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.C new file mode 100644 index 0000000000000000000000000000000000000000..86956517456ec7dfd5b3b282b892c6fb70473090 --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.C @@ -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/>. + +\*---------------------------------------------------------------------------*/ + +#include "waveDisplacementPointPatchVectorField.H" +#include "pointPatchFields.H" +#include "addToRunTimeSelectionTable.H" +#include "Time.H" +#include "polyMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +waveDisplacementPointPatchVectorField:: +waveDisplacementPointPatchVectorField +( + const pointPatch& p, + const DimensionedField<vector, pointMesh>& iF +) +: + fixedValuePointPatchField<vector>(p, iF), + amplitude_(vector::zero), + omega_(0.0), + waveNumber_(vector::zero) +{} + + +waveDisplacementPointPatchVectorField:: +waveDisplacementPointPatchVectorField +( + const pointPatch& p, + const DimensionedField<vector, pointMesh>& iF, + const dictionary& dict +) +: + fixedValuePointPatchField<vector>(p, iF, dict), + amplitude_(dict.lookup("amplitude")), + omega_(readScalar(dict.lookup("omega"))), + waveNumber_(dict.lookupOrDefault<vector>("waveLength", vector::zero)) +{ + if (!dict.found("value")) + { + updateCoeffs(); + } +} + + +waveDisplacementPointPatchVectorField:: +waveDisplacementPointPatchVectorField +( + const waveDisplacementPointPatchVectorField& ptf, + const pointPatch& p, + const DimensionedField<vector, pointMesh>& iF, + const pointPatchFieldMapper& mapper +) +: + fixedValuePointPatchField<vector>(ptf, p, iF, mapper), + amplitude_(ptf.amplitude_), + omega_(ptf.omega_), + waveNumber_(ptf.waveNumber_) +{} + + +waveDisplacementPointPatchVectorField:: +waveDisplacementPointPatchVectorField +( + const waveDisplacementPointPatchVectorField& ptf, + const DimensionedField<vector, pointMesh>& iF +) +: + fixedValuePointPatchField<vector>(ptf, iF), + amplitude_(ptf.amplitude_), + omega_(ptf.omega_), + waveNumber_(ptf.waveNumber_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void waveDisplacementPointPatchVectorField::updateCoeffs() +{ + if (this->updated()) + { + return; + } + + const polyMesh& mesh = this->dimensionedInternalField().mesh()(); + const Time& t = mesh.time(); + + const scalarField points( waveNumber_ & patch().localPoints()); + + Field<vector>::operator= + ( + amplitude_*cos(omega_*t.value() - points) + ); + + fixedValuePointPatchField<vector>::updateCoeffs(); +} + + +void waveDisplacementPointPatchVectorField::write(Ostream& os) const +{ + pointPatchField<vector>::write(os); + os.writeKeyword("amplitude") + << amplitude_ << token::END_STATEMENT << nl; + os.writeKeyword("omega") + << omega_ << token::END_STATEMENT << nl; + os.writeKeyword("waveNumber") + << waveNumber_ << token::END_STATEMENT << nl; + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePointPatchTypeField +( + pointPatchVectorField, + waveDisplacementPointPatchVectorField +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/fvMotionSolver/pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.H b/src/fvMotionSolver/pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.H new file mode 100644 index 0000000000000000000000000000000000000000..70d2f6c7f76be1d1b61d59b14dc50b601eedd2d2 --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.H @@ -0,0 +1,149 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::waveDisplacementPointPatchVectorField + +Description + Foam::waveDisplacementPointPatchVectorField + +SourceFiles + waveDisplacementPointPatchVectorField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef waveDisplacementPointPatchVectorField_H +#define waveDisplacementPointPatchVectorField_H + +#include "fixedValuePointPatchField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class waveDisplacementPointPatchVectorField Declaration +\*---------------------------------------------------------------------------*/ + +class waveDisplacementPointPatchVectorField +: + public fixedValuePointPatchField<vector> +{ + // Private data + + vector amplitude_; + scalar omega_; + vector waveNumber_; + + +public: + + //- Runtime type information + TypeName("waveDisplacement"); + + + // Constructors + + //- Construct from patch and internal field + waveDisplacementPointPatchVectorField + ( + const pointPatch&, + const DimensionedField<vector, pointMesh>& + ); + + //- Construct from patch, internal field and dictionary + waveDisplacementPointPatchVectorField + ( + const pointPatch&, + const DimensionedField<vector, pointMesh>&, + const dictionary& + ); + + //- Construct by mapping given patchField<vector> onto a new patch + waveDisplacementPointPatchVectorField + ( + const waveDisplacementPointPatchVectorField&, + const pointPatch&, + const DimensionedField<vector, pointMesh>&, + const pointPatchFieldMapper& + ); + + //- Construct and return a clone + virtual autoPtr<pointPatchField<vector> > clone() const + { + return autoPtr<pointPatchField<vector> > + ( + new waveDisplacementPointPatchVectorField + ( + *this + ) + ); + } + + //- Construct as copy setting internal field reference + waveDisplacementPointPatchVectorField + ( + const waveDisplacementPointPatchVectorField&, + const DimensionedField<vector, pointMesh>& + ); + + //- Construct and return a clone setting internal field reference + virtual autoPtr<pointPatchField<vector> > clone + ( + const DimensionedField<vector, pointMesh>& iF + ) const + { + return autoPtr<pointPatchField<vector> > + ( + new waveDisplacementPointPatchVectorField + ( + *this, + iF + ) + ); + } + + + // Member functions + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // 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/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C b/src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWaveName.C similarity index 55% rename from applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C rename to src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWaveName.C index ea8f509a292d1761a0e9665f17153aa0171c87fb..5d3f3e1d581a7219bb0f1b9ce2f8f39b86753fc5 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C +++ b/src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWaveName.C @@ -23,45 +23,10 @@ License \*---------------------------------------------------------------------------*/ -#include "patchPointEdgeCirculator.H" +#include "PatchEdgeFaceWave.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 -) -{ - 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; -} - +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/postProcessing/functionObjects/field/streamLine/controlDict b/src/postProcessing/functionObjects/field/streamLine/controlDict new file mode 100644 index 0000000000000000000000000000000000000000..4d51dd82e1f758f5af36f68e927088bb87b14d90 --- /dev/null +++ b/src/postProcessing/functionObjects/field/streamLine/controlDict @@ -0,0 +1,99 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: 2.0.0 | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application icoFoam; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 0.5; + +deltaT 0.005; + +writeControl timeStep; + +writeInterval 20; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression off; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable true; + +functions +{ + streamLines + { + type streamLine; + + // Where to load it from (if not already in solver) + functionObjectLibs ("libfieldFunctionObjects.so"); + + // Output every + outputControl outputTime; + // outputInterval 10; + + setFormat vtk; //gnuplot; //xmgr; //raw; //jplot; + + // Velocity field to use for tracking. + UName U; + + // Interpolation method. Default is cellPoint. See sampleDict. + //interpolationScheme pointMVC; + + // Tracked forwards (+U) or backwards (-U) + trackForward true; + + // Names of fields to sample. Should contain above velocity field! + fields (p U); + + // Steps particles can travel before being removed + lifeTime 10000; + + // Number of steps per cell (estimate). Set to 1 to disable subcycling. + nSubCycle 5; + + // Cloud name to use + cloudName particleTracks; + + // Seeding method. See the sampleSets in sampleDict. + seedSampleSet uniform; //cloud;//triSurfaceMeshPointSet; + + uniformCoeffs + { + type uniform; + axis x; //distance; + + start (-0.0205 0.0001 0.00001); + end (-0.0205 0.0005 0.00001); + nPoints 100; + } + } +} + + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLine.C b/src/postProcessing/functionObjects/field/streamLine/streamLine.C index fd57762a9f3ba73ab9e801ea05667fb898a75250..03a5b8b539c4b803a9b050dda20a0f9cb3e43ef6 100644 --- a/src/postProcessing/functionObjects/field/streamLine/streamLine.C +++ b/src/postProcessing/functionObjects/field/streamLine/streamLine.C @@ -33,6 +33,7 @@ License #include "sampledSet.H" #include "globalIndex.H" #include "mapDistribute.H" +#include "interpolationCellPoint.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -75,9 +76,9 @@ void Foam::streamLine::track() // Read or lookup fields PtrList<volScalarField> vsFlds; - PtrList<interpolationCellPoint<scalar> > vsInterp; + PtrList<interpolation<scalar> > vsInterp; PtrList<volVectorField> vvFlds; - PtrList<interpolationCellPoint<vector> > vvInterp; + PtrList<interpolation<vector> > vvInterp; label UIndex = -1; @@ -95,13 +96,29 @@ void Foam::streamLine::track() vsInterp.setSize(vsFlds.size()); forAll(vsFlds, i) { - vsInterp.set(i, new interpolationCellPoint<scalar>(vsFlds[i])); + vsInterp.set + ( + i, + interpolation<scalar>::New + ( + interpolationScheme_, + vsFlds[i] + ) + ); } ReadFields(mesh, objects, vvFlds); vvInterp.setSize(vvFlds.size()); forAll(vvFlds, i) { - vvInterp.set(i, new interpolationCellPoint<vector>(vvFlds[i])); + vvInterp.set + ( + i, + interpolation<vector>::New + ( + interpolationScheme_, + vvFlds[i] + ) + ); } } else @@ -146,7 +163,12 @@ void Foam::streamLine::track() vsInterp.set ( nScalar++, - new interpolationCellPoint<scalar>(f) + //new interpolationCellPoint<scalar>(f) + interpolation<scalar>::New + ( + interpolationScheme_, + f + ) ); } else if (mesh.foundObject<volVectorField>(fields_[i])) @@ -164,7 +186,12 @@ void Foam::streamLine::track() vvInterp.set ( nVector++, - new interpolationCellPoint<vector>(f) + //new interpolationCellPoint<vector>(f) + interpolation<vector>::New + ( + interpolationScheme_, + f + ) ); } } @@ -291,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) @@ -307,6 +349,15 @@ void Foam::streamLine::read(const dictionary& dict) nSubCycle_ = 1; } + interpolationScheme_ = dict.lookupOrDefault + ( + "interpolationScheme", + interpolationCellPoint<scalar>::typeName + ); + + //Info<< typeName << " using interpolation " << interpolationScheme_ + // << endl; + cloudName_ = dict.lookupOrDefault<word>("cloudName", "streamLine"); dict.lookup("seedSampleSet") >> seedSet_; diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLine.H b/src/postProcessing/functionObjects/field/streamLine/streamLine.H index e43fdd7936784446878190ad507b03180fdae526..9fc981034c4f0b6a8a3cbad0332d5cfcd1c1742d 100644 --- a/src/postProcessing/functionObjects/field/streamLine/streamLine.H +++ b/src/postProcessing/functionObjects/field/streamLine/streamLine.H @@ -86,6 +86,9 @@ class streamLine //- Field to transport particle with word UName_; + //- Interpolation scheme to use + word interpolationScheme_; + //- Whether to use +u or -u bool trackForward_; diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H index d01353c46f7312f5565619731c1ce8467e5ebd9c..b1ee9acf8847e5f54b6856f2247aff972ae21311 100644 --- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H +++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H @@ -38,7 +38,7 @@ SourceFiles #include "particle.H" #include "autoPtr.H" -#include "interpolationCellPoint.H" +#include "interpolation.H" #include "vectorList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -68,8 +68,8 @@ public: public: - const PtrList<interpolationCellPoint<scalar> >& vsInterp_; - const PtrList<interpolationCellPoint<vector> >& vvInterp_; + const PtrList<interpolation<scalar> >& vsInterp_; + const PtrList<interpolation<vector> >& vvInterp_; const label UIndex_; const bool trackForward_; const label nSubCycle_; @@ -84,8 +84,8 @@ public: trackingData ( Cloud<streamLineParticle>& cloud, - const PtrList<interpolationCellPoint<scalar> >& vsInterp, - const PtrList<interpolationCellPoint<vector> >& vvInterp, + const PtrList<interpolation<scalar> >& vsInterp, + const PtrList<interpolation<vector> >& vvInterp, const label UIndex, const bool trackForward, const label nSubCycle, diff --git a/src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModelCollection.C b/src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModelCollection.C index 4ec1871f88875a2983dcb521c3f9a803450a3784..4bc3ca977e8fa08ec1f1ed0af079f93058c746dd 100644 --- a/src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModelCollection.C +++ b/src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModelCollection.C @@ -26,6 +26,16 @@ License #include "pyrolysisModelCollection.H" #include "volFields.H" +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTemplateTypeNameAndDebug +( + Foam::IOPtrList<Foam::regionModels::pyrolysisModels::pyrolysisModel>, + 0 +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + namespace Foam { namespace regionModels @@ -33,11 +43,6 @@ namespace regionModels namespace pyrolysisModels { -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -defineTemplateTypeNameAndDebug(IOPtrList<pyrolysisModel>, 0); - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // pyrolysisModelCollection::pyrolysisModelCollection diff --git a/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C b/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C index 394e35e536d7ff79ddf707e81590ff28275aae6c..d5707fe8fd5fba5deea7d145356688aa18b60f55 100644 --- a/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C +++ b/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C @@ -49,7 +49,8 @@ const ThermoType& Foam::multiComponentMixture<ThermoType>::constructSpeciesData template<class ThermoType> void Foam::multiComponentMixture<ThermoType>::correctMassFractions() { - volScalarField Yt("Yt", Y_[0]); + // It changes Yt patches to "calculated" + volScalarField Yt("Yt", 1.0*Y_[0]); for (label n=1; n<Y_.size(); n++) { diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C index 6509ccd1b71fa5cf1e3fd1cfb5a19523196c7c22..a3d313e30285ffb9eaf10ac5c4f658518f6abec0 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C @@ -86,7 +86,7 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField z_(dict.lookup("z")), z0_(readScalar(dict.lookup("z0"))), kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)), - zGround_(readScalar(dict.lookup("zGround"))) + zGround_("zGround", dict, p.size()) { if (mag(z_) < SMALL) { diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H index 6b4f7a8d453c04f666dec43019ce0d961114eaec..ab08c621c95da1aa092c751da6192598e066a7d8 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H @@ -90,7 +90,7 @@ class atmBoundaryLayerInletEpsilonFvPatchScalarField const scalar kappa_; //- Minimum corrdinate value in z direction - const scalar zGround_; + const scalarField zGround_; public: diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C index db71b28d2d1d9a75b382f06b1ebb6ef2dc33f748..94eefec0f673c647ada1e3391b19aebede30280e 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C @@ -115,7 +115,7 @@ atmBoundaryLayerInletVelocityFvPatchVectorField n_ /= mag(n_); z_ /= mag(z_); - Ustar_ = kappa_*Uref_/(log((Href_ + z0_)/min(z0_ , 0.001))); + Ustar_ = kappa_*Uref_/(log((Href_ + z0_)/max(z0_ , 0.001))); evaluate(); } @@ -150,9 +150,11 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::updateCoeffs() forAll(coord, i) { - if ((coord[i] - zGround_) < Href_) + if ((coord[i] - zGround_[i]) < Href_) { - Un[i] = (Ustar_/kappa_)*log((coord[i] - zGround_ + z0_)/z0_); + Un[i] = + (Ustar_/kappa_) + * log((coord[i] - zGround_[i] + z0_)/max(z0_, 0.001)); } else { @@ -181,8 +183,7 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::write(Ostream& os) const << Uref_ << token::END_STATEMENT << nl; os.writeKeyword("Href") << Href_ << token::END_STATEMENT << nl; - os.writeKeyword("zGround") - << zGround_ << token::END_STATEMENT << nl; + zGround_.writeEntry("zGround", os) ; writeEntry("value", os); } diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H index c2f4e4c2f153a0130c454281c62442574d15d4ed..aa383742b05372ca025a49e7d7a05783fd7dff40 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H @@ -113,7 +113,7 @@ class atmBoundaryLayerInletVelocityFvPatchVectorField const scalar Href_; //- Minimum corrdinate value in z direction - const scalar zGround_; + const scalarField zGround_; public: diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/N2 b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/N2 index db64580e610171ec989eff31d2a2be8a12222e2b..af44210d0c860629e25902f15079216ca7dfcc32 100644 --- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/N2 +++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/N2 @@ -21,32 +21,27 @@ boundaryField { top { - type inletOutlet; - inletValue $internalField; - value $internalField; + type calculated; } ground { - type zeroGradient; + type calculated; } sides { - type inletOutlet; - inletValue $internalField; - value $internalField; + type calculated; } burner { - type fixedValue; - value uniform 0; + type calculated; } "(region0_to.*)" { - type zeroGradient; + type calculated; } } diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/N2 b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/N2 index 5dd08b790392ab39c5cb8386c1668cd9346c762d..ff1a3b1c4207e80ead1c2554cd39958ebdc88854 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/N2 +++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/N2 @@ -23,24 +23,19 @@ boundaryField { outlet { - type inletOutlet; - inletValue $internalField; - value $internalField; + type calculated; } sides { - type inletOutlet; - inletValue $internalField; - value $internalField; + type calculated; } base { - type zeroGradient; + type calculated; } inlet { - type fixedValue; - value uniform 0; + type calculated; } frontBack { diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/N2 b/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/N2 index 774ae2ad307ee28cf1e8f37141c857230844a817..f44c142482410b3000f84ccb49f6ab243c6fae2e 100644 --- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/N2 +++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/N2 @@ -23,24 +23,19 @@ boundaryField { outlet { - type inletOutlet; - inletValue $internalField; - value $internalField; + type calculated; } sides { - type inletOutlet; - inletValue $internalField; - value $internalField; + type calculated; } base { - type zeroGradient; + type calculated; } inlet { - type fixedValue; - value uniform 0; + type calculated; } } diff --git a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/0/N2 b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/0/N2 index 88bb6fb5aac25e209a51777957ff333ef25e0ce2..ac188ad5051b369db8900495ef79d1cd4ae3e12f 100644 --- a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/0/N2 +++ b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/0/N2 @@ -23,20 +23,15 @@ boundaryField { fuel { - type fixedValue; - value uniform 0.0; + type calculated; } air { - type fixedValue; - value uniform 0.77; + type calculated; } outlet { - type inletOutlet; - inletValue uniform 0.77; - value uniform 0.77; - + type calculated; } frontAndBack { diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution index 346b9ea433a2ef3f6dadaa5f0185d45f62d09baf..8b5f9ccd5ff8c138a1d5a1251807e82e0483a2db 100644 --- a/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution +++ b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution @@ -190,7 +190,7 @@ relaxationFactors h 0.95; k 0.9; epsilon 0.9; - }w + } } relaxationFactors0 diff --git a/tutorials/incompressible/icoFoam/cavity/system/controlDict b/tutorials/incompressible/icoFoam/cavity/system/controlDict index d2fa09cc526b96181407595f79c0968a81be7c44..8016db385c2528d5d285b268fea33f6a968f984c 100644 --- a/tutorials/incompressible/icoFoam/cavity/system/controlDict +++ b/tutorials/incompressible/icoFoam/cavity/system/controlDict @@ -31,6 +31,10 @@ writeControl timeStep; writeInterval 20; +//- Additional dump every hour of cpuTime for e.g. restarts +//secondaryWriteControl cpuTime; +//secondaryWriteInterval 3600; + purgeWrite 0; writeFormat ascii; 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; diff --git a/tutorials/incompressible/windSimpleFoam/turbineSiting/0/include/ABLConditions b/tutorials/incompressible/windSimpleFoam/turbineSiting/0/include/ABLConditions index 495de0f2b8e0f06ac9abd33b8e82a3c8daccc008..be005f46a0764e366af7149e6670fb6ba8826638 100644 --- a/tutorials/incompressible/windSimpleFoam/turbineSiting/0/include/ABLConditions +++ b/tutorials/incompressible/windSimpleFoam/turbineSiting/0/include/ABLConditions @@ -13,5 +13,5 @@ z0 0.1; turbulentKE 1.3; windDirection (1 0 0); zDirection (0 0 1); -zGround 935.0; +zGround uniform 935.0; // ************************************************************************* // diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution index 9d4bfea5c024637160060544272ed4ca6a108e46..d09194f954dde0ac71cc6201ac1cb1f475455c51 100644 --- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution +++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution @@ -94,6 +94,7 @@ PIMPLE nCorrectors 2; nNonOrthogonalCorrectors 0; nAlphaCorr 2; + correctAlpha no; } relaxationFactors diff --git a/wmake/rules/linux64Icc/c++Opt b/wmake/rules/linux64Icc/c++Opt index 28a17f783cae1846508385ed35e278b34e6cea64..66638ffe984a81a2145ace8c5e44491fdb8fe2ae 100644 --- a/wmake/rules/linux64Icc/c++Opt +++ b/wmake/rules/linux64Icc/c++Opt @@ -1,3 +1,2 @@ c++DBUG = -#c++OPT = -xSSE3 -O3 -no-prec-div -c++OPT = -xSSE3 -O1 -no-prec-div +c++OPT = -xSSE3 -O2 -no-prec-div diff --git a/wmake/rules/linuxIcc/c++Opt b/wmake/rules/linuxIcc/c++Opt index 62f12c3eb533c5902776749d3e3916f78fe95790..66638ffe984a81a2145ace8c5e44491fdb8fe2ae 100644 --- a/wmake/rules/linuxIcc/c++Opt +++ b/wmake/rules/linuxIcc/c++Opt @@ -1,5 +1,2 @@ c++DBUG = -#c++OPT = -O3 -xP -no-prec-div -c++OPT = -ansi-alias -O3 -ftz -fno-alias \ - -fargument-noalias-global \ - -unroll0 +c++OPT = -xSSE3 -O2 -no-prec-div diff --git a/wmake/wmakeScheduler b/wmake/wmakeScheduler index aa6fb1ec81bd963412192fff74fcbda92d50a46e..f614e73758b03de205de3f69022053b56cf04b13 100755 --- a/wmake/wmakeScheduler +++ b/wmake/wmakeScheduler @@ -151,7 +151,8 @@ set -o pipefail # colourPipe() { - if [ "$1" ] + + if tty -s <&1 # [ "$1" ] then ( while read line