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.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C deleted file mode 100644 index ea8f509a292d1761a0e9665f17153aa0171c87fb..0000000000000000000000000000000000000000 --- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C +++ /dev/null @@ -1,67 +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/>. - -\*---------------------------------------------------------------------------*/ - -#include "patchPointEdgeCirculator.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; -} - - -// ************************************************************************* // 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; -} - - -// ************************************************************************* //