From ba10c424463d4a1ceb05017452af99a4e3199113 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 14 Nov 2008 14:50:51 +0000 Subject: [PATCH 01/19] isoSurface instead of isoSurfaceCell --- .../distanceSurface/distanceSurface.C | 170 +++++++++++------- .../distanceSurface/distanceSurface.H | 37 ++-- .../distanceSurfaceTemplates.C | 42 ++--- 3 files changed, 142 insertions(+), 107 deletions(-) diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C index 8c4acea996..f337d14253 100644 --- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C +++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C @@ -30,7 +30,8 @@ License #include "volPointInterpolation.H" #include "addToRunTimeSelectionTable.H" #include "fvMesh.H" -#include "isoSurfaceCell.H" +#include "isoSurface.H" +//#include "isoSurfaceCell.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -42,19 +43,46 @@ namespace Foam // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::distanceSurface::createGeometry() const +void Foam::distanceSurface::createGeometry() { // Clear any stored topo facesPtr_.clear(); + + const fvMesh& fvm = static_cast(mesh()); + // Distance to cell centres - scalarField cellDistance(mesh().nCells()); + // ~~~~~~~~~~~~~~~~~~~~~~~~ + + cellDistancePtr_.reset + ( + new volScalarField + ( + IOobject + ( + "cellDistance", + fvm.time().timeName(), + fvm.time(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + fvm, + dimensionedScalar("zero", dimless/dimTime, 0) + ) + ); + volScalarField& cellDistance = cellDistancePtr_(); + + // Internal field { + const pointField& cc = fvm.C(); + scalarField& fld = cellDistance.internalField(); + List nearest; surfPtr_().findNearest ( - mesh().cellCentres(), - scalarField(mesh().nCells(), GREAT), + cc, + scalarField(cc.size(), GREAT), nearest ); @@ -63,98 +91,109 @@ void Foam::distanceSurface::createGeometry() const vectorField normal; surfPtr_().getNormal(nearest, normal); - forAll(cellDistance, cellI) + forAll(nearest, i) { - vector d(mesh().cellCentres()[cellI]-nearest[cellI].hitPoint()); + vector d(cc[i]-nearest[i].hitPoint()); - if ((d&normal[cellI]) > 0) + if ((d&normal[i]) > 0) { - cellDistance[cellI] = Foam::mag(d); + fld[i] = Foam::mag(d); } else { - cellDistance[cellI] = -Foam::mag(d); + fld[i] = -Foam::mag(d); } } } else { - forAll(cellDistance, cellI) + forAll(nearest, i) { - cellDistance[cellI] = Foam::mag - ( - nearest[cellI].hitPoint() - - mesh().cellCentres()[cellI] - ); + fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint()); } } } + // Patch fields + { + forAll(fvm.C().boundaryField(), patchI) + { + const pointField& cc = fvm.C().boundaryField()[patchI]; + fvPatchScalarField& fld = cellDistance.boundaryField()[patchI]; + + List nearest; + surfPtr_().findNearest + ( + cc, + scalarField(cc.size(), GREAT), + nearest + ); + + if (signed_) + { + vectorField normal; + surfPtr_().getNormal(nearest, normal); + + forAll(nearest, i) + { + vector d(cc[i]-nearest[i].hitPoint()); + + if ((d&normal[i]) > 0) + { + fld[i] = Foam::mag(d); + } + else + { + fld[i] = -Foam::mag(d); + } + } + } + else + { + forAll(nearest, i) + { + fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint()); + } + } + } + } + + // On processor patches the mesh.C() will already be the cell centre + // on the opposite side so no need to swap cellDistance. + + // Distance to points - scalarField pointDistance(mesh().nPoints()); + pointDistance_.setSize(fvm.nPoints()); { List nearest; surfPtr_().findNearest ( - mesh().points(), - scalarField(mesh().nPoints(), GREAT), + fvm.points(), + scalarField(fvm.nPoints(), GREAT), nearest ); - forAll(pointDistance, pointI) + forAll(pointDistance_, pointI) { - pointDistance[pointI] = Foam::mag + pointDistance_[pointI] = Foam::mag ( nearest[pointI].hitPoint() - - mesh().points()[pointI] + - fvm.points()[pointI] ); } } //- Direct from cell field and point field. - const isoSurfaceCell iso + isoSurfPtr_.reset ( - mesh(), - cellDistance, - pointDistance, - distance_, - regularise_ + new isoSurface + ( + cellDistance, + pointDistance_, + distance_, + regularise_ + ) ); - ////- From point field and interpolated cell. - //scalarField cellAvg(mesh().nCells(), scalar(0.0)); - //labelField nPointCells(mesh().nCells(), 0); - //{ - // for (label pointI = 0; pointI < mesh().nPoints(); pointI++) - // { - // const labelList& pCells = mesh().pointCells(pointI); - // - // forAll(pCells, i) - // { - // label cellI = pCells[i]; - // - // cellAvg[cellI] += pointDistance[pointI]; - // nPointCells[cellI]++; - // } - // } - //} - //forAll(cellAvg, cellI) - //{ - // cellAvg[cellI] /= nPointCells[cellI]; - //} - // - //const isoSurface iso - //( - // mesh(), - // cellAvg, - // pointDistance, - // distance_, - // regularise_ - //); - - - const_cast(*this).triSurface::operator=(iso); - meshCells_ = iso.meshCells(); - if (debug) { print(Pout); @@ -193,9 +232,8 @@ Foam::distanceSurface::distanceSurface signed_(readBool(dict.lookup("signed"))), regularise_(dict.lookupOrDefault("regularise", true)), zoneName_(word::null), - facesPtr_(NULL), - storedTimeIndex_(-1), - meshCells_(0) + isoSurfPtr_(NULL), + facesPtr_(NULL) { // label zoneId = -1; // if (dict.readIfPresent("zone", zoneName_)) diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H index 82fba8a83f..64aebb1a9a 100644 --- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H +++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H @@ -37,8 +37,9 @@ SourceFiles #define distanceSurface_H #include "sampledSurface.H" -#include "triSurface.H" #include "searchableSurface.H" +//#include "isoSurfaceCell.H" +#include "isoSurface.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -51,8 +52,7 @@ namespace Foam class distanceSurface : - public sampledSurface, - public triSurface + public sampledSurface { // Private data @@ -71,23 +71,24 @@ class distanceSurface //- zone name (if restricted to zones) word zoneName_; + + //- Distance to cell centres + autoPtr cellDistancePtr_; + + //- Distance to points + scalarField pointDistance_; + + //- Constructed iso surface + autoPtr isoSurfPtr_; + //- triangles converted to faceList mutable autoPtr facesPtr_; - // Recreated for every isoSurface - - //- Time at last call - mutable label storedTimeIndex_; - - //- For every triangle the original cell in mesh - mutable labelList meshCells_; - - // Private Member Functions - //- Create iso surface (if time has changed) - void createGeometry() const; + //- Create iso surface + void createGeometry(); //- sample field on faces template @@ -129,7 +130,7 @@ public: //- Points of surface virtual const pointField& points() const { - return triSurface::points(); + return surface().points(); } //- Faces of surface @@ -137,7 +138,7 @@ public: { if (!facesPtr_.valid()) { - const triSurface& s = *this; + const triSurface& s = surface(); facesPtr_.reset(new faceList(s.size())); @@ -153,6 +154,10 @@ public: //- Correct for mesh movement and/or field changes virtual void correct(const bool meshChanged); + const isoSurface& surface() const + { + return isoSurfPtr_(); + } //- sample field on surface virtual tmp sample diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C index ce983c7f5f..fd1447af78 100644 --- a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C +++ b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C @@ -39,7 +39,7 @@ Foam::distanceSurface::sampleField const GeometricField& vField ) const { - return tmp >(new Field(vField, meshCells_)); + return tmp >(new Field(vField, surface().meshCells())); } @@ -50,33 +50,25 @@ Foam::distanceSurface::interpolateField const interpolation& interpolator ) const { - // One value per point - tmp > tvalues(new Field(points().size())); - Field& values = tvalues(); + const fvMesh& fvm = static_cast(mesh()); - boolList pointDone(points().size(), false); + // Get fields to sample. Assume volPointInterpolation! + const GeometricField& volFld = + interpolator.psi(); - forAll(faces(), cutFaceI) - { - const face& f = faces()[cutFaceI]; + tmp > pointFld + ( + volPointInterpolation::New(fvm).interpolate(volFld) + ); - forAll(f, faceVertI) - { - label pointI = f[faceVertI]; - - if (!pointDone[pointI]) - { - values[pointI] = interpolator.interpolate - ( - points()[pointI], - meshCells_[cutFaceI] - ); - pointDone[pointI] = true; - } - } - } - - return tvalues; + // Sample. + return surface().interpolate + ( + cellDistancePtr_(), + pointDistance_, + volFld, + pointFld() + ); } From ef19e4aab4545792630f1b25a59bc50693bc49e8 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 18 Nov 2008 12:35:43 +0000 Subject: [PATCH 02/19] better dualisation --- .../mesh/conversion/polyDualMesh/Make/files | 3 +- .../mesh/conversion/polyDualMesh/Make/options | 7 +- .../polyDualMesh/makePolyDualMesh.C | 515 ++++++ .../conversion/polyDualMesh/meshDualiser.C | 1489 +++++++++++++++++ .../conversion/polyDualMesh/meshDualiser.H | 262 +++ .../conversion/polyDualMesh/polyDualMeshApp.C | 79 - 6 files changed, 2273 insertions(+), 82 deletions(-) create mode 100644 applications/utilities/mesh/conversion/polyDualMesh/makePolyDualMesh.C create mode 100644 applications/utilities/mesh/conversion/polyDualMesh/meshDualiser.C create mode 100644 applications/utilities/mesh/conversion/polyDualMesh/meshDualiser.H delete mode 100644 applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C diff --git a/applications/utilities/mesh/conversion/polyDualMesh/Make/files b/applications/utilities/mesh/conversion/polyDualMesh/Make/files index c43f79e8e1..752da5cfdd 100644 --- a/applications/utilities/mesh/conversion/polyDualMesh/Make/files +++ b/applications/utilities/mesh/conversion/polyDualMesh/Make/files @@ -1,3 +1,4 @@ -polyDualMeshApp.C +meshDualiser.C +makePolyDualMesh.C EXE = $(FOAM_APPBIN)/polyDualMesh diff --git a/applications/utilities/mesh/conversion/polyDualMesh/Make/options b/applications/utilities/mesh/conversion/polyDualMesh/Make/options index 6dc63a7a98..dc318df998 100644 --- a/applications/utilities/mesh/conversion/polyDualMesh/Make/options +++ b/applications/utilities/mesh/conversion/polyDualMesh/Make/options @@ -1,6 +1,9 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/conversion/lnInclude + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude EXE_LIBS = \ - -lmeshTools -lconversion + -lfiniteVolume \ + -ldynamicMesh \ + -lmeshTools diff --git a/applications/utilities/mesh/conversion/polyDualMesh/makePolyDualMesh.C b/applications/utilities/mesh/conversion/polyDualMesh/makePolyDualMesh.C new file mode 100644 index 0000000000..26eb18e26d --- /dev/null +++ b/applications/utilities/mesh/conversion/polyDualMesh/makePolyDualMesh.C @@ -0,0 +1,515 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Description + Calculate the dual of a polyMesh. Adheres to all the feature&patch edges. + + Feature angle: + convex features : point becomes single boundary cell with multiple + boundary faces. + concave features: point becomes multiple boundary cells. + + -splitAllFaces: + Normally only constructs a single face between two cells. This single face + might be too distorted. splitAllFaces will create a single face for every + original cell the face passes through. The mesh will thus have + multiple faces inbetween two cells! (so is not strictly upper-triangular + anymore - checkMesh will complain) + -doNotPreserveFaceZones: + By default all faceZones are preserved by marking all faces, edges and + points on them as features. The -doNotPreserveFaceZones disables this + behaviour. + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "Time.H" +#include "fvMesh.H" +#include "mathematicalConstants.H" +#include "polyTopoChange.H" +#include "mapPolyMesh.H" +#include "PackedList.H" +#include "meshTools.H" +#include "OFstream.H" +#include "meshDualiser.H" +#include "ReadFields.H" +#include "volFields.H" +#include "surfaceFields.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Naive feature detection. All boundary edges with angle > featureAngle become +// feature edges. All points on feature edges become feature points. All +// boundary faces become feature faces. +void simpleMarkFeatures +( + const polyMesh& mesh, + const PackedList<1>& isBoundaryEdge, + const scalar featureAngle, + const bool doNotPreserveFaceZones, + + labelList& featureFaces, + labelList& featureEdges, + labelList& singleCellFeaturePoints, + labelList& multiCellFeaturePoints +) +{ + scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0); + + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + + // Working sets + labelHashSet featureEdgeSet; + labelHashSet singleCellFeaturePointSet; + labelHashSet multiCellFeaturePointSet; + + + // 1. Mark all edges between patches + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + const labelList& meshEdges = pp.meshEdges(); + + // All patch corner edges. These need to be feature points & edges! + for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++) + { + label meshEdgeI = meshEdges[edgeI]; + featureEdgeSet.insert(meshEdgeI); + singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][0]); + singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][1]); + } + } + + + + // 2. Mark all geometric feature edges + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Make distinction between convex features where the boundary point becomes + // a single cell and concave features where the boundary point becomes + // multiple 'half' cells. + + // Addressing for all outside faces + primitivePatch allBoundary + ( + SubList + ( + mesh.faces(), + mesh.nFaces()-mesh.nInternalFaces(), + mesh.nInternalFaces() + ), + mesh.points() + ); + + // Check for non-manifold points (surface pinched at point) + allBoundary.checkPointManifold(false, &singleCellFeaturePointSet); + + // Check for non-manifold edges (surface pinched at edge) + const labelListList& edgeFaces = allBoundary.edgeFaces(); + const labelList& meshPoints = allBoundary.meshPoints(); + + forAll(edgeFaces, edgeI) + { + const labelList& eFaces = edgeFaces[edgeI]; + + if (eFaces.size() > 2) + { + const edge& e = allBoundary.edges()[edgeI]; + + //Info<< "Detected non-manifold boundary edge:" << edgeI + // << " coords:" + // << allBoundary.points()[meshPoints[e[0]]] + // << allBoundary.points()[meshPoints[e[1]]] << endl; + + singleCellFeaturePointSet.insert(meshPoints[e[0]]); + singleCellFeaturePointSet.insert(meshPoints[e[1]]); + } + } + + // Check for features. + forAll(edgeFaces, edgeI) + { + const labelList& eFaces = edgeFaces[edgeI]; + + if (eFaces.size() == 2) + { + label f0 = eFaces[0]; + label f1 = eFaces[1]; + + // check angle + const vector& n0 = allBoundary.faceNormals()[f0]; + const vector& n1 = allBoundary.faceNormals()[f1]; + + if ((n0 & n1) < minCos) + { + const edge& e = allBoundary.edges()[edgeI]; + label v0 = meshPoints[e[0]]; + label v1 = meshPoints[e[1]]; + + label meshEdgeI = meshTools::findEdge(mesh, v0, v1); + featureEdgeSet.insert(meshEdgeI); + + // Check if convex or concave by looking at angle + // between face centres and normal + vector c1c0 + ( + allBoundary[f1].centre(allBoundary.points()) + - allBoundary[f0].centre(allBoundary.points()) + ); + + if ((c1c0 & n0) > SMALL) + { + // Found concave edge. Make into multiCell features + Info<< "Detected concave feature edge:" << edgeI + << " cos:" << (c1c0 & n0) + << " coords:" + << allBoundary.points()[v0] + << allBoundary.points()[v1] + << endl; + + singleCellFeaturePointSet.erase(v0); + multiCellFeaturePointSet.insert(v0); + singleCellFeaturePointSet.erase(v1); + multiCellFeaturePointSet.insert(v1); + } + else + { + // Convex. singleCell feature. + if (!multiCellFeaturePointSet.found(v0)) + { + singleCellFeaturePointSet.insert(v0); + } + if (!multiCellFeaturePointSet.found(v1)) + { + singleCellFeaturePointSet.insert(v1); + } + } + } + } + } + + + // 3. Mark all feature faces + // ~~~~~~~~~~~~~~~~~~~~~~~~~ + + // Face centres that need inclusion in the dual mesh + labelHashSet featureFaceSet(mesh.nFaces()-mesh.nInternalFaces()); + // A. boundary faces. + for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++) + { + featureFaceSet.insert(faceI); + } + + // B. face zones. + const faceZoneMesh& faceZones = mesh.faceZones(); + + if (doNotPreserveFaceZones) + { + if (faceZones.size() > 0) + { + WarningIn("simpleMarkFeatures(..)") + << "Detected " << faceZones.size() + << " faceZones. These will not be preserved." + << endl; + } + } + else + { + if (faceZones.size() > 0) + { + Info<< "Detected " << faceZones.size() + << " faceZones. Preserving these by marking their" + << " points, edges and faces as features." << endl; + } + + forAll(faceZones, zoneI) + { + const faceZone& fz = faceZones[zoneI]; + + Info<< "Inserting all faces in faceZone " << fz.name() + << " as features." << endl; + + forAll(fz, i) + { + label faceI = fz[i]; + const face& f = mesh.faces()[faceI]; + const labelList& fEdges = mesh.faceEdges()[faceI]; + + featureFaceSet.insert(faceI); + forAll(f, fp) + { + // Mark point as multi cell point (since both sides of + // face should have different cells) + singleCellFeaturePointSet.erase(f[fp]); + multiCellFeaturePointSet.insert(f[fp]); + + // Make sure there are points on the edges. + featureEdgeSet.insert(fEdges[fp]); + } + } + } + } + + // Transfer to arguments + featureFaces = featureFaceSet.toc(); + featureEdges = featureEdgeSet.toc(); + singleCellFeaturePoints = singleCellFeaturePointSet.toc(); + multiCellFeaturePoints = multiCellFeaturePointSet.toc(); +} + + +// Dump features to .obj files +void dumpFeatures +( + const polyMesh& mesh, + const labelList& featureFaces, + const labelList& featureEdges, + const labelList& singleCellFeaturePoints, + const labelList& multiCellFeaturePoints +) +{ + { + OFstream str("featureFaces.obj"); + Info<< "Dumping centres of featureFaces to obj file " << str.name() + << endl; + forAll(featureFaces, i) + { + meshTools::writeOBJ(str, mesh.faceCentres()[featureFaces[i]]); + } + } + { + OFstream str("featureEdges.obj"); + Info<< "Dumping featureEdges to obj file " << str.name() << endl; + label vertI = 0; + + forAll(featureEdges, i) + { + const edge& e = mesh.edges()[featureEdges[i]]; + meshTools::writeOBJ(str, mesh.points()[e[0]]); + vertI++; + meshTools::writeOBJ(str, mesh.points()[e[1]]); + vertI++; + str<< "l " << vertI-1 << ' ' << vertI << nl; + } + } + { + OFstream str("singleCellFeaturePoints.obj"); + Info<< "Dumping featurePoints that become a single cell to obj file " + << str.name() << endl; + forAll(singleCellFeaturePoints, i) + { + meshTools::writeOBJ(str, mesh.points()[singleCellFeaturePoints[i]]); + } + } + { + OFstream str("multiCellFeaturePoints.obj"); + Info<< "Dumping featurePoints that become multiple cells to obj file " + << str.name() << endl; + forAll(multiCellFeaturePoints, i) + { + meshTools::writeOBJ(str, mesh.points()[multiCellFeaturePoints[i]]); + } + } +} + + +int main(int argc, char *argv[]) +{ + argList::noParallel(); +# include "addTimeOptions.H" + argList::validArgs.append("feature angle[0-180]"); + argList::validOptions.insert("splitAllFaces", ""); + argList::validOptions.insert("doNotPreserveFaceZones", ""); + argList::validOptions.insert("overwrite", ""); + +# include "setRootCase.H" +# include "createTime.H" + // Get times list + instantList Times = runTime.times(); +# include "checkTimeOptions.H" + runTime.setTime(Times[startTime], startTime); + +# include "createMesh.H" + + // Mark boundary edges and points. + // (Note: in 1.4.2 we can use the built-in mesh point ordering + // facility instead) + PackedList<1> isBoundaryEdge(mesh.nEdges()); + for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++) + { + const labelList& fEdges = mesh.faceEdges()[faceI]; + + forAll(fEdges, i) + { + isBoundaryEdge.set(fEdges[i], 1); + } + } + + scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])())); + + scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0); + + Info<< "Feature:" << featureAngle << endl + << "minCos :" << minCos << endl + << endl; + + + const bool splitAllFaces = args.options().found("splitAllFaces"); + const bool overwrite = args.options().found("overwrite"); + const bool doNotPreserveFaceZones = args.options().found + ( + "doNotPreserveFaceZones" + ); + + // Face(centre)s that need inclusion in the dual mesh + labelList featureFaces; + // Edge(centre)s ,, + labelList featureEdges; + // Points (that become a single cell) that need inclusion in the dual mesh + labelList singleCellFeaturePoints; + // Points (that become a mulitple cells) ,, + labelList multiCellFeaturePoints; + + // Sample implementation of feature detection. + simpleMarkFeatures + ( + mesh, + isBoundaryEdge, + featureAngle, + doNotPreserveFaceZones, + + featureFaces, + featureEdges, + singleCellFeaturePoints, + multiCellFeaturePoints + ); + + // If we want to split all polyMesh faces into one dualface per cell + // we are passing through we also need a point + // at the polyMesh facecentre and edgemid of the faces we want to + // split. + if (splitAllFaces) + { + featureEdges = identity(mesh.nEdges()); + featureFaces = identity(mesh.nFaces()); + } + + // Write obj files for debugging + dumpFeatures + ( + mesh, + featureFaces, + featureEdges, + singleCellFeaturePoints, + multiCellFeaturePoints + ); + + + + // Read objects in time directory + IOobjectList objects(mesh, runTime.timeName()); + + // Read vol fields. + PtrList vsFlds; + ReadFields(mesh, objects, vsFlds); + + PtrList vvFlds; + ReadFields(mesh, objects, vvFlds); + + PtrList vstFlds; + ReadFields(mesh, objects, vstFlds); + + PtrList vsymtFlds; + ReadFields(mesh, objects, vsymtFlds); + + PtrList vtFlds; + ReadFields(mesh, objects, vtFlds); + + // Read surface fields. + PtrList ssFlds; + ReadFields(mesh, objects, ssFlds); + + PtrList svFlds; + ReadFields(mesh, objects, svFlds); + + PtrList sstFlds; + ReadFields(mesh, objects, sstFlds); + + PtrList ssymtFlds; + ReadFields(mesh, objects, ssymtFlds); + + PtrList stFlds; + ReadFields(mesh, objects, stFlds); + + + // Topo change container + polyTopoChange meshMod(mesh.boundaryMesh().size()); + + // Mesh dualiser engine + meshDualiser dualMaker(mesh); + + // Insert all commands into polyTopoChange to create dual of mesh. This does + // all the hard work. + dualMaker.setRefinement + ( + splitAllFaces, + featureFaces, + featureEdges, + singleCellFeaturePoints, + multiCellFeaturePoints, + meshMod + ); + + // Create mesh, return map from old to new mesh. + autoPtr map = meshMod.changeMesh(mesh, false); + + // Update fields + mesh.updateMesh(map); + + // Optionally inflate mesh + if (map().hasMotionPoints()) + { + mesh.movePoints(map().preMotionPoints()); + } + + if (!overwrite) + { + runTime++; + mesh.setInstance(runTime.timeName()); + } + + Info<< "Writing dual mesh to " << runTime.timeName() << endl; + + mesh.write(); + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/conversion/polyDualMesh/meshDualiser.C b/applications/utilities/mesh/conversion/polyDualMesh/meshDualiser.C new file mode 100644 index 0000000000..4492a6f6b6 --- /dev/null +++ b/applications/utilities/mesh/conversion/polyDualMesh/meshDualiser.C @@ -0,0 +1,1489 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + meshDualiser + +\*---------------------------------------------------------------------------*/ + +#include "meshDualiser.H" +#include "meshTools.H" +#include "polyMesh.H" +#include "polyTopoChange.H" +#include "mapPolyMesh.H" +#include "edgeFaceCirculator.H" +#include "mergePoints.H" +#include "OFstream.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(meshDualiser, 0); +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::meshDualiser::checkPolyTopoChange(const polyTopoChange& meshMod) +{ + // Assume no removed points + pointField points(meshMod.points().size()); + forAll(meshMod.points(), i) + { + points[i] = meshMod.points()[i]; + } + + labelList oldToNew; + pointField newPoints; + bool hasMerged = mergePoints + ( + points, + 1E-6, + false, + oldToNew, + newPoints + ); + + if (hasMerged) + { + labelListList newToOld(invertOneToMany(newPoints.size(), oldToNew)); + + forAll(newToOld, newI) + { + if (newToOld[newI].size() != 1) + { + FatalErrorIn + ( + "meshDualiser::checkPolyTopoChange(const polyTopoChange&)" + ) << "duplicate verts:" << newToOld[newI] + << " coords:" + << IndirectList(points, newToOld[newI])() + << abort(FatalError); + } + } + } +} + + +// Dump state so far. +void Foam::meshDualiser::dumpPolyTopoChange +( + const polyTopoChange& meshMod, + const fileName& prefix +) +{ + OFstream str1(prefix + "Faces.obj"); + OFstream str2(prefix + "Edges.obj"); + + Info<< "Dumping current polyTopoChange. Faces to " << str1.name() + << " , points and edges to " << str2.name() << endl; + + const DynamicList& points = meshMod.points(); + + forAll(points, pointI) + { + meshTools::writeOBJ(str1, points[pointI]); + meshTools::writeOBJ(str2, points[pointI]); + } + + const DynamicList& faces = meshMod.faces(); + + forAll(faces, faceI) + { + const face& f = faces[faceI]; + + str1<< 'f'; + forAll(f, fp) + { + str1<< ' ' << f[fp]+1; + } + str1<< nl; + + str2<< 'l'; + forAll(f, fp) + { + str2<< ' ' << f[fp]+1; + } + str2<< ' ' << f[0]+1 << nl; + } +} + + +//- Given cell and point on mesh finds the corresponding dualCell. Most +// points become only one cell but the feature points might be split. +Foam::label Foam::meshDualiser::findDualCell +( + const label cellI, + const label pointI +) const +{ + const labelList& dualCells = pointToDualCells_[pointI]; + + if (dualCells.size() == 1) + { + return dualCells[0]; + } + else + { + label index = findIndex(mesh_.pointCells()[pointI], cellI); + + return dualCells[index]; + } +} + + +// Helper function to generate dualpoints on all boundary edges emanating +// from (boundary & feature) point +void Foam::meshDualiser::generateDualBoundaryEdges +( + const PackedList<1>& isBoundaryEdge, + const label pointI, + polyTopoChange& meshMod +) +{ + const labelList& pEdges = mesh_.pointEdges()[pointI]; + + forAll(pEdges, pEdgeI) + { + label edgeI = pEdges[pEdgeI]; + + if (edgeToDualPoint_[edgeI] == -1 && isBoundaryEdge.get(edgeI) == 1) + { + const edge& e = mesh_.edges()[edgeI]; + + edgeToDualPoint_[edgeI] = meshMod.addPoint + ( + e.centre(mesh_.points()), + pointI, // masterPoint + -1, // zoneID + true // inCell + ); + } + } +} + + +// Return true if point on face has same dual cells on both owner and neighbour +// sides. +bool Foam::meshDualiser::sameDualCell +( + const label faceI, + const label pointI +) const +{ + if (!mesh_.isInternalFace(faceI)) + { + FatalErrorIn("sameDualCell(const label, const label)") + << "face:" << faceI << " is not internal face." + << abort(FatalError); + } + + label own = mesh_.faceOwner()[faceI]; + label nei = mesh_.faceNeighbour()[faceI]; + + return findDualCell(own, pointI) == findDualCell(nei, pointI); +} + + +Foam::label Foam::meshDualiser::addInternalFace +( + const label masterPointI, + const label masterEdgeI, + const label masterFaceI, + + const bool edgeOrder, + const label dualCell0, + const label dualCell1, + const DynamicList