diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.H b/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.H index c89d8dc6c9..7721cadaeb 100644 --- a/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.H +++ b/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.H @@ -29,6 +29,36 @@ Description Basic sub-grid obstacle flame-wrinking enhancement factor model. Details supplied by J Puttock 2/7/06. + Sub-grid flame area generation <\b> + + \f$ n = N - \hat{\dwea{\vec{U}}}.n_{s}.\hat{\dwea{\vec{U}}} \f$ + \f$ n_{r} = \sqrt{n} \f$ + + where: + + \f$ \hat{\dwea{\vec{U}}} = \dwea{\vec{U}} / \vert \dwea{\vec{U}} + \vert \f$ + + \f$ b = \hat{\dwea{\vec{U}}}.B.\hat{\dwea{\vec{U}}} / n_{r} \f$ + + where: + + \f$ B \f$ is the file "B". + + \f$ N \f$ is the file "N". + + \f$ n_{s} \f$ is the file "ns". + + The flame area enhancement factor \f$ \Xi_{sub} \f$ is expected to + approach: + + \f[ + \Xi_{{sub}_{eq}} = + 1 + max(2.2 \sqrt{b}, min(0.34 \frac{\vert \dwea{\vec{U}} + \vert}{{\vec{U}}^{'}}, 1.6)) \times min(\frac{n}{4}, 1) + \f] + + SourceFiles basicSubGrid.C diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.H b/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.H index 01a6d7cb11..a455fd28b4 100644 --- a/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.H +++ b/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.H @@ -25,10 +25,28 @@ License Class basicSubGrid + Description + Basic sub-grid obstacle flame-wrinking generation rate coefficient model. Details supplied by J Puttock 2/7/06. + \f$ G_{sub} \f$ denotes the generation coefficient and it is given by + + \f[ + G_{sub} = k_{1} /frac{\vert \dwea{\vec{U}} \vert}{L_{obs}} + \frac{/Xi_{{sub}_{eq}}-1}{/Xi_{sub}} + \f] + + and the removal: + + \f[ - k_{1} /frac{\vert \dwea{\vec{U}} \vert}{L_{sub}} + \frac{\Xi_{sub}-1}{\Xi_{sub}} \f] + + Finally, \f$ G_{sub} \f$ is added to generation rate \f$ G_{in} \f$ + due to the turbulence. + + SourceFiles basicSubGrid.C diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.H b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.H index ceb75ff827..d843214ab3 100644 --- a/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.H +++ b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.H @@ -29,6 +29,50 @@ Description Basic sub-grid obstacle drag model. Details supplied by J Puttock 2/7/06. + Sub-grid drag term <\b> + + The resistance term (force per unit of volume) is given by: + + \f[ + R = -\frac{1}{2} \rho \vert \dwea{\vec{U}} \vert \dwea{\vec{U}}.D + \f] + + where: + + \f$ D \f$ is the tensor field "CR" in \f$ m^{-1} \f$ + + This is term is treated implicitly in UEqn.H + + Sub-grid turbulence generation <\b> + + The turbulence source term \f$ G_{R} \f$ occurring in the + \f$ \kappa-\epsilon \f$ equations for the generation of turbulence due + to interaction with unresolved obstacles : + + \f$ G_{R} = C_{s}\beta_{\nu} + \mu_{eff} A_{w}^{2}(\dwea{\vec{U}}-\dwea{\vec{U}_{s}})^2 + \frac{1}{2} + \rho \vert \dwea{\vec{U}} \vert \dwea{\vec{U}}.T.\dwea{\vec{U}} \f$ + + where: + + \f$ C_{s} \f$ = 1 + + \f$ \beta_{\nu} \f$ is the volume porosity (file "betav"). + + \f$ \mu_{eff} \f$ is the effective viscosity. + + \f$ A_{w}^{2}\f$ is the obstacle surface area per unit of volume + (file "Aw"). + + \f$ \dwea{\vec{U}_{s}} \f$ is the slip velocity and is considered + \f$ \frac{1}{2}. \dwea{\vec{U}} \f$. + + \f$ T \f$ is a tensor in the file CT. + + The term \f$ G_{R} \f$ is treated explicitly in the \f$ \kappa-\epsilon + \f$ Eqs in the PDRkEpsilon.C file. + + SourceFiles basic.C @@ -40,7 +84,6 @@ SourceFiles #include "PDRDragModel.H" #include "XiEqModel.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H index 15d52e7269..b35ef9151e 100644 --- a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H +++ b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H @@ -26,7 +26,17 @@ Class PDRkEpsilon Description - Standard k-epsilon turbulence model. + Standard k-epsilon turbulence model with additional source terms + corresponding to PDR basic drag model (basic.H) + + The turbulence source term \f$ G_{R} \f$ appears in the + \f$ \kappa-\epsilon \f$ equation for the generation of turbulence due to + interaction with unresolved obstacles. + + In the \f$ \epsilon \f$ equation \f$ C_{1} G_{R} \f$ is added as a source + term. + + In the \f$ \kappa \f$ equation \f$ G_{R} \f$ is added as a source term. SourceFiles PDRkEpsilon.C diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModel.H b/applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModel.H index 13fa58a53e..2b76f9bc92 100644 --- a/applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModel.H +++ b/applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModel.H @@ -27,6 +27,57 @@ Class Description Base-class for all Xi models used by the b-Xi combustion model. + See Technical Report SH/RE/01R for details on the PDR modelling. + + Xi is given through an algebraic expression (algebraic.H), + by solving a transport equation (transport.H) or a fixed value (fixed.H). + See report TR/HGW/10 for details on the Weller two equations model. + + In the algebraic and transport methods \f$\Xi_{eq}\f$ is calculated in + similar way. In the algebraic approach, \f$\Xi_{eq}\f$ is the value used in + the \f$ b \f$ transport equation. + + \f$\Xi_{eq}\f$ is calculated as follows: + + \f$\Xi_{eq} = 1 + (1 + 2\Xi_{coeff}(0.5 - \dwea{b}))(\Xi^* - 1)\f$ + + where: + + \f$ \dwea{b} \f$ is the regress variable. + + \f$ \Xi^* \f$ is the total equilibrium wrinkling combining the effects + of the flame inestability and turbulence interaction and is given by + + \f[ + \Xi^* = \frac {R}{R - G_\eta - G_{in}} + \f] + + where: + + \f$ G_\eta \f$ is the generation rate of wrinkling due to turbulence + interaction. + + \f$ G_{in} = \kappa \rho_{u}/\rho_{b} \f$ is the generation + rate due to the flame inestability. + + By adding the removal rates of the two effects: + + \f[ + R = G_\eta \frac{\Xi_{\eta_{eq}}}{\Xi_{\eta_{eq}} - 1} + + G_{in} \frac{\Xi_{{in}_{eq}}}{\Xi_{{in}_{eq}} - 1} + \f] + + where: + + \f$ R \f$ is the total removal. + + \f$ G_\eta \f$ is a model constant. + + \f$ \Xi_{\eta_{eq}} \f$ is the flame wrinkling due to turbulence. + + \f$ \Xi_{{in}_{eq}} \f$ is the equilibrium level of the flame wrinkling + generated by inestability. It is a constant (default 2.5). + SourceFiles XiModel.C @@ -51,6 +102,8 @@ namespace Foam Class XiModel Declaration \*---------------------------------------------------------------------------*/ + + class XiModel { diff --git a/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H b/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H index e837fefa8b..8c65453b94 100644 --- a/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H +++ b/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H @@ -28,6 +28,33 @@ Class Description Laminar flame speed obtained from the SCOPE correlation. + Seven parameters are specified in terms of polynomial functions of + stoichiometry. Two polynomials are fitted, covering different parts of the + flammable range. If the mixture is outside the fitted range, linear + interpolation is used between the extreme of the polynomio and the upper or + lower flammable limit with the Markstein number constant. + + Variations of pressure and temperature from the reference values are taken + into account through \f$ pexp \f$ and \f$ texp \f$ + + The laminar burning velocity fitting polynomial is: + + \f$ Su = a_{0}(1+a_{1}x+K+..a_{i}x^{i}..+a_{6}x^{6}) (p/p_{ref})^{pexp} + (T/T_{ref})^{texp} \f$ + + where: + + \f$ a_{i} \f$ are the polinomial coefficients. + + \f$ pexp \f$ and \f$ texp \f$ are the pressure and temperature factors + respectively. + + \f$ x \f$ is the equivalence ratio. + + \f$ T_{ref} \f$ and \f$ p_{ref} \f$ are the temperature and pressure + references for the laminar burning velocity. + + SourceFiles SCOPELaminarFlameSpeed.C @@ -125,7 +152,7 @@ class SCOPE // corrected for temperature and pressure dependence inline scalar Su0pTphi(scalar p, scalar Tu, scalar phi) const; - //- Laminar flame speed evaluated from the given uniform + //- Laminar flame speed evaluated from the given uniform // equivalence ratio corrected for temperature and pressure dependence tmp Su0pTphi ( @@ -134,7 +161,7 @@ class SCOPE scalar phi ) const; - //- Laminar flame speed evaluated from the given equivalence ratio + //- Laminar flame speed evaluated from the given equivalence ratio // distribution corrected for temperature and pressure dependence tmp Su0pTphi ( @@ -144,7 +171,7 @@ class SCOPE ) const; //- Return the Markstein number - // evaluated from the given equivalence ratio + // evaluated from the given equivalence ratio tmp Ma(const volScalarField& phi) const; //- Construct as copy (not implemented) diff --git a/applications/test/dictionary/dictionaryTest.C b/applications/test/dictionary/dictionaryTest.C index 6064102781..37ae9bd90c 100644 --- a/applications/test/dictionary/dictionaryTest.C +++ b/applications/test/dictionary/dictionaryTest.C @@ -45,6 +45,27 @@ int main(int argc, char *argv[]) Info<< testDict << endl; + { + dictionary someDict; + someDict.add(keyType("a.*", true), "subdictValue"); + + dictionary dict; + dict.add("someDict", someDict); + dict.add(keyType(".*", true), "parentValue"); + + Info<< "dict:" << dict << endl; + + // Wildcard find. + Info<< "Wildcard find \"abc\" in top directory : " + << dict.lookup("abc") << endl; + Info<< "Wildcard find \"abc\" in sub directory : " + << dict.subDict("someDict").lookup("abc") + << endl; + Info<< "Recursive wildcard find \"def\" in sub directory : " + << dict.subDict("someDict").lookup("def", true) + << endl; + } + return 0; } 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