diff --git a/src/dynamicMesh/Make/files b/src/dynamicMesh/Make/files index a7f72ef9fa..d7310029a7 100644 --- a/src/dynamicMesh/Make/files +++ b/src/dynamicMesh/Make/files @@ -37,6 +37,7 @@ polyTopoChange/polyTopoChange/removePoints.C polyTopoChange/polyTopoChange/combineFaces.C polyTopoChange/polyTopoChange/localPointRegion.C polyTopoChange/polyTopoChange/duplicatePoints.C +polyTopoChange/polyTopoChange/tetDecomposer.C slidingInterface/slidingInterface.C slidingInterface/slidingInterfaceProjectPoints.C diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C new file mode 100644 index 0000000000..e8948c4fad --- /dev/null +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C @@ -0,0 +1,713 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 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 . + +\*---------------------------------------------------------------------------*/ + +#include "tetDecomposer.H" +#include "meshTools.H" +#include "polyMesh.H" +#include "polyTopoChange.H" +#include "mapPolyMesh.H" +#include "OFstream.H" +#include "EdgeMap.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(Foam::tetDecomposer, 0); + +namespace Foam +{ + template<> + const char* NamedEnum::names[] = + { + "faceCentre", + "faceDiagonal" + }; + + const NamedEnum + tetDecomposer::decompositionTypeNames; +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::tetDecomposer::modifyFace +( + polyTopoChange& meshMod, + const face& f, + const label faceI, + const label own, + const label nei, + const label patchI, + const label zoneI, + const bool zoneFlip +) const +{ + // First usage of face. Modify. + if (nei == -1 || own < nei) + { + meshMod.modifyFace + ( + f, // modified face + faceI, // label of face + own, // owner + nei, // neighbour + false, // face flip + patchI, // patch for face + zoneI, // zone for face + zoneFlip // face flip in zone + ); + } + else + { + meshMod.modifyFace + ( + f.reverseFace(), // modified face + faceI, // label of face + nei, // owner + own, // neighbour + true, // face flip + patchI, // patch for face + zoneI, // zone for face + !zoneFlip // face flip in zone + ); + } +} + + +void Foam::tetDecomposer::addFace +( + polyTopoChange& meshMod, + const face& f, + const label own, + const label nei, + const label masterPointID, + const label masterEdgeID, + const label masterFaceID, + const label patchI, + const label zoneI, + const bool zoneFlip +) const +{ + // Second or more usage of face. Add. + if (nei == -1 || own < nei) + { + meshMod.addFace + ( + f, // modified face + own, // owner + nei, // neighbour + masterPointID, // master point + masterEdgeID, // master edge + masterFaceID, // master face + false, // face flip + patchI, // patch for face + zoneI, // zone for face + zoneFlip // face flip in zone + ); + } + else + { + meshMod.addFace + ( + f.reverseFace(), // modified face + nei, // owner + own, // neighbour + masterPointID, // master point + masterEdgeID, // master edge + masterFaceID, // master face + true, // face flip + patchI, // patch for face + zoneI, // zone for face + !zoneFlip // face flip in zone + ); + } +} + + +// Work out triangle index given the starting vertex in the face +Foam::label Foam::tetDecomposer::triIndex(const label faceI, const label fp) +const +{ + const face& f = mesh_.faces()[faceI]; + const label fp0 = mesh_.tetBasePtIs()[faceI]; + + // Work out triangle index on this face + label thisTriI; + if (fp == fp0) + { + thisTriI = 0; + } + else if (fp == f.rcIndex(fp0)) + { + thisTriI = f.size()-3; + } + else + { + thisTriI = (fp-fp0-1) % (f.size()-2); + } + return thisTriI; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::tetDecomposer::tetDecomposer(const polyMesh& mesh) +: + mesh_(mesh) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::tetDecomposer::setRefinement +( + const decompositionType decomposeType, + polyTopoChange& meshMod +) +{ + cellToPoint_.setSize(mesh_.nCells()); + forAll(mesh_.cellCentres(), cellI) + { + // Any point on the cell + label masterPointI = mesh_.faces()[mesh_.cells()[cellI][0]][0]; + + cellToPoint_[cellI] = meshMod.addPoint + ( + mesh_.cellCentres()[cellI], + masterPointI, + -1, + true + ); + } + + + // Add face centre points + if (decomposeType == FACECENTRETETS) + { + faceToPoint_.setSize(mesh_.nFaces()); + forAll(mesh_.faceCentres(), faceI) + { + // Any point on the face + const label masterPointI = mesh_.faces()[faceI][0]; + + faceToPoint_[faceI] = meshMod.addPoint + ( + mesh_.faceCentres()[faceI], + masterPointI, + -1, + true + ); + } + } + + + // Per face, per point (faceCentre) or triangle (faceDiag) the added cell + faceOwnerCells_.setSize(mesh_.nFaces()); + faceNeighbourCells_.setSize(mesh_.nFaces()); + + if (decomposeType == FACECENTRETETS) + { + forAll(faceOwnerCells_, faceI) + { + const face& f = mesh_.faces()[faceI]; + faceOwnerCells_[faceI].setSize(f.size(), -1); + faceNeighbourCells_[faceI].setSize(f.size(), -1); + } + } + else + { + // Force construction of diagonal decomposition + (void)mesh_.tetBasePtIs(); + + forAll(faceOwnerCells_, faceI) + { + const face& f = mesh_.faces()[faceI]; + faceOwnerCells_[faceI].setSize(f.size()-2, -1); + faceNeighbourCells_[faceI].setSize(f.size()-2, -1); + } + } + + + forAll(mesh_.cells(), cellI) + { + const cell& cFaces = mesh_.cells()[cellI]; + + EdgeMap