diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.H b/applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.H
index c89d8dc6c9..04cc2a83ca 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
+
+ \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..05ec433919 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
+
+ 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
+
+ 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/DynamicList/DynamicListTest.C b/applications/test/DynamicList/DynamicListTest.C
index 6cdbb84890..28960444a2 100644
--- a/applications/test/DynamicList/DynamicListTest.C
+++ b/applications/test/DynamicList/DynamicListTest.C
@@ -43,16 +43,21 @@ int main(int argc, char *argv[])
ldl[0](3) = 3;
ldl[0](1) = 1;
- ldl[0].setSize(5); // increase allocated size
- ldl[1].setSize(10); // increase allocated size
- ldl[1](2) = 2;
+ ldl[0].setCapacity(5); // increase allocated size
+ ldl[1].setCapacity(10); // increase allocated size
+ ldl[0].reserve(15); // should increase allocated size
+ ldl[1].reserve(5); // should not decrease allocated size
+ ldl[1](3) = 2; // allocates space and sets value
+
+ // this works without a segfault, but doesn't change the list size
+ ldl[0][4] = 4;
ldl[1] = 3;
Info<< "" << ldl << " " << nl << "sizes: ";
forAll(ldl, i)
{
- Info<< " " << ldl[i].size() << "/" << ldl[i].allocSize();
+ Info<< " " << ldl[i].size() << "/" << ldl[i].capacity();
}
Info<< endl;
@@ -63,7 +68,7 @@ int main(int argc, char *argv[])
Info<< "" << ldl << " " << nl << "sizes: ";
forAll(ldl, i)
{
- Info<< " " << ldl[i].size() << "/" << ldl[i].allocSize();
+ Info<< " " << ldl[i].size() << "/" << ldl[i].capacity();
}
Info<< endl;
@@ -78,10 +83,10 @@ int main(int argc, char *argv[])
{
dlA.append(i);
}
- dlA.setSize(10);
+ dlA.setCapacity(10);
Info<< "" << dlA << " " << nl << "sizes: "
- << " " << dlA.size() << "/" << dlA.allocSize() << endl;
+ << " " << dlA.size() << "/" << dlA.capacity() << endl;
dlB.transfer(dlA);
@@ -91,10 +96,54 @@ int main(int argc, char *argv[])
Info<< "Transferred to dlB" << endl;
Info<< "" << dlA << " " << nl << "sizes: "
- << " " << dlA.size() << "/" << dlA.allocSize() << endl;
+ << " " << dlA.size() << "/" << dlA.capacity() << endl;
Info<< "" << dlB << " " << nl << "sizes: "
- << " " << dlB.size() << "/" << dlB.allocSize() << endl;
+ << " " << dlB.size() << "/" << dlB.capacity() << endl;
+ // try with a normal list:
+ List lstA;
+ lstA.transfer(dlB);
+ Info<< "Transferred to normal list" << endl;
+ Info<< "" << lstA << " " << nl << "sizes: "
+ << " " << lstA.size() << endl;
+ Info<< "" << dlB << " " << nl << "sizes: "
+ << " " << dlB.size() << "/" << dlB.capacity() << endl;
+
+ // Copy back and append a few time
+ for (label i=0; i < 3; i++)
+ {
+ dlB.append(lstA);
+ }
+
+ Info<< "appended list a few times" << endl;
+ Info<< "" << dlB << " " << nl << "sizes: "
+ << " " << dlB.size() << "/" << dlB.capacity() << endl;
+
+ // assign the list (should maintain allocated space)
+ dlB = lstA;
+ Info<< "assigned list" << endl;
+ Info<< "" << dlB << " " << nl << "sizes: "
+ << " " << dlB.size() << "/" << dlB.capacity() << endl;
+
+
+ // Copy back and append a few time
+ for (label i=0; i < 3; i++)
+ {
+ dlB.append(lstA);
+ }
+
+
+ // check allocation granularity
+ DynamicList dlC;
+
+ Info<< "" << dlC << " " << nl << "sizes: "
+ << " " << dlC.size() << "/" << dlC.capacity() << endl;
+
+ dlC.reserve(dlB.size());
+ dlC = dlB;
+
+ Info<< "" << dlC << " " << nl << "sizes: "
+ << " " << dlC.size() << "/" << dlC.capacity() << endl;
return 0;
}
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/test/sort/sortListTest.C b/applications/test/sort/sortListTest.C
index 89252ae1f3..06d5b01186 100644
--- a/applications/test/sort/sortListTest.C
+++ b/applications/test/sort/sortListTest.C
@@ -27,6 +27,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "SortableList.H"
+#include "ListOps.H"
using namespace Foam;
@@ -35,42 +36,101 @@ using namespace Foam;
int main(int argc, char *argv[])
{
- labelList orig(5);
+ labelList orig(8);
orig[0] = 7;
- orig[1] = 4;
+ orig[1] = 9;
orig[2] = 1;
orig[3] = 2;
- orig[4] = 9;
+ orig[4] = 4;
+ orig[5] = 7;
+ orig[6] = 4;
+ orig[7] = 0;
+
+ labelList order;
labelList a(orig);
- Info << "before: " << a << endl;
+ sortedOrder(a, order);
+
+ Info<< "unsorted: " << a << endl;
sort(a);
- Info << "after: " << a << endl;
+ Info<< "sorted: " << a << endl;
+ Info<< "indices: " << order << endl;
SortableList b(orig);
- Info << "sorted: " << b << endl;
- Info << "indices: " << b.indices() << endl;
+ Info<< "unsorted: " << orig << endl;
+ Info<< "sorted: " << b << endl;
+ Info<< "indices: " << b.indices() << endl;
- Info << "shrunk: " << b.shrink() << endl;
- Info << "indices: " << b.indices() << endl;
+ Info<< "shrunk: " << b.shrink() << endl;
+ Info<< "indices: " << b.indices() << endl;
// repeat by assignment
b = orig;
- Info << "unsorted: " << b << endl;
+ Info<< "unsorted: " << b << endl;
b.sort();
- Info << "sorted: " << b << endl;
- Info << "indices: " << b.indices() << endl;
+ Info<< "sorted: " << b << endl;
+ Info<< "indices: " << b.indices() << endl;
+
+ // find unique/duplicate values
+ b = orig;
+
+ Info<< "unsorted: " << b << endl;
+ uniqueOrder(b, order);
+ Info<< "unique: " << order << endl;
+ duplicateOrder(b, order);
+ Info<< "duplicate:" << order << endl;
+
+ // sort reverse
+ Info<< "unsorted: " << b << endl;
+ b.reverseSort();
+ Info<< "rsort: " << b << endl;
+ Info<< "indices: " << b.indices() << endl;
// transfer assignment
- b.transfer(orig);
- Info << "unsorted: " << b << endl;
+ a = orig;
+ b.transfer(a);
+ Info<< "unsorted: " << b << endl;
b.sort();
- Info << "sorted: " << b << endl;
- Info << "indices: " << b.indices() << endl;
+ Info<< "sorted: " << b << endl;
+ Info<< "indices: " << b.indices() << endl;
- Info << "End\n" << endl;
+ a.transfer(b);
+
+ Info<< "plain: " << a << endl;
+ Info<< "sorted: " << b << endl;
+ Info<< "indices: " << b.indices() << endl;
+
+ // sort/duplicate/unique with identical values
+ b.setSize(8);
+ b = 5;
+
+ Info<< "unsorted: " << b << endl;
+
+ uniqueOrder(b, order);
+ Info<< "unique: " << order << endl;
+ duplicateOrder(b, order);
+ Info<< "duplicate:" << order << endl;
+ b.sort();
+
+ Info<< "sorted: " << b << endl;
+ Info<< "indices: " << b.indices() << endl;
+
+ // with a single value
+ b.setSize(1);
+
+ Info<< "unsorted: " << b << endl;
+ uniqueOrder(b, order);
+ Info<< "unique: " << order << endl;
+ duplicateOrder(b, order);
+ Info<< "duplicate:" << order << endl;
+ b.sort();
+
+ Info<< "sorted: " << b << endl;
+ Info<< "indices: " << b.indices() << endl;
+
+ Info<< "\nEnd\n" << endl;
return 0;
}
diff --git a/applications/test/xfer/xferListTest.C b/applications/test/xfer/xferListTest.C
index 6f2af0932a..25d0d9f3db 100644
--- a/applications/test/xfer/xferListTest.C
+++ b/applications/test/xfer/xferListTest.C
@@ -108,13 +108,13 @@ int main(int argc, char *argv[])
face f1(dl);
face f2(xferCopy(dl));
- Info<< "dl[" << dl.size() << "/" << dl.allocSize() << "] " << dl << endl;
+ Info<< "dl[" << dl.size() << "/" << dl.capacity() << "] " << dl << endl;
Info<< "f1: " << f1 << endl;
Info<< "f2: " << f2 << endl;
// note: using xferMoveTo to ensure the correct transfer() method is called
face f3( xferMoveTo(dl) );
- Info<< "dl[" << dl.size() << "/" << dl.allocSize() << "] " << dl << endl;
+ Info<< "dl[" << dl.size() << "/" << dl.capacity() << "] " << dl << endl;
Info<< "f3: " << f3 << 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& verts,
+ polyTopoChange& meshMod
+) const
+{
+ face newFace(verts);
+
+ if (edgeOrder != (dualCell0 < dualCell1))
+ {
+ reverse(newFace);
+ }
+
+ if (debug)
+ {
+ pointField facePoints
+ (
+ IndirectList(meshMod.points(), newFace)()
+ );
+
+ labelList oldToNew;
+ pointField newPoints;
+ bool hasMerged = mergePoints
+ (
+ facePoints,
+ 1E-6,
+ false,
+ oldToNew,
+ newPoints
+ );
+
+ if (hasMerged)
+ {
+ FatalErrorIn("addInternalFace(..)")
+ << "verts:" << verts << " newFace:" << newFace
+ << " face points:" << facePoints
+ << abort(FatalError);
+ }
+ }
+
+
+ label zoneID = -1;
+ bool zoneFlip = false;
+ if (masterFaceI != -1)
+ {
+ zoneID = mesh_.faceZones().whichZone(masterFaceI);
+
+ if (zoneID != -1)
+ {
+ const faceZone& fZone = mesh_.faceZones()[zoneID];
+
+ zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
+ }
+ }
+
+ label dualFaceI;
+
+ if (dualCell0 < dualCell1)
+ {
+ dualFaceI = meshMod.addFace
+ (
+ newFace,
+ dualCell0, // own
+ dualCell1, // nei
+ masterPointI, // masterPointID
+ masterEdgeI, // masterEdgeID
+ masterFaceI, // masterFaceID
+ false, // flipFaceFlux
+ -1, // patchID
+ zoneID, // zoneID
+ zoneFlip // zoneFlip
+ );
+
+ //pointField dualPoints(meshMod.points());
+ //vector n(newFace.normal(dualPoints));
+ //n /= mag(n);
+ //Pout<< "Generated internal dualFace:" << dualFaceI
+ // << " verts:" << newFace
+ // << " points:" << IndirectList(meshMod.points(), newFace)()
+ // << " n:" << n
+ // << " between dualowner:" << dualCell0
+ // << " dualneigbour:" << dualCell1
+ // << endl;
+ }
+ else
+ {
+ dualFaceI = meshMod.addFace
+ (
+ newFace,
+ dualCell1, // own
+ dualCell0, // nei
+ masterPointI, // masterPointID
+ masterEdgeI, // masterEdgeID
+ masterFaceI, // masterFaceID
+ false, // flipFaceFlux
+ -1, // patchID
+ zoneID, // zoneID
+ zoneFlip // zoneFlip
+ );
+
+ //pointField dualPoints(meshMod.points());
+ //vector n(newFace.normal(dualPoints));
+ //n /= mag(n);
+ //Pout<< "Generated internal dualFace:" << dualFaceI
+ // << " verts:" << newFace
+ // << " points:" << IndirectList(meshMod.points(), newFace)()
+ // << " n:" << n
+ // << " between dualowner:" << dualCell1
+ // << " dualneigbour:" << dualCell0
+ // << endl;
+ }
+ return dualFaceI;
+}
+
+
+Foam::label Foam::meshDualiser::addBoundaryFace
+(
+ const label masterPointI,
+ const label masterEdgeI,
+ const label masterFaceI,
+
+ const label dualCellI,
+ const label patchI,
+ const DynamicList& verts,
+ polyTopoChange& meshMod
+) const
+{
+ face newFace(verts);
+
+ label zoneID = -1;
+ bool zoneFlip = false;
+ if (masterFaceI != -1)
+ {
+ zoneID = mesh_.faceZones().whichZone(masterFaceI);
+
+ if (zoneID != -1)
+ {
+ const faceZone& fZone = mesh_.faceZones()[zoneID];
+
+ zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
+ }
+ }
+
+ label dualFaceI = meshMod.addFace
+ (
+ newFace,
+ dualCellI, // own
+ -1, // nei
+ masterPointI, // masterPointID
+ masterEdgeI, // masterEdgeID
+ masterFaceI, // masterFaceID
+ false, // flipFaceFlux
+ patchI, // patchID
+ zoneID, // zoneID
+ zoneFlip // zoneFlip
+ );
+
+ //pointField dualPoints(meshMod.points());
+ //vector n(newFace.normal(dualPoints));
+ //n /= mag(n);
+ //Pout<< "Generated boundary dualFace:" << dualFaceI
+ // << " verts:" << newFace
+ // << " points:" << IndirectList(meshMod.points(), newFace)()
+ // << " n:" << n
+ // << " on dualowner:" << dualCellI
+ // << endl;
+ return dualFaceI;
+}
+
+
+// Walks around edgeI.
+// splitFace=true : creates multiple faces
+// splitFace=false: creates single face if same dual cells on both sides,
+// multiple faces otherwise.
+void Foam::meshDualiser::createFacesAroundEdge
+(
+ const bool splitFace,
+ const PackedList<1>& isBoundaryEdge,
+ const label edgeI,
+ const label startFaceI,
+ polyTopoChange& meshMod,
+ boolList& doneEFaces
+) const
+{
+ const edge& e = mesh_.edges()[edgeI];
+ const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+
+ label fp = edgeFaceCirculator::getMinIndex
+ (
+ mesh_.faces()[startFaceI],
+ e[0],
+ e[1]
+ );
+
+ edgeFaceCirculator ie
+ (
+ mesh_,
+ startFaceI, // face
+ true, // ownerSide
+ fp, // fp
+ isBoundaryEdge.get(edgeI) == 1 // isBoundaryEdge
+ );
+ ie.setCanonical();
+
+ bool edgeOrder = ie.sameOrder(e[0], e[1]);
+ label startFaceLabel = ie.faceLabel();
+
+ //Pout<< "At edge:" << edgeI << " verts:" << e
+ // << " points:" << mesh_.points()[e[0]] << mesh_.points()[e[1]]
+ // << " started walking at face:" << ie.faceLabel()
+ // << " verts:" << mesh_.faces()[ie.faceLabel()]
+ // << " edgeOrder:" << edgeOrder
+ // << " in direction of cell:" << ie.cellLabel()
+ // << endl;
+
+ // Walk and collect face.
+ DynamicList verts(100);
+
+ if (edgeToDualPoint_[edgeI] != -1)
+ {
+ verts.append(edgeToDualPoint_[edgeI]);
+ }
+ if (faceToDualPoint_[ie.faceLabel()] != -1)
+ {
+ doneEFaces[findIndex(eFaces, ie.faceLabel())] = true;
+ verts.append(faceToDualPoint_[ie.faceLabel()]);
+ }
+ if (cellToDualPoint_[ie.cellLabel()] != -1)
+ {
+ verts.append(cellToDualPoint_[ie.cellLabel()]);
+ }
+
+ label currentDualCell0 = findDualCell(ie.cellLabel(), e[0]);
+ label currentDualCell1 = findDualCell(ie.cellLabel(), e[1]);
+
+ ++ie;
+
+ while (true)
+ {
+ label faceI = ie.faceLabel();
+
+ // Mark face as visited.
+ doneEFaces[findIndex(eFaces, faceI)] = true;
+
+ if (faceToDualPoint_[faceI] != -1)
+ {
+ verts.append(faceToDualPoint_[faceI]);
+ }
+
+ label cellI = ie.cellLabel();
+
+ if (cellI == -1)
+ {
+ // At ending boundary face. We've stored the face point above
+ // so this is the whole face.
+ break;
+ }
+
+
+ label dualCell0 = findDualCell(cellI, e[0]);
+ label dualCell1 = findDualCell(cellI, e[1]);
+
+ // Generate face. (always if splitFace=true; only if needed to
+ // separate cells otherwise)
+ if
+ (
+ splitFace
+ || (
+ dualCell0 != currentDualCell0
+ || dualCell1 != currentDualCell1
+ )
+ )
+ {
+ // Close current face.
+ addInternalFace
+ (
+ -1, // masterPointI
+ edgeI, // masterEdgeI
+ -1, // masterFaceI
+ edgeOrder,
+ currentDualCell0,
+ currentDualCell1,
+ verts.shrink(),
+ meshMod
+ );
+
+ // Restart
+ currentDualCell0 = dualCell0;
+ currentDualCell1 = dualCell1;
+
+ verts.clear();
+ if (edgeToDualPoint_[edgeI] != -1)
+ {
+ verts.append(edgeToDualPoint_[edgeI]);
+ }
+ if (faceToDualPoint_[faceI] != -1)
+ {
+ verts.append(faceToDualPoint_[faceI]);
+ }
+ }
+
+ if (cellToDualPoint_[cellI] != -1)
+ {
+ verts.append(cellToDualPoint_[cellI]);
+ }
+
+ ++ie;
+
+ if (ie == ie.end())
+ {
+ // Back at start face (for internal edge only). See if this needs
+ // adding.
+ if (isBoundaryEdge.get(edgeI) == 0)
+ {
+ label startDual = faceToDualPoint_[startFaceLabel];
+
+ if (startDual != -1 && findIndex(verts, startDual) == -1)
+ {
+ verts.append(startDual);
+ }
+ }
+ break;
+ }
+ }
+
+ verts.shrink();
+ addInternalFace
+ (
+ -1, // masterPointI
+ edgeI, // masterEdgeI
+ -1, // masterFaceI
+ edgeOrder,
+ currentDualCell0,
+ currentDualCell1,
+ verts,
+ meshMod
+ );
+}
+
+
+// Walks around circumference of faceI. Creates single face. Gets given
+// starting (feature) edge to start from. Returns ending edge. (all edges
+// in form of index in faceEdges)
+void Foam::meshDualiser::createFaceFromInternalFace
+(
+ const label faceI,
+ label& fp,
+ polyTopoChange& meshMod
+) const
+{
+ const face& f = mesh_.faces()[faceI];
+ const labelList& fEdges = mesh_.faceEdges()[faceI];
+ label own = mesh_.faceOwner()[faceI];
+ label nei = mesh_.faceNeighbour()[faceI];
+
+ //Pout<< "createFaceFromInternalFace : At face:" << faceI
+ // << " verts:" << f
+ // << " points:" << IndirectList(mesh_.points(), f)()
+ // << " started walking at edge:" << fEdges[fp]
+ // << " verts:" << mesh_.edges()[fEdges[fp]]
+ // << endl;
+
+
+ // Walk and collect face.
+ DynamicList verts(100);
+
+ verts.append(faceToDualPoint_[faceI]);
+ verts.append(edgeToDualPoint_[fEdges[fp]]);
+
+ // Step to vertex after edge mid
+ fp = f.fcIndex(fp);
+
+ // Get cells on either side of face at that point
+ label currentDualCell0 = findDualCell(own, f[fp]);
+ label currentDualCell1 = findDualCell(nei, f[fp]);
+
+ forAll(f, i)
+ {
+ // Check vertex
+ if (pointToDualPoint_[f[fp]] != -1)
+ {
+ verts.append(pointToDualPoint_[f[fp]]);
+ }
+
+ // Edge between fp and fp+1
+ label edgeI = fEdges[fp];
+
+ if (edgeToDualPoint_[edgeI] != -1)
+ {
+ verts.append(edgeToDualPoint_[edgeI]);
+ }
+
+ // Next vertex on edge
+ label nextFp = f.fcIndex(fp);
+
+ // Get dual cells on nextFp to check whether face needs closing.
+ label dualCell0 = findDualCell(own, f[nextFp]);
+ label dualCell1 = findDualCell(nei, f[nextFp]);
+
+ if (dualCell0 != currentDualCell0 || dualCell1 != currentDualCell1)
+ {
+ // Check: make sure that there is a midpoint on the edge.
+ if (edgeToDualPoint_[edgeI] == -1)
+ {
+ FatalErrorIn("createFacesFromInternalFace(..)")
+ << "face:" << faceI << " verts:" << f
+ << " points:" << IndirectList(mesh_.points(), f)()
+ << " no feature edge between " << f[fp]
+ << " and " << f[nextFp] << " although have different"
+ << " dual cells." << endl
+ << "point " << f[fp] << " has dual cells "
+ << currentDualCell0 << " and " << currentDualCell1
+ << " ; point "<< f[nextFp] << " has dual cells "
+ << dualCell0 << " and " << dualCell1
+ << abort(FatalError);
+ }
+
+
+ // Close current face.
+ verts.shrink();
+ addInternalFace
+ (
+ -1, // masterPointI
+ -1, // masterEdgeI
+ faceI, // masterFaceI
+ true, // edgeOrder,
+ currentDualCell0,
+ currentDualCell1,
+ verts,
+ meshMod
+ );
+ break;
+ }
+
+ fp = nextFp;
+ }
+}
+
+
+// Given a point on a face converts the faces around the point.
+// (pointFaces()). Gets starting face and marks off visited faces in donePFaces.
+void Foam::meshDualiser::createFacesAroundBoundaryPoint
+(
+ const label patchI,
+ const label patchPointI,
+ const label startFaceI,
+ polyTopoChange& meshMod,
+ boolList& donePFaces // pFaces visited
+) const
+{
+ const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+ const polyPatch& pp = patches[patchI];
+ const labelList& pFaces = pp.pointFaces()[patchPointI];
+ const labelList& own = mesh_.faceOwner();
+
+ label pointI = pp.meshPoints()[patchPointI];
+
+ if (pointToDualPoint_[pointI] == -1)
+ {
+ // Not a feature point. Loop over all connected
+ // pointFaces.
+
+ // Starting face
+ label faceI = startFaceI;
+
+ DynamicList verts(4);
+
+ while (true)
+ {
+ label index = findIndex(pFaces, faceI-pp.start());
+
+ // Has face been visited already?
+ if (donePFaces[index])
+ {
+ break;
+ }
+ donePFaces[index] = true;
+
+ // Insert face centre
+ verts.append(faceToDualPoint_[faceI]);
+
+ label dualCellI = findDualCell(own[faceI], pointI);
+
+ // Get the edge before the patchPointI
+ const face& f = mesh_.faces()[faceI];
+ label fp = findIndex(f, pointI);
+ label prevFp = f.rcIndex(fp);
+ label edgeI = mesh_.faceEdges()[faceI][prevFp];
+
+ if (edgeToDualPoint_[edgeI] != -1)
+ {
+ verts.append(edgeToDualPoint_[edgeI]);
+ }
+
+ // Get next boundary face (whilst staying on edge).
+ edgeFaceCirculator circ
+ (
+ mesh_,
+ faceI,
+ true, // ownerSide
+ prevFp, // index of edge in face
+ true // isBoundaryEdge
+ );
+
+ do
+ {
+ ++circ;
+ }
+ while (mesh_.isInternalFace(circ.faceLabel()));
+
+ // Step to next face
+ faceI = circ.faceLabel();
+
+ if (faceI < pp.start() || faceI >= pp.start()+pp.size())
+ {
+ FatalErrorIn
+ (
+ "meshDualiser::createFacesAroundBoundaryPoint(..)"
+ ) << "Walked from face on patch:" << patchI
+ << " to face:" << faceI
+ << " fc:" << mesh_.faceCentres()[faceI]
+ << " on patch:" << patches.whichPatch(faceI)
+ << abort(FatalError);
+ }
+
+ // Check if different cell.
+ if (dualCellI != findDualCell(own[faceI], pointI))
+ {
+ FatalErrorIn
+ (
+ "meshDualiser::createFacesAroundBoundaryPoint(..)"
+ ) << "Different dual cells but no feature edge"
+ << " inbetween point:" << pointI
+ << " coord:" << mesh_.points()[pointI]
+ << abort(FatalError);
+ }
+ }
+
+ verts.shrink();
+
+ label dualCellI = findDualCell(own[faceI], pointI);
+
+ //Bit dodgy: create dualface from the last face (instead of from
+ // the central point). This will also use the original faceZone to
+ // put the new face (which might span multiple original faces) in.
+
+ addBoundaryFace
+ (
+ //pointI, // masterPointI
+ -1, // masterPointI
+ -1, // masterEdgeI
+ faceI, // masterFaceI
+ dualCellI,
+ patchI,
+ verts,
+ meshMod
+ );
+ }
+ else
+ {
+ label faceI = startFaceI;
+
+ // Storage for face
+ DynamicList verts(mesh_.faces()[faceI].size());
+
+ // Starting point.
+ verts.append(pointToDualPoint_[pointI]);
+
+ // Find edge between pointI and next point on face.
+ const labelList& fEdges = mesh_.faceEdges()[faceI];
+ label nextEdgeI = fEdges[findIndex(mesh_.faces()[faceI], pointI)];
+ if (edgeToDualPoint_[nextEdgeI] != -1)
+ {
+ verts.append(edgeToDualPoint_[nextEdgeI]);
+ }
+
+ do
+ {
+ label index = findIndex(pFaces, faceI-pp.start());
+
+ // Has face been visited already?
+ if (donePFaces[index])
+ {
+ break;
+ }
+ donePFaces[index] = true;
+
+ // Face centre
+ verts.append(faceToDualPoint_[faceI]);
+
+ // Find edge before pointI on faceI
+ const labelList& fEdges = mesh_.faceEdges()[faceI];
+ const face& f = mesh_.faces()[faceI];
+ label prevFp = f.rcIndex(findIndex(f, pointI));
+ label edgeI = fEdges[prevFp];
+
+ if (edgeToDualPoint_[edgeI] != -1)
+ {
+ // Feature edge. Close any face so far. Note: uses face to
+ // create dualFace from. Could use pointI instead.
+ verts.append(edgeToDualPoint_[edgeI]);
+ addBoundaryFace
+ (
+ -1, // masterPointI
+ -1, // masterEdgeI
+ faceI, // masterFaceI
+ findDualCell(own[faceI], pointI),
+ patchI,
+ verts.shrink(),
+ meshMod
+ );
+ verts.clear();
+
+ verts.append(pointToDualPoint_[pointI]);
+ verts.append(edgeToDualPoint_[edgeI]);
+ }
+
+ // Cross edgeI to next boundary face
+ edgeFaceCirculator circ
+ (
+ mesh_,
+ faceI,
+ true, // ownerSide
+ prevFp, // index of edge in face
+ true // isBoundaryEdge
+ );
+
+ do
+ {
+ ++circ;
+ }
+ while (mesh_.isInternalFace(circ.faceLabel()));
+
+ // Step to next face. Quit if not on same patch.
+ faceI = circ.faceLabel();
+ }
+ while
+ (
+ faceI != startFaceI
+ && faceI >= pp.start()
+ && faceI < pp.start()+pp.size()
+ );
+
+ if (verts.size() > 2)
+ {
+ // Note: face created from face, not from pointI
+ addBoundaryFace
+ (
+ -1, // masterPointI
+ -1, // masterEdgeI
+ startFaceI, // masterFaceI
+ findDualCell(own[faceI], pointI),
+ patchI,
+ verts.shrink(),
+ meshMod
+ );
+ }
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::meshDualiser::meshDualiser(const polyMesh& mesh)
+:
+ mesh_(mesh),
+ pointToDualCells_(mesh_.nPoints()),
+ pointToDualPoint_(mesh_.nPoints(), -1),
+ cellToDualPoint_(mesh_.nCells()),
+ faceToDualPoint_(mesh_.nFaces(), -1),
+ edgeToDualPoint_(mesh_.nEdges(), -1)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::meshDualiser::setRefinement
+(
+ const bool splitFace,
+ const labelList& featureFaces,
+ const labelList& featureEdges,
+ const labelList& singleCellFeaturePoints,
+ const labelList& multiCellFeaturePoints,
+ polyTopoChange& meshMod
+)
+{
+ const labelList& own = mesh_.faceOwner();
+ const labelList& nei = mesh_.faceNeighbour();
+ const vectorField& cellCentres = mesh_.cellCentres();
+
+ // 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);
+ }
+ }
+
+
+ if (splitFace)
+ {
+ // This is a special mode where whenever we are walking around an edge
+ // every area through a cell becomes a separate dualface. So two
+ // dual cells will probably have more than one dualface between them!
+ // This mode implies that
+ // - all faces have to be feature faces since there has to be a
+ // dualpoint at the face centre.
+ // - all edges have to be feature edges ,,
+ boolList featureFaceSet(mesh_.nFaces(), false);
+ forAll(featureFaces, i)
+ {
+ featureFaceSet[featureFaces[i]] = true;
+ }
+ label faceI = findIndex(featureFaceSet, false);
+
+ if (faceI != -1)
+ {
+ FatalErrorIn
+ (
+ "meshDualiser::setRefinement"
+ "(const labelList&, const labelList&, const labelList&"
+ ", const labelList, polyTopoChange&)"
+ ) << "In split-face-mode (splitFace=true) but not all faces"
+ << " marked as feature faces." << endl
+ << "First conflicting face:" << faceI
+ << " centre:" << mesh_.faceCentres()[faceI]
+ << abort(FatalError);
+ }
+
+ boolList featureEdgeSet(mesh_.nEdges(), false);
+ forAll(featureEdges, i)
+ {
+ featureEdgeSet[featureEdges[i]] = true;
+ }
+ label edgeI = findIndex(featureEdgeSet, false);
+
+ if (edgeI != -1)
+ {
+ const edge& e = mesh_.edges()[edgeI];
+ FatalErrorIn
+ (
+ "meshDualiser::setRefinement"
+ "(const labelList&, const labelList&, const labelList&"
+ ", const labelList, polyTopoChange&)"
+ ) << "In split-face-mode (splitFace=true) but not all edges"
+ << " marked as feature edges." << endl
+ << "First conflicting edge:" << edgeI
+ << " verts:" << e
+ << " coords:" << mesh_.points()[e[0]] << mesh_.points()[e[1]]
+ << abort(FatalError);
+ }
+ }
+ else
+ {
+ // Check that all boundary faces are feature faces.
+
+ boolList featureFaceSet(mesh_.nFaces(), false);
+ forAll(featureFaces, i)
+ {
+ featureFaceSet[featureFaces[i]] = true;
+ }
+ for
+ (
+ label faceI = mesh_.nInternalFaces();
+ faceI < mesh_.nFaces();
+ faceI++
+ )
+ {
+ if (!featureFaceSet[faceI])
+ {
+ FatalErrorIn
+ (
+ "meshDualiser::setRefinement"
+ "(const labelList&, const labelList&, const labelList&"
+ ", const labelList, polyTopoChange&)"
+ ) << "Not all boundary faces marked as feature faces."
+ << endl
+ << "First conflicting face:" << faceI
+ << " centre:" << mesh_.faceCentres()[faceI]
+ << abort(FatalError);
+ }
+ }
+ }
+
+
+
+
+ // Start creating cells, points, and faces (in that order)
+
+
+ // 1. Mark which cells to create
+ // Mostly every point becomes one cell but sometimes (for feature points)
+ // all cells surrounding a feature point become cells. Also a non-manifold
+ // point can create two cells! So a dual cell is uniquely defined by a
+ // mesh point + cell (as in pointCells index)
+
+ // 2. Mark which face centres to create
+
+ // 3. Internal faces can now consist of
+ // - only cell centres of walk around edge
+ // - cell centres + face centres of walk around edge
+ // - same but now other side is not a single cell
+
+ // 4. Boundary faces (or internal faces between cell zones!) now consist of
+ // - walk around boundary point.
+
+
+
+ autoPtr dualCcStr;
+ if (debug)
+ {
+ dualCcStr.reset(new OFstream("dualCc.obj"));
+ Pout<< "Dumping centres of dual cells to " << dualCcStr().name()
+ << endl;
+ }
+
+
+ // Dual cells (from points)
+ // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+ // pointToDualCells_[pointI]
+ // - single entry : all cells surrounding point all become the same
+ // cell.
+ // - multiple entries: in order of pointCells.
+
+
+ // feature points that become single cell
+ forAll(singleCellFeaturePoints, i)
+ {
+ label pointI = singleCellFeaturePoints[i];
+
+ pointToDualPoint_[pointI] = meshMod.addPoint
+ (
+ mesh_.points()[pointI],
+ pointI, // masterPoint
+ mesh_.pointZones().whichZone(pointI), // zoneID
+ true // inCell
+ );
+
+ // Generate single cell
+ pointToDualCells_[pointI].setSize(1);
+ pointToDualCells_[pointI][0] = meshMod.addCell
+ (
+ pointI, //masterPointID,
+ -1, //masterEdgeID,
+ -1, //masterFaceID,
+ -1, //masterCellID,
+ -1 //zoneID
+ );
+ if (dualCcStr.valid())
+ {
+ meshTools::writeOBJ(dualCcStr(), mesh_.points()[pointI]);
+ }
+ }
+
+ // feature points that become multiple cells
+ forAll(multiCellFeaturePoints, i)
+ {
+ label pointI = multiCellFeaturePoints[i];
+
+ if (pointToDualCells_[pointI].size() > 0)
+ {
+ FatalErrorIn
+ (
+ "meshDualiser::setRefinement"
+ "(const labelList&, const labelList&, const labelList&"
+ ", const labelList, polyTopoChange&)"
+ ) << "Point " << pointI << " at:" << mesh_.points()[pointI]
+ << " is both in singleCellFeaturePoints"
+ << " and multiCellFeaturePoints."
+ << abort(FatalError);
+ }
+
+ pointToDualPoint_[pointI] = meshMod.addPoint
+ (
+ mesh_.points()[pointI],
+ pointI, // masterPoint
+ mesh_.pointZones().whichZone(pointI), // zoneID
+ true // inCell
+ );
+
+ // Create dualcell for every cell connected to dual point
+
+ const labelList& pCells = mesh_.pointCells()[pointI];
+
+ pointToDualCells_[pointI].setSize(pCells.size());
+
+ forAll(pCells, pCellI)
+ {
+ pointToDualCells_[pointI][pCellI] = meshMod.addCell
+ (
+ pointI, //masterPointID
+ -1, //masterEdgeID
+ -1, //masterFaceID
+ -1, //masterCellID
+ mesh_.cellZones().whichZone(pCells[pCellI]) //zoneID
+ );
+ if (dualCcStr.valid())
+ {
+ meshTools::writeOBJ
+ (
+ dualCcStr(),
+ 0.5*(mesh_.points()[pointI]+cellCentres[pCells[pCellI]])
+ );
+ }
+ }
+ }
+ // Normal points
+ forAll(mesh_.points(), pointI)
+ {
+ if (pointToDualCells_[pointI].size() == 0)
+ {
+ pointToDualCells_[pointI].setSize(1);
+ pointToDualCells_[pointI][0] = meshMod.addCell
+ (
+ pointI, //masterPointID,
+ -1, //masterEdgeID,
+ -1, //masterFaceID,
+ -1, //masterCellID,
+ -1 //zoneID
+ );
+
+ if (dualCcStr.valid())
+ {
+ meshTools::writeOBJ(dualCcStr(), mesh_.points()[pointI]);
+ }
+ }
+ }
+
+
+ // Dual points (from cell centres, feature faces, feature edges)
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ forAll(cellToDualPoint_, cellI)
+ {
+ cellToDualPoint_[cellI] = meshMod.addPoint
+ (
+ cellCentres[cellI],
+ mesh_.faces()[mesh_.cells()[cellI][0]][0], // masterPoint
+ -1, // zoneID
+ true // inCell
+ );
+ }
+
+ // From face to dual point
+
+ forAll(featureFaces, i)
+ {
+ label faceI = featureFaces[i];
+
+ faceToDualPoint_[faceI] = meshMod.addPoint
+ (
+ mesh_.faceCentres()[faceI],
+ mesh_.faces()[faceI][0], // masterPoint
+ -1, // zoneID
+ true // inCell
+ );
+ }
+ // Detect whether different dual cells on either side of a face. This
+ // would neccesitate having a dual face built from the face and thus a
+ // dual point at the face centre.
+ for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+ {
+ if (faceToDualPoint_[faceI] == -1)
+ {
+ const face& f = mesh_.faces()[faceI];
+
+ forAll(f, fp)
+ {
+ label ownDualCell = findDualCell(own[faceI], f[fp]);
+ label neiDualCell = findDualCell(nei[faceI], f[fp]);
+
+ if (ownDualCell != neiDualCell)
+ {
+ faceToDualPoint_[faceI] = meshMod.addPoint
+ (
+ mesh_.faceCentres()[faceI],
+ f[fp], // masterPoint
+ -1, // zoneID
+ true // inCell
+ );
+
+ break;
+ }
+ }
+ }
+ }
+
+ // From edge to dual point
+
+ forAll(featureEdges, i)
+ {
+ label edgeI = featureEdges[i];
+
+ const edge& e = mesh_.edges()[edgeI];
+
+ edgeToDualPoint_[edgeI] = meshMod.addPoint
+ (
+ e.centre(mesh_.points()),
+ e[0], // masterPoint
+ -1, // zoneID
+ true // inCell
+ );
+ }
+
+ // Detect whether different dual cells on either side of an edge. This
+ // would neccesitate having a dual face built perpendicular to the edge
+ // and thus a dual point at the mid of the edge.
+ // Note: not really true - the face can be built without the edge centre!
+ const labelListList& edgeCells = mesh_.edgeCells();
+
+ forAll(edgeCells, edgeI)
+ {
+ if (edgeToDualPoint_[edgeI] == -1)
+ {
+ const edge& e = mesh_.edges()[edgeI];
+
+ // We need a point on the edge if not all cells on both sides
+ // are the same.
+
+ const labelList& eCells = mesh_.edgeCells()[edgeI];
+
+ label dualE0 = findDualCell(eCells[0], e[0]);
+ label dualE1 = findDualCell(eCells[0], e[1]);
+
+ for (label i = 1; i < eCells.size(); i++)
+ {
+ label newDualE0 = findDualCell(eCells[i], e[0]);
+
+ if (dualE0 != newDualE0)
+ {
+ edgeToDualPoint_[edgeI] = meshMod.addPoint
+ (
+ e.centre(mesh_.points()),
+ e[0], // masterPoint
+ mesh_.pointZones().whichZone(e[0]), // zoneID
+ true // inCell
+ );
+
+ break;
+ }
+
+ label newDualE1 = findDualCell(eCells[i], e[1]);
+
+ if (dualE1 != newDualE1)
+ {
+ edgeToDualPoint_[edgeI] = meshMod.addPoint
+ (
+ e.centre(mesh_.points()),
+ e[1], // masterPoint
+ mesh_.pointZones().whichZone(e[1]), // zoneID
+ true // inCell
+ );
+
+ break;
+ }
+ }
+ }
+ }
+
+ // Make sure all boundary edges emanating from feature points are
+ // feature edges as well.
+ forAll(singleCellFeaturePoints, i)
+ {
+ generateDualBoundaryEdges
+ (
+ isBoundaryEdge,
+ singleCellFeaturePoints[i],
+ meshMod
+ );
+ }
+ forAll(multiCellFeaturePoints, i)
+ {
+ generateDualBoundaryEdges
+ (
+ isBoundaryEdge,
+ multiCellFeaturePoints[i],
+ meshMod
+ );
+ }
+
+
+ // Check for duplicate points
+ if (debug)
+ {
+ dumpPolyTopoChange(meshMod, "generatedPoints_");
+ checkPolyTopoChange(meshMod);
+ }
+
+
+ // Now we have all points and cells
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // - pointToDualCells_ : per point a single dualCell or multiple dualCells
+ // - pointToDualPoint_ : per point -1 or the dual point at the coordinate
+ // - edgeToDualPoint_ : per edge -1 or the edge centre
+ // - faceToDualPoint_ : per face -1 or the face centre
+ // - cellToDualPoint_ : per cell the cell centre
+ // Now we have to walk all edges and construct faces. Either single face
+ // per edge or multiple (-if nonmanifold edge -if different dualcells)
+
+ const edgeList& edges = mesh_.edges();
+
+ forAll(edges, edgeI)
+ {
+ const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+
+ boolList doneEFaces(eFaces.size(), false);
+
+ forAll(eFaces, i)
+ {
+ if (!doneEFaces[i])
+ {
+ // We found a face that hasn't yet been visited. This might
+ // happen for non-manifold edges where a single edge can
+ // become multiple faces.
+
+ label startFaceI = eFaces[i];
+
+ //Pout<< "Walking edge:" << edgeI
+ // << " points:" << mesh_.points()[e[0]]
+ // << mesh_.points()[e[1]]
+ // << " startFace:" << startFaceI
+ // << " at:" << mesh_.faceCentres()[startFaceI]
+ // << endl;
+
+ createFacesAroundEdge
+ (
+ splitFace,
+ isBoundaryEdge,
+ edgeI,
+ startFaceI,
+ meshMod,
+ doneEFaces
+ );
+ }
+ }
+ }
+
+ if (debug)
+ {
+ dumpPolyTopoChange(meshMod, "generatedFacesFromEdges_");
+ }
+
+ // Create faces from feature faces. These can be internal or external faces.
+ // - feature face : centre needs to be included.
+ // - single cells on either side: triangulate
+ // - multiple cells: create single face between unique cell pair. Only
+ // create face where cells differ on either side.
+ // - non-feature face : inbetween cell zones.
+ forAll(faceToDualPoint_, faceI)
+ {
+ if (faceToDualPoint_[faceI] != -1 && mesh_.isInternalFace(faceI))
+ {
+ const face& f = mesh_.faces()[faceI];
+ const labelList& fEdges = mesh_.faceEdges()[faceI];
+
+ // Starting edge
+ label fp = 0;
+
+ do
+ {
+ // Find edge that is in dual mesh and where the
+ // next point (fp+1) has different dual cells on either side.
+ bool foundStart = false;
+
+ do
+ {
+ if
+ (
+ edgeToDualPoint_[fEdges[fp]] != -1
+ && !sameDualCell(faceI, f.nextLabel(fp))
+ )
+ {
+ foundStart = true;
+ break;
+ }
+ fp = f.fcIndex(fp);
+ }
+ while (fp != 0);
+
+ if (!foundStart)
+ {
+ break;
+ }
+
+ // Walk from edge fp and generate a face.
+ createFaceFromInternalFace
+ (
+ faceI,
+ fp,
+ meshMod
+ );
+ }
+ while(fp != 0);
+ }
+ }
+
+ if (debug)
+ {
+ dumpPolyTopoChange(meshMod, "generatedFacesFromFeatFaces_");
+ }
+
+
+ // Create boundary faces. Every boundary point has one or more dualcells.
+ // These need to be closed.
+ const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+ forAll(patches, patchI)
+ {
+ const polyPatch& pp = patches[patchI];
+
+ const labelListList& pointFaces = pp.pointFaces();
+
+ forAll(pointFaces, patchPointI)
+ {
+ const labelList& pFaces = pointFaces[patchPointI];
+
+ boolList donePFaces(pFaces.size(), false);
+
+ forAll(pFaces, i)
+ {
+ if (!donePFaces[i])
+ {
+ // Starting face
+ label startFaceI = pp.start()+pFaces[i];
+
+ //Pout<< "Walking around point:" << pointI
+ // << " coord:" << mesh_.points()[pointI]
+ // << " on patch:" << patchI
+ // << " startFace:" << startFaceI
+ // << " at:" << mesh_.faceCentres()[startFaceI]
+ // << endl;
+
+ createFacesAroundBoundaryPoint
+ (
+ patchI,
+ patchPointI,
+ startFaceI,
+ meshMod,
+ donePFaces // pFaces visited
+ );
+ }
+ }
+ }
+ }
+
+ if (debug)
+ {
+ dumpPolyTopoChange(meshMod, "generatedFacesFromBndFaces_");
+ }
+}
+
+
+
+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/conversion/polyDualMesh/meshDualiser.H b/applications/utilities/mesh/conversion/polyDualMesh/meshDualiser.H
new file mode 100644
index 0000000000..ee5599816c
--- /dev/null
+++ b/applications/utilities/mesh/conversion/polyDualMesh/meshDualiser.H
@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+
+Description
+ Creates dual of polyMesh. Every point becomes a cell (or multiple cells
+ for feature points), a walk around every edge creates faces between them.
+
+ Put all points you want in the final mesh into featurePoints; all edge(mid)s
+ you want in the final mesh into featureEdges; all face(centre)s in
+ faceFaces.
+
+ Usually to preserve boundaries:
+ - all boundary faces are featureFaces
+ - all edges and points inbetween different patches are
+ featureEdges/points.
+
+ In same way you can also preserve internal faces (e.g. faceZones)
+
+SourceFiles
+ meshDualiser.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef meshDualiser_H
+#define meshDualiser_H
+
+#include "DynamicList.H"
+#include "PackedList.H"
+#include "boolList.H"
+#include "typeInfo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+class polyTopoChange;
+
+/*---------------------------------------------------------------------------*\
+ Class meshDualiser Declaration
+\*---------------------------------------------------------------------------*/
+
+class meshDualiser
+{
+ // Private data
+
+ const polyMesh& mesh_;
+
+ //- From point on cell to dual cell. Either single entry or
+ // one entry per pointCells
+ labelListList pointToDualCells_;
+
+ //- From point to dual point (or -1 if not feature point).
+ labelList pointToDualPoint_;
+
+ //- From cell to dual point. All cells become point
+ labelList cellToDualPoint_;
+
+ //- From face to dual point (or -1 if not feature face)
+ labelList faceToDualPoint_;
+
+ //- From edge to dual point (or -1 if not feature edge)
+ labelList edgeToDualPoint_;
+
+
+ // Private Member Functions
+
+ static void checkPolyTopoChange(const polyTopoChange&);
+
+ static void dumpPolyTopoChange(const polyTopoChange&, const fileName&);
+
+ //- Find dual cell given point and cell
+ label findDualCell(const label cellI, const label pointI) const;
+
+ //- Helper function to generate dualpoints on all boundary edges
+ // emanating from (boundary & feature) point
+ void generateDualBoundaryEdges
+ (
+ const PackedList<1>&,
+ const label pointI,
+ polyTopoChange&
+ );
+
+ //- Check that owner and neighbour of face have same dual cell
+ bool sameDualCell
+ (
+ const label faceI,
+ const label pointI
+ ) const;
+
+ //- Add internal face
+ label addInternalFace
+ (
+ const label masterPointI,
+ const label masterEdgeI,
+ const label masterFaceI,
+
+ const bool edgeOrder,
+ const label dualCell0,
+ const label dualCell1,
+ const DynamicList& verts,
+ polyTopoChange& meshMod
+ ) const;
+
+ //- Add boundary face
+ label addBoundaryFace
+ (
+ const label masterPointI,
+ const label masterEdgeI,
+ const label masterFaceI,
+
+ const label dualCellI,
+ const label patchI,
+ const DynamicList& verts,
+ polyTopoChange& meshMod
+ ) const;
+
+ //- Create internal faces walking around edge
+ void createFacesAroundEdge
+ (
+ const bool splitFace,
+ const PackedList<1>&,
+ const label edgeI,
+ const label startFaceI,
+ polyTopoChange&,
+ boolList& doneEFaces
+ ) const;
+
+ //- Create single internal face from internal face
+ void createFaceFromInternalFace
+ (
+ const label faceI,
+ label& fp,
+ polyTopoChange&
+ ) const;
+
+ //- Creates boundary faces walking around point on patchI.
+ void createFacesAroundBoundaryPoint
+ (
+ const label patchI,
+ const label patchPointI,
+ const label startFaceI,
+ polyTopoChange&,
+ boolList& donePFaces // pFaces visited
+ ) const;
+
+ //- Disallow default bitwise copy construct
+ meshDualiser(const meshDualiser&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const meshDualiser&);
+
+
+public:
+
+ //- Runtime type information
+ ClassName("meshDualiser");
+
+
+ // Constructors
+
+ //- Construct from mesh
+ meshDualiser(const polyMesh&);
+
+
+ // Member Functions
+
+ // Access
+
+ //- From point on cell to dual cell. Either single entry or
+ // one entry per pointCells.
+ const labelListList& pointToDualCells() const
+ {
+ return pointToDualCells_;
+ }
+
+ //- From point to dual point (or -1 if not feature point).
+ const labelList& pointToDualPoint() const
+ {
+ return pointToDualPoint_;
+ }
+
+ //- From cell to dual point (at cell centre). All cells become
+ // points.
+ const labelList& cellToDualPoint() const
+ {
+ return cellToDualPoint_;
+ }
+
+ //- From face to dual point (at face centre; or -1 if not
+ // feature face).
+ const labelList& faceToDualPoint() const
+ {
+ return faceToDualPoint_;
+ }
+
+ //- From edge to dual point (at edge mid; or -1 if not feature
+ // edge).
+ const labelList& edgeToDualPoint() const
+ {
+ return edgeToDualPoint_;
+ }
+
+
+ // Edit
+
+
+ //- Insert all changes into meshMod to convert the polyMesh into
+ // its dual.
+ // featureFaces : faces where we want a point at the face centre
+ // featureEdges : edges ,, edge mid
+ // featurePoints : points ,, point. Two variants:
+ // singleCellFeaturePoints : point becomes one dualcell.
+ // Use this for e.g. convex boundary points.
+ // multiCellFeaturePoints : one dualcell per original cell
+ // around point. Use this for e.g. concave boundary points
+ // since it prevents big concave boundary cells.
+ void setRefinement
+ (
+ const bool splitFace,
+ const labelList& featureFaces,
+ const labelList& featureEdges,
+ const labelList& singleCellFeaturePoints,
+ const labelList& multiCellFeaturePoints,
+ polyTopoChange& meshMod
+ );
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C b/applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C
deleted file mode 100644
index ffa8314952..0000000000
--- a/applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C
+++ /dev/null
@@ -1,79 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 1991-2008 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
-
-Application
- polyDualMesh
-
-Description
- Calculate the dual of a polyMesh. Adheres to all the feature&patch edges
-
-\*---------------------------------------------------------------------------*/
-
-#include "argList.H"
-#include "Time.H"
-#include "polyDualMesh.H"
-#include "mathematicalConstants.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
- argList::noParallel();
- argList::validArgs.append("feature angle[0-180]");
- argList::validOptions.insert("overwrite", "");
-
-# include "setRootCase.H"
-# include "createTime.H"
- runTime.functionObjects().off();
-# include "createPolyMesh.H"
-
- scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
- bool overwrite = args.options().found("overwrite");
-
- scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0);
-
- Info<< "Feature:" << featureAngle << endl
- << "minCos :" << minCos << endl
- << endl;
-
- polyDualMesh dualMesh(mesh, minCos);
-
- if (!overwrite)
- {
- runTime++;
- }
-
- Pout<< "Writing dualMesh to " << runTime.timeName() << endl;
-
- dualMesh.write();
-
- Info << "End\n" << endl;
-
- return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index ac924136e1..2596ae6150 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -109,7 +109,7 @@ Foam::label Foam::checkTopology
{
Info<< " Number of regions: " << rs.nRegions() << " (OK)."
<< endl;
-
+
}
else
{
@@ -214,7 +214,7 @@ Foam::label Foam::checkTopology
const pointField& pts = pp.points();
const labelList& mp = pp.meshPoints();
- boundBox bb(vector::zero, vector::zero);
+ boundBox bb; // zero-sized
if (returnReduce(mp.size(), sumOp()) > 0)
{
bb.min() = pts[mp[0]];
diff --git a/applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C b/applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C
index bcd1c5c9ce..99eea4a2ca 100644
--- a/applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C
+++ b/applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C
@@ -144,7 +144,7 @@ Foam::mergePolyMesh::mergePolyMesh(const IOobject& io)
wordList curPointZoneNames = pointZones().names();
if (curPointZoneNames.size() > 0)
{
- pointZoneNames_.setSize(2*curPointZoneNames.size());
+ pointZoneNames_.setCapacity(2*curPointZoneNames.size());
}
forAll (curPointZoneNames, zoneI)
@@ -157,7 +157,7 @@ Foam::mergePolyMesh::mergePolyMesh(const IOobject& io)
if (curFaceZoneNames.size() > 0)
{
- faceZoneNames_.setSize(2*curFaceZoneNames.size());
+ faceZoneNames_.setCapacity(2*curFaceZoneNames.size());
}
forAll (curFaceZoneNames, zoneI)
{
@@ -169,7 +169,7 @@ Foam::mergePolyMesh::mergePolyMesh(const IOobject& io)
if (curCellZoneNames.size() > 0)
{
- cellZoneNames_.setSize(2*curCellZoneNames.size());
+ cellZoneNames_.setCapacity(2*curCellZoneNames.size());
}
forAll (curCellZoneNames, zoneI)
{
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
index 548f8247dc..c24e779a83 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
@@ -354,7 +354,7 @@ bool domainDecomposition::writeDecomposition()
// Estimate size
forAll(zonePoints, zoneI)
{
- zonePoints[zoneI].setSize(pz[zoneI].size() / nProcs_);
+ zonePoints[zoneI].setCapacity(pz[zoneI].size() / nProcs_);
}
// Use the pointToZone map to find out the single zone (if any),
@@ -423,8 +423,8 @@ bool domainDecomposition::writeDecomposition()
{
label procSize = fz[zoneI].size() / nProcs_;
- zoneFaces[zoneI].setSize(procSize);
- zoneFaceFlips[zoneI].setSize(procSize);
+ zoneFaces[zoneI].setCapacity(procSize);
+ zoneFaceFlips[zoneI].setCapacity(procSize);
}
// Go through all the zoned faces and find out if they
@@ -514,7 +514,7 @@ bool domainDecomposition::writeDecomposition()
// Estimate size
forAll(zoneCells, zoneI)
{
- zoneCells[zoneI].setSize(cz[zoneI].size() / nProcs_);
+ zoneCells[zoneI].setCapacity(cz[zoneI].size() / nProcs_);
}
forAll (curCellLabels, celli)
diff --git a/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C b/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
index 834b75dfe6..87b134f6e9 100644
--- a/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
+++ b/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
@@ -273,7 +273,7 @@ autoPtr mergeSharedPoints
}
}
- return map;
+ return map;
}
@@ -418,11 +418,7 @@ int main(int argc, char *argv[])
// Read point on individual processors to determine merge tolerance
// (otherwise single cell domains might give problems)
- boundBox bb
- (
- point(GREAT, GREAT, GREAT),
- point(-GREAT, -GREAT, -GREAT)
- );
+ boundBox bb = boundBox::invertedBox;
for (label procI = 0; procI < nProcs; procI++)
{
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
index 7780ad3cfd..7ae5853896 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
@@ -617,8 +617,10 @@ void Foam::vtkPV3Foam::addPatchNames(vtkRenderer* renderer)
{
const labelList& eFaces = edgeFaces[edgeI];
- if (eFaces.size() != 2)
+ if (eFaces.size() == 1)
{
+ // Note: could also do ones with > 2 faces but this gives
+ // too many zones for baffles
featEdge[edgeI] = true;
}
else if (mag(n[eFaces[0]] & n[eFaces[1]]) < 0.5)
diff --git a/applications/utilities/postProcessing/wall/wallHeatFlux/wallHeatFlux.C b/applications/utilities/postProcessing/wall/wallHeatFlux/wallHeatFlux.C
index 1eb405cf29..49ac7d7431 100644
--- a/applications/utilities/postProcessing/wall/wallHeatFlux/wallHeatFlux.C
+++ b/applications/utilities/postProcessing/wall/wallHeatFlux/wallHeatFlux.C
@@ -64,7 +64,7 @@ int main(int argc, char *argv[])
Info<< "\nWall heat fluxes [W]" << endl;
forAll(patchHeatFlux, patchi)
{
- if (typeid(mesh.boundary()) == typeid(wallFvPatch))
+ if (typeid(mesh.boundary()[patchi]) == typeid(wallFvPatch))
{
Info<< mesh.boundary()[patchi].name()
<< " "
diff --git a/applications/utilities/preProcessing/mapFields/MapVolFields.H b/applications/utilities/preProcessing/mapFields/MapVolFields.H
index d72ae2948d..affbbb68f0 100644
--- a/applications/utilities/preProcessing/mapFields/MapVolFields.H
+++ b/applications/utilities/preProcessing/mapFields/MapVolFields.H
@@ -29,7 +29,6 @@ License
#include "GeometricField.H"
#include "meshToMesh.H"
-#include "cuttingPlane.H"
#include "IOobjectList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/surface/surfaceCoordinateSystemTransform/surfaceCoordinateSystemTransform.C b/applications/utilities/surface/surfaceCoordinateSystemTransform/surfaceCoordinateSystemTransform.C
index 8453d6f94a..a8b2f980ca 100644
--- a/applications/utilities/surface/surfaceCoordinateSystemTransform/surfaceCoordinateSystemTransform.C
+++ b/applications/utilities/surface/surfaceCoordinateSystemTransform/surfaceCoordinateSystemTransform.C
@@ -199,7 +199,6 @@ int main(int argc, char *argv[])
if (args.options().found("clean"))
{
surf.cleanup(true);
- surf.checkOrientation(true);
}
if (fromCsys.valid())
diff --git a/applications/utilities/surface/surfaceMeshConvert/surfaceMeshConvert.C b/applications/utilities/surface/surfaceMeshConvert/surfaceMeshConvert.C
index a907c7208c..0f18b5ac86 100644
--- a/applications/utilities/surface/surfaceMeshConvert/surfaceMeshConvert.C
+++ b/applications/utilities/surface/surfaceMeshConvert/surfaceMeshConvert.C
@@ -104,13 +104,19 @@ int main(int argc, char *argv[])
{
triSurface surf(importName);
+ Info<< "Read surface:" << endl;
+ surf.writeStats(Info);
+ Info<< endl;
+
if (args.options().found("clean"))
{
+ Info<< "Cleaning up surface" << endl;
surf.cleanup(true);
- surf.checkOrientation(true);
+ surf.writeStats(Info);
+ Info<< endl;
}
- Info << "writing " << exportName;
+ Info<< "writing " << exportName;
if (scaleFactor <= 0)
{
Info<< " without scaling" << endl;
@@ -119,6 +125,8 @@ int main(int argc, char *argv[])
{
Info<< " with scaling " << scaleFactor << endl;
surf.scalePoints(scaleFactor);
+ surf.writeStats(Info);
+ Info<< endl;
}
// write sorted by region
@@ -128,34 +136,16 @@ int main(int argc, char *argv[])
{
UnsortedMeshedSurface surf(importName);
- if (args.options().found("clean"))
- {
- surf.cleanup(true);
- surf.checkOrientation(true);
- }
-
- Info << "writing " << exportName;
- if (scaleFactor <= 0)
- {
- Info<< " without scaling" << endl;
- }
- else
- {
- Info<< " with scaling " << scaleFactor << endl;
- surf.scalePoints(scaleFactor);
- }
-
- surf.write(exportName);
- }
-#if 1
- else if (args.options().found("triFace"))
- {
- MeshedSurface surf(importName);
+ Info<< "Read surface:" << endl;
+ surf.writeStats(Info);
+ Info<< endl;
if (args.options().found("clean"))
{
+ Info<< "Cleaning up surface" << endl;
surf.cleanup(true);
- surf.checkOrientation(true);
+ surf.writeStats(Info);
+ Info<< endl;
}
Info<< "writing " << exportName;
@@ -167,6 +157,39 @@ int main(int argc, char *argv[])
{
Info<< " with scaling " << scaleFactor << endl;
surf.scalePoints(scaleFactor);
+ surf.writeStats(Info);
+ Info<< endl;
+ }
+ surf.write(exportName);
+ }
+#if 1
+ else if (args.options().found("triFace"))
+ {
+ MeshedSurface surf(importName);
+
+ Info<< "Read surface:" << endl;
+ surf.writeStats(Info);
+ Info<< endl;
+
+ if (args.options().found("clean"))
+ {
+ Info<< "Cleaning up surface" << endl;
+ surf.cleanup(true);
+ surf.writeStats(Info);
+ Info<< endl;
+ }
+
+ Info<< "writing " << exportName;
+ if (scaleFactor <= 0)
+ {
+ Info<< " without scaling" << endl;
+ }
+ else
+ {
+ Info<< " with scaling " << scaleFactor << endl;
+ surf.scalePoints(scaleFactor);
+ surf.writeStats(Info);
+ Info<< endl;
}
surf.write(exportName);
}
@@ -175,10 +198,16 @@ int main(int argc, char *argv[])
{
MeshedSurface surf(importName);
+ Info<< "Read surface:" << endl;
+ surf.writeStats(Info);
+ Info<< endl;
+
if (args.options().found("clean"))
{
+ Info<< "Cleaning up surface" << endl;
surf.cleanup(true);
- surf.checkOrientation(true);
+ surf.writeStats(Info);
+ Info<< endl;
}
Info<< "writing " << exportName;
@@ -190,6 +219,8 @@ int main(int argc, char *argv[])
{
Info<< " with scaling " << scaleFactor << endl;
surf.scalePoints(scaleFactor);
+ surf.writeStats(Info);
+ Info<< endl;
}
surf.write(exportName);
}
diff --git a/bin/doxyFilt b/bin/doxyFilt
index e2f8ee0a1e..c08851c53e 100755
--- a/bin/doxyFilt
+++ b/bin/doxyFilt
@@ -50,9 +50,9 @@ then
*/applications/solvers/*.C | */applications/utilities/*.C )
awkScript=$WM_PROJECT_DIR/bin/tools/doxyFilt-top.awk
;;
- */applications/solvers/*.H | */applications/utilities/*.H )
- awkScript=$WM_PROJECT_DIR/bin/tools/doxyFilt-ignore.awk
- ;;
+# */applications/solvers/*.H | */applications/utilities/*.H )
+# awkScript=$WM_PROJECT_DIR/bin/tools/doxyFilt-ignore.awk
+# ;;
esac
awk -f $awkScript $1 | \
diff --git a/bin/foamCheckSourceDeps b/bin/foamCheckSourceDeps
deleted file mode 100755
index 5653f050d1..0000000000
--- a/bin/foamCheckSourceDeps
+++ /dev/null
@@ -1,70 +0,0 @@
-#!/bin/sh
-#------------------------------------------------------------------------------
-# ========= |
-# \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
-# \\ / O peration |
-# \\ / A nd | Copyright (C) 1991-2008 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
-#
-# Script
-# foamCheckSourceDeps
-#
-# Description
-# Usage: foamCheckSourceDeps [dir1 .. dirN]
-#
-# Search for *.dep files that are without a corresponding .C or .L file.
-# These could indicate a directory that has been moved.
-# - print questionable directory and dep file
-#------------------------------------------------------------------------------
-if [ "$1" = "-h" -o "$1" = "-help" ]
-then
- cat <&2
-Usage: ${0##*/} [dir1 .. dirN]
-
- Search for .dep files that are without a corresponding .C or .L file.
- This could indicate a directory that has been moved.
- - print questionable directory and file
-USAGE
- exit 1
-fi
-
-if [ "$#" -eq 0 ]
-then
- set -- .
-fi
-
-for checkDir
-do
- if [ -d $checkDir ]
- then
- find $checkDir -name '*.dep' -print | while read depFile;
- do
- Cfile=$(echo $depFile | sed -e 's/\.dep$/.C/')
- # also check flex files
- Lfile=$(echo $depFile | sed -e 's/\.C$/.L/')
- if [ ! -f $Cfile -a ! -f $Lfile ]
- then
- echo "$(dirname $Cfile) ${depFile##*/}"
- fi
- done
- fi
-done
-
-# -----------------------------------------------------------------------------
diff --git a/bin/rmdepold b/bin/rmdepold
new file mode 100755
index 0000000000..72123a3cce
--- /dev/null
+++ b/bin/rmdepold
@@ -0,0 +1,109 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# ========= |
+# \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+# \\ / O peration |
+# \\ / A nd | Copyright (C) 1991-2008 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
+#
+# Script
+# rmdepold
+#
+# Description
+# Usage: rmdepold [dir1 .. dirN]
+#
+# Remove *.dep files that are without a corresponding .C or .L file.
+# This often occurs when a directory has been moved.
+# - print questionable directory and the *.dep file
+# - optionally remove empty directories
+#------------------------------------------------------------------------------
+usage() {
+ while [ "$#" -ge 1 ]; do echo "$1"; shift; done
+ cat<&2
+Usage: ${0##*/} [OPTION] [dir1 .. dirN]
+options:
+ -rmdir find and remove empty directories (recursively)
+
+Remove *.dep files that are without a corresponding .C or .L file.
+This often occurs when a directory has been moved.
+ - print questionable directory and file
+ - optionally remove empty directories
+
+USAGE
+ exit 1
+}
+
+unset optRmdir
+
+# parse options
+while [ "$#" -gt 0 ]
+do
+ case "$1" in
+ -h | -help)
+ usage
+ ;;
+ -rmdir)
+ optRmdir=true
+ shift
+ ;;
+ -*)
+ usage "unknown option: '$*'"
+ ;;
+ *)
+ break
+ ;;
+ esac
+done
+
+# default is the current directory
+[ "$#" -gt 0 ] || set -- .
+
+for checkDir
+do
+ if [ -d $checkDir ]
+ then
+ echo "searching: $checkDir"
+ else
+ echo "skipping non-dir: $checkDir"
+ continue
+ fi
+
+ find $checkDir -name '*.dep' -print | while read depFile
+ do
+ # check C++ and Flex files
+ if [ ! -r "${depFile%dep}C" -a ! -r "${depFile%dep}L" ];
+ then
+ echo "rm $depFile"
+ rm -f $depFile 2>/dev/null
+ fi
+ done
+
+ # remove empty dirs
+ if [ "$optRmdir" ]
+ then
+ # get subdirs ourselves so we can avoid particular directories
+ for dir in $(find $checkDir -mindepth 1 -maxdepth 1 -type d \( -name .git -prune -o -print \) )
+ do
+ echo "check dir: $dir"
+ find $dir -depth -type d -empty -exec rmdir {} \; -print
+ done
+ fi
+done
+# -----------------------------------------------------------------------------
diff --git a/doc/Doxygen/Doxyfile b/doc/Doxygen/Doxyfile
index f5ab120424..908fa543b9 100644
--- a/doc/Doxygen/Doxyfile
+++ b/doc/Doxygen/Doxyfile
@@ -538,6 +538,7 @@ EXCLUDE_SYMBOLS =
EXAMPLE_PATH =
+
# If the value of the EXAMPLE_PATH tag contains directories, you can use the
# EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
# and *.h) to filter out the source-files in the directories. If left
@@ -787,6 +788,10 @@ TREEVIEW_WIDTH = 250
# configuration options related to the LaTeX output
#---------------------------------------------------------------------------
+# Path for OpenCFD LaTeX macros
+
+@INCLUDE_PATH = $(WM_PROJECT_DIR)/doc/Doxygen/Macros/
+
# If the GENERATE_LATEX tag is set to YES (the default) Doxygen will
# generate Latex output.
@@ -824,7 +829,7 @@ PAPER_TYPE = a4wide
# The EXTRA_PACKAGES tag can be to specify one or more names of LaTeX
# packages that should be included in the LaTeX output.
-EXTRA_PACKAGES =
+EXTRA_PACKAGES = $(WM_PROJECT_DIR)/doc/Doxygen/Macros/tensorOperator
# The LATEX_HEADER tag can be used to specify a personal LaTeX header for
# the generated latex document. The header should contain everything until
diff --git a/doc/Doxygen/Macros/tensorOperator.sty b/doc/Doxygen/Macros/tensorOperator.sty
new file mode 100644
index 0000000000..2b027493a8
--- /dev/null
+++ b/doc/Doxygen/Macros/tensorOperator.sty
@@ -0,0 +1,129 @@
+%-----------------------------------------------------------------------------
+% ========= |
+% \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+% \\ / O peration |
+% \\ / A nd | Copyright (C) 1991-2008 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
+%
+% LaTeX Style file
+% tensorOperator.sty
+%
+% Description
+% Standard OpenCFD LaTeX macros for typesetting tensor algebra.
+%
+%------------------------------------------------------------------------------
+
+% tensor style
+% ~~~~~~~~~~~~
+\renewcommand{\vec}[1] {\ensuremath{\mathbf #1}}
+\newcommand{\gvec}[1] {\ensuremath{\mbox{\boldmath$\bf#1$}}}
+
+% products
+% ~~~~~~~~
+\newcommand{\anyprod}{\star}
+\newcommand{\cprod} {\times}
+\newcommand{\dprod} {\,{\scriptscriptstyle \stackrel{\bullet}{{}}}\,}
+\newcommand{\ddprod} {\,{\scriptscriptstyle \stackrel{\bullet}{\bullet}}\,}
+\newcommand{\tdprod} {\,{\scriptscriptstyle \stackrel{3}{\bullet}}\,}
+\newcommand{\tprod} {\,{\scriptscriptstyle \stackrel{\otimes}{{}}}\,}
+
+% operations
+% ~~~~~~~~~~
+\newcommand{\adj} {\ensuremath{\operatorname{adj}}}
+\newcommand{\cof} {\ensuremath{\operatorname{cof}}}
+\newcommand{\diag} {\ensuremath{\operatorname{diag}}}
+\newcommand{\dev} {\ensuremath{\operatorname{dev}}}
+
+\newcommand{\Hodge} {\ensuremath{\operatorname{\stackrel{\displaystyle \ast}{}}}}
+\newcommand{\hyd} {\ensuremath{\operatorname{hyd}}}
+\renewcommand{\max} {\ensuremath{\operatorname{max}}}
+\renewcommand{\min} {\ensuremath{\operatorname{min}}}
+\newcommand{\inv} {\ensuremath{\operatorname{inv}}}
+\newcommand{\sym} {\ensuremath{\operatorname{symm}}} % symm ?
+\newcommand{\skw} {\ensuremath{\operatorname{skew}}} % skew already defined
+\newcommand{\tr} {\ensuremath{\operatorname{tr}}}
+\newcommand{\trans}[1] {\ensuremath{#1^{\operatorname{T}}}}
+
+% alternative tensor operators for hypersonics etc.
+% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+\newcommand{\devs}[1] {\overset{\scriptscriptstyle\circ}{#1}}
+%\newcommand{\trans}[1] {\ensuremath{#1^{\operatorname{T}}}}
+\newcommand{\symms}[1] {\overleftrightarrow{#1}}
+\newlength{\skewslength}
+\newlength{\skewsheight}
+\newcommand{\skews}[1]{
+ \settowidth{\skewslength}{#1}%
+ \settoheight{\skewsheight}{#1}%
+ \addtolength{\skewsheight}{0.4mm}%
+ {\overleftrightarrow{#1}\hspace{-.5\skewslength}%
+ \rule[\skewsheight]{.4pt}{1.4mm}
+ \hspace{.5\skewslength}%
+}}
+%\newcommand{\skew}[1] {\ensuremath{#1^{\operatorname{A}}}}
+
+% spatial derivatives
+% ~~~~~~~~~~~~~~~~~~~
+\newcommand{\curl}{\ensuremath{\nabla\cprod}}
+\renewcommand{\div} {\ensuremath{\nabla\dprod}}
+\newcommand{\grad}{\ensuremath{\nabla}}
+\newcommand{\laplacian}{\ensuremath{\nabla^{2}}}
+
+% temporal derivatives
+% ~~~~~~~~~~~~~~~~~~~~
+\newcommand{\ddt}[1] {\ensuremath{\frac{\partial #1}{\partial t }}}
+\newcommand{\DDt}[1] {\ensuremath{\frac{D #1}{D t}}}
+\newcommand{\DpDt}[2] {\ensuremath{\frac{d_{#1} #2}{d t }}}
+\newcommand{\dsdts}[1] {\ensuremath{\frac{\partial ^2 #1}{\partial t^2}}}
+\newcommand{\rate}[1] {\ensuremath{\dot{#1}}}
+
+\newcommand{\genDer}{\mathcal{L}}
+
+% time average symbols
+% ~~~~~~~~~~~~~~~~~~~~
+\newcommand{\av}[1] {\ensuremath{\overline{#1}}}
+\newcommand{\corrtwo}[2] {{\dwea{\dprime{#1} \dprime{#2}}}}
+\newcommand{\curly}[1] {{\cal #1}}
+\newcommand{\dprime}[1] {\ensuremath{{#1}^{^{\prime \prime}}}}
+\newcommand{\dwea}[1] {\ensuremath{\widetilde{#1}}}
+\newcommand{\dweafluc}[1] {\ensuremath{\dprime{#1}}}
+\newcommand{\fluc}[1] {\ensuremath{#1^{\prime}}}
+
+% index style
+% ~~~~~~~~~~~
+\newcommand{\veci}[2][i] {\ensuremath{#2_{#1}}}
+\newcommand{\teni}[2][ij] {\ensuremath{#2_{#1}}}
+\newcommand{\tenTi}[2][ji] {\ensuremath{#2_{#1}}}
+
+% index operations
+% ~~~~~~~~~~~~~~~~
+\newcommand{\deltai}[1] {\ensuremath{\partial_{#1}}}
+
+% Sub-subscripts
+% ~~~~~~~~~~~~~~
+\newcommand{\eff} {{\scriptscriptstyle e\!f\!\!f\!}}
+
+% unknown use
+% ~~~~~~~~~~~
+%\font\bigtenrm=cmr12 scaled 1200
+%\newcommand{\eexp}[1]{{\hbox{$\textfont1=\bigtenrm e$}}^{\raise3pt
+%\hbox{$#1$}}}
+
+
+% ------------------------------------------------------------------------------
diff --git a/etc/settings.csh b/etc/settings.csh
index 2c9a0eb656..4bcf3eb510 100644
--- a/etc/settings.csh
+++ b/etc/settings.csh
@@ -187,7 +187,7 @@ case MPICH-GM:
setenv FOAM_MPI_LIBBIN $FOAM_LIBBIN/mpich-gm
breaksw
-case MPICH-GM:
+case HPMPI:
setenv MPI_HOME /opt/hpmpi
setenv MPI_ARCH_PATH $MPI_HOME
setenv MPICH_ROOT=$MPI_ARCH_PATH
diff --git a/src/OSspecific/Unix/regularExpression.H b/src/OSspecific/Unix/regularExpression.H
index 9924caef28..fed14379ad 100644
--- a/src/OSspecific/Unix/regularExpression.H
+++ b/src/OSspecific/Unix/regularExpression.H
@@ -37,6 +37,8 @@ SourceFiles
#include
#include
+#include "string.H"
+#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -73,7 +75,7 @@ public:
{
preg_ = new regex_t;
- if (regcomp(preg_, s.c_str(), REG_EXTENDED|REG_NOSUB) != 0)
+ if (regcomp(preg_, s.c_str(), REG_EXTENDED) != 0)
{
FatalErrorIn
(
@@ -102,10 +104,12 @@ public:
//- Matches?
inline bool matches(const string& s) const
{
- size_t nmatch = 0;
- regmatch_t *pmatch = NULL;
+ size_t nmatch = 1;
+ regmatch_t pmatch[1];
- return regexec(preg_, s.c_str(), nmatch, pmatch, 0) == 0;
+ int errVal = regexec(preg_, s.c_str(), nmatch, pmatch, 0);
+
+ return (errVal == 0 && pmatch[0].rm_eo == label(s.size()));
}
};
diff --git a/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.C b/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.C
index 24d1edea53..96159327c1 100644
--- a/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.C
+++ b/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.C
@@ -28,23 +28,23 @@ License
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
-// Construct from Istream
+
template
Foam::DynamicList::DynamicList(Istream& is)
:
List(is),
- allocSize_(List::size())
+ capacity_(List::size())
{}
template
Foam::Ostream& Foam::operator<<
(
- Foam::Ostream& os,
- const Foam::DynamicList& DL
+ Ostream& os,
+ const DynamicList& lst
)
{
- os << static_cast&>(DL);
+ os << static_cast&>(lst);
return os;
}
@@ -52,12 +52,12 @@ Foam::Ostream& Foam::operator<<
template
Foam::Istream& Foam::operator>>
(
- Foam::Istream& is,
- Foam::DynamicList& DL
+ Istream& is,
+ DynamicList& lst
)
{
- is >> static_cast&>(DL);
- DL.allocSize_ = DL.List::size();
+ is >> static_cast&>(lst);
+ lst.capacity_ = lst.List::size();
return is;
}
diff --git a/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H b/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H
index 214c780620..847b94f647 100644
--- a/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H
+++ b/src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H
@@ -32,7 +32,7 @@ Description
Internal storage is a compact array and the list can be shrunk to compact
storage. The increase of list size is controlled by three template
parameters, which allows the list storage to either increase by the given
- increment or the given multiplier and divider (allowing non-integer
+ increment or by the given multiplier and divider (allowing non-integer
multiples).
SourceFiles
@@ -81,14 +81,11 @@ class DynamicList
{
// Private data
- //- Allocated size for underlying List.
- label allocSize_;
+ //- The capacity (allocated size) of the underlying list.
+ label capacity_;
// Private Member Functions
- // Disabled, since the usefulness and semantics are not quite clear
- void setSize(const label, const T&);
-
public:
// Related types
@@ -116,20 +113,30 @@ public:
// Access
//- Size of the underlying storage.
- inline label allocSize() const;
-
+ inline label capacity() const;
// Edit
- //- Alter the list size.
- // When the new size is greater than the addressed list size, the
- // allocated list sizes is adjusted and the
- // addressed size does not change.
- // Otherwise the addressed list size is just reduced and the
- // allocated size does not change.
+ //- Alter the size of the underlying storage.
+ // The addressed size will be truncated if needed to fit, but will
+ // remain otherwise untouched.
+ // Use this or reserve() in combination with append().
+ inline void setCapacity(const label);
+
+ //- Alter the addressed list size.
+ // New space will be allocated if required.
+ // Use this to resize the list prior to using the operator[] for
+ // setting values (as per List usage).
inline void setSize(const label);
- //- Clear the list, i.e. set the size to zero.
+ //- Alter the addressed list size and fill new space with a constant.
+ inline void setSize(const label, const T&);
+
+ //- Reserve allocation space for at least this size.
+ // Never shrinks the allocated size, use setCapacity() for that.
+ inline void reserve(const label);
+
+ //- Clear the addressed list, i.e. set the size to zero.
// Allocated size does not change
inline void clear();
@@ -152,6 +159,9 @@ public:
//- Append an element at the end of the list
inline void append(const T& e);
+ //- Append a List at the end of this list
+ inline void append(const UList&);
+
//- Remove and return the top element
inline T remove();
@@ -162,7 +172,7 @@ public:
inline void operator=(const T&);
//- Assignment from List. Also handles assignment from DynamicList.
- inline void operator=(const List&);
+ inline void operator=(const UList&);
// IOstream operators
diff --git a/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H b/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
index 1a40bdd0ac..1f7b8d4043 100644
--- a/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
+++ b/src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
@@ -30,7 +30,7 @@ template
inline Foam::DynamicList::DynamicList()
:
List(SizeInc),
- allocSize_(SizeInc)
+ capacity_(SizeInc)
{
List::size(0);
}
@@ -39,11 +39,11 @@ inline Foam::DynamicList::DynamicList()
template
inline Foam::DynamicList::DynamicList
(
- const label s
+ const label nElem
)
:
- List(s),
- allocSize_(s)
+ List(nElem),
+ capacity_(nElem)
{
List::size(0);
}
@@ -56,39 +56,101 @@ inline Foam::DynamicList::DynamicList
)
:
List(lst),
- allocSize_(lst.size())
+ capacity_(lst.size())
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template
-inline Foam::label Foam::DynamicList::allocSize()
+inline Foam::label Foam::DynamicList::capacity()
const
{
- return allocSize_;
+ return capacity_;
+}
+
+
+template
+inline void Foam::DynamicList::setCapacity
+(
+ const label nElem
+)
+{
+ label nextFree = List::size();
+ capacity_ = nElem;
+
+ if (nextFree > capacity_)
+ {
+ // truncate addressed sizes too
+ nextFree = capacity_;
+ }
+
+ List::setSize(capacity_);
+ List::size(nextFree);
+}
+
+
+template
+inline void Foam::DynamicList::reserve
+(
+ const label nElem
+)
+{
+ // allocate more capacity?
+ if (nElem > capacity_)
+ {
+ capacity_ = max
+ (
+ nElem,
+ label(SizeInc + capacity_ * SizeMult / SizeDiv)
+ );
+
+ // adjust allocated size, leave addressed size untouched
+ label nextFree = List::size();
+ List::setSize(capacity_);
+ List::size(nextFree);
+ }
}
template
inline void Foam::DynamicList::setSize
(
- const label s
+ const label nElem
+)
+{
+ // allocate more capacity?
+ if (nElem > capacity_)
+ {
+ capacity_ = max
+ (
+ nElem,
+ label(SizeInc + capacity_ * SizeMult / SizeDiv)
+ );
+
+ List::setSize(capacity_);
+ }
+
+ // adjust addressed size
+ List::size(nElem);
+}
+
+
+template
+inline void Foam::DynamicList::setSize
+(
+ const label nElem,
+ const T& t
)
{
label nextFree = List::size();
- if (s <= nextFree)
+ setSize(nElem);
+
+ // set new elements to constant value
+ while (nextFree < nElem)
{
- // adjust addressed size, leave allocated size untouched
- nextFree = s;
+ this->operator[](nextFree++) = t;
}
- else
- {
- // adjust allocated size, leave addressed size untouched
- allocSize_ = s;
- List::setSize(allocSize_);
- }
- List::size(nextFree);
}
@@ -103,7 +165,7 @@ template
inline void Foam::DynamicList::clearStorage()
{
List::clear();
- allocSize_ = 0;
+ capacity_ = 0;
}
@@ -111,13 +173,15 @@ template
inline Foam::DynamicList&
Foam::DynamicList::shrink()
{
- if (allocSize_ > List::size())
+ label nextFree = List::size();
+ if (capacity_ > nextFree)
{
- allocSize_ = List::size();
- // force re-allocation/copying in List::setSize() by temporarily
- // faking a larger list size that will be truncated
- List::size(allocSize_+1);
- List::setSize(allocSize_);
+ // use the full list when resizing
+ List::size(capacity_);
+ // the new size
+ capacity_ = nextFree;
+ List::setSize(capacity_);
+ List::size(nextFree);
}
return *this;
}
@@ -127,7 +191,7 @@ template
inline void
Foam::DynamicList::transfer(List& lst)
{
- allocSize_ = lst.size();
+ capacity_ = lst.size();
List::transfer(lst); // take over storage, clear addressing for lst.
}
@@ -140,41 +204,58 @@ Foam::DynamicList::transfer
)
{
// take over storage as-is (without shrink), clear addressing for lst.
- allocSize_ = lst.allocSize_;
- lst.allocSize_ = 0;
+ capacity_ = lst.capacity_;
+ lst.capacity_ = 0;
List::transfer(static_cast&>(lst));
}
template
-inline void Foam::DynamicList::append(const T& e)
+inline void Foam::DynamicList::append
+(
+ const T& t
+)
+{
+ label elemI = List::size();
+ setSize(elemI + 1);
+
+ this->operator[](elemI) = t;
+}
+
+
+template
+inline void Foam::DynamicList::append
+(
+ const UList& lst
+)
{
- // Work on copy free index since gets overwritten by setSize
label nextFree = List::size();
- nextFree++;
-
- if (nextFree > allocSize_)
+ if (this == &lst)
{
- allocSize_ = max
+ FatalErrorIn
(
- nextFree,
- label(SizeMult*allocSize_/SizeDiv + SizeInc)
- );
- List::setSize(allocSize_);
+ "DynamicList::append"
+ "(const UList&)"
+ ) << "attempted appending to self" << abort(FatalError);
}
- List::size(nextFree);
+ setSize(nextFree + lst.size());
- this->operator[](nextFree - 1) = e;
+ forAll(lst, elemI)
+ {
+ this->operator[](nextFree++) = lst[elemI];
+ }
}
template
inline T Foam::DynamicList::remove()
{
- if (List::size() == 0)
+ label elemI = List::size() - 1;
+
+ if (elemI < 0)
{
FatalErrorIn
(
@@ -182,11 +263,9 @@ inline T Foam::DynamicList::remove()
) << "List is empty" << abort(FatalError);
}
- label nextFree = List::size()-1;
+ const T& val = List::operator[](elemI);
- const T& val = List::operator[](nextFree);
-
- List::size(nextFree);
+ List::size(elemI);
return val;
}
@@ -197,26 +276,15 @@ inline T Foam::DynamicList::remove()
template
inline T& Foam::DynamicList::operator()
(
- const label i
+ const label elemI
)
{
- label nextFree = List::size();
-
- nextFree = max(nextFree, i + 1);
-
- if (nextFree > allocSize_)
+ if (elemI >= List::size())
{
- allocSize_ = max
- (
- nextFree,
- label(SizeMult*allocSize_/SizeDiv + SizeInc)
- );
- List::setSize(allocSize_);
+ setSize(elemI + 1);
}
- List::size(nextFree);
-
- return this->operator[](i);
+ return this->operator[](elemI);
}
@@ -226,21 +294,39 @@ inline void Foam::DynamicList::operator=
const T& t
)
{
- List::operator=(t);
+ UList::operator=(t);
}
template
inline void Foam::DynamicList::operator=
(
- const List& lst
+ const UList& lst
)
{
- // make the entire storage available for the copy operation:
- List::size(allocSize_);
+ if (this == &lst)
+ {
+ FatalErrorIn
+ (
+ "DynamicList::operator="
+ "(const UList&)"
+ ) << "attempted assignment to self" << abort(FatalError);
+ }
- List::operator=(lst);
- allocSize_ = List::size();
+ if (capacity_ >= lst.size())
+ {
+ // can copy w/o reallocating, match initial size to avoid reallocation
+ List::size(lst.size());
+ List::operator=(lst);
+ }
+ else
+ {
+ // make everything available for the copy operation
+ List::size(capacity_);
+
+ List::operator=(lst);
+ capacity_ = List::size();
+ }
}
diff --git a/src/OpenFOAM/containers/Lists/List/List.C b/src/OpenFOAM/containers/Lists/List/List.C
index 5d6c1608f7..83f3464cca 100644
--- a/src/OpenFOAM/containers/Lists/List/List.C
+++ b/src/OpenFOAM/containers/Lists/List/List.C
@@ -432,14 +432,19 @@ void Foam::List::transfer(DynamicList& a)
{
// shrink the allocated space to the number of elements used
a.shrink();
+ transfer(static_cast&>(a));
+ a.clearStorage();
+}
- if (this->v_) delete[] this->v_;
- this->size_ = a.size_;
- this->v_ = a.v_;
- a.size_ = 0;
- a.v_ = 0;
- a.allocSize_ = 0;
+// Transfer the contents of the argument SortableList into this List
+// and anull the argument list
+template
+void Foam::List::transfer(SortableList& a)
+{
+ // shrink away the sort indices
+ a.shrink();
+ transfer(static_cast&>(a));
}
diff --git a/src/OpenFOAM/containers/Lists/List/List.H b/src/OpenFOAM/containers/Lists/List/List.H
index 4f51878d45..c3fada1d62 100644
--- a/src/OpenFOAM/containers/Lists/List/List.H
+++ b/src/OpenFOAM/containers/Lists/List/List.H
@@ -63,7 +63,8 @@ template class FixedList;
template class PtrList;
template class SLList;
template
- class DynamicList;
+ class DynamicList;
+template class SortableList;
template class IndirectList;
template class BiIndirectList;
@@ -173,6 +174,10 @@ public:
template
void transfer(DynamicList&);
+ //- Transfer the contents of the argument List into this List
+ // and annull the argument list.
+ void transfer(SortableList&);
+
//- Return subscript-checked element of UList.
inline T& newElmt(const label);
diff --git a/src/OpenFOAM/containers/Lists/ListOps/ListOps.H b/src/OpenFOAM/containers/Lists/ListOps/ListOps.H
index 889dd57d7e..67135f957c 100644
--- a/src/OpenFOAM/containers/Lists/ListOps/ListOps.H
+++ b/src/OpenFOAM/containers/Lists/ListOps/ListOps.H
@@ -82,6 +82,13 @@ void inplaceMapKey(const UList& oldToNew, Container&);
template
void sortedOrder(const UList&, labelList& order);
+//- Generate (sorted) indices corresponding to duplicate list values
+template
+void duplicateOrder(const UList&, labelList& order);
+
+//- Generate (sorted) indices corresponding to unique list values
+template
+void uniqueOrder(const UList&, labelList& order);
//- Extract elements of List whose region is certain value.
// Use e.g. to extract all selected elements:
diff --git a/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C b/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C
index b048ef1bb8..942cb6e8c0 100644
--- a/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C
+++ b/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C
@@ -174,16 +174,65 @@ void Foam::sortedOrder
)
{
order.setSize(lst.size());
-
forAll(order, elemI)
{
order[elemI] = elemI;
}
-
Foam::stableSort(order, typename UList::less(lst));
}
+template
+void Foam::duplicateOrder
+(
+ const UList& lst,
+ labelList& order
+)
+{
+ if (lst.size() < 2)
+ {
+ order.clear();
+ return;
+ }
+
+ sortedOrder(lst, order);
+
+ label n = 0;
+ for (label i = 0; i < order.size() - 1; ++i)
+ {
+ if (lst[order[i]] == lst[order[i+1]])
+ {
+ order[n++] = order[i];
+ }
+ }
+ order.setSize(n);
+}
+
+
+template
+void Foam::uniqueOrder
+(
+ const UList& lst,
+ labelList& order
+)
+{
+ sortedOrder(lst, order);
+
+ if (order.size() > 1)
+ {
+ label n = 0;
+ for (label i = 0; i < order.size() - 1; ++i)
+ {
+ if (lst[order[i]] != lst[order[i+1]])
+ {
+ order[n++] = order[i];
+ }
+ }
+ order.setSize(n);
+ }
+}
+
+
template
ListType Foam::subset
(
diff --git a/src/OpenFOAM/containers/Lists/SortableList/SortableList.C b/src/OpenFOAM/containers/Lists/SortableList/SortableList.C
index cb6c1571d3..0ae74df8af 100644
--- a/src/OpenFOAM/containers/Lists/SortableList/SortableList.C
+++ b/src/OpenFOAM/containers/Lists/SortableList/SortableList.C
@@ -24,6 +24,23 @@ License
\*---------------------------------------------------------------------------*/
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+template
+void Foam::SortableList::sortIndices(List& ind) const
+{
+ // list lengths must be identical
+ ind.setSize(this->size());
+
+ forAll(ind, i)
+ {
+ ind[i] = i;
+ }
+
+ Foam::stableSort(ind, typename UList::less(*this));
+}
+
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template
@@ -73,13 +90,6 @@ Foam::SortableList::SortableList(const SortableList& lst)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-template
-void Foam::SortableList::setSize(const label newSize)
-{
- List::setSize(newSize);
- indices_.setSize(newSize, -1);
-}
-
template
void Foam::SortableList::clear()
@@ -100,15 +110,7 @@ Foam::List& Foam::SortableList::shrink()
template
void Foam::SortableList::sort()
{
- // list lengths must be identical
- indices_.setSize(this->size());
-
- forAll(indices_, i)
- {
- indices_[i] = i;
- }
-
- Foam::stableSort(indices_, typename UList::less(*this));
+ sortIndices(indices_);
List lst(this->size());
forAll(indices_, i)
@@ -120,10 +122,33 @@ void Foam::SortableList::sort()
}
+template
+void Foam::SortableList::reverseSort()
+{
+ sortIndices(indices_);
+
+ List lst(this->size());
+ label endI = indices_.size();
+ forAll(indices_, i)
+ {
+ lst[--endI] = this->operator[](indices_[i]);
+ }
+
+ List::transfer(lst);
+}
+
+
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template
-void Foam::SortableList::operator=(const UList& rhs)
+inline void Foam::SortableList::operator=(const Type& t)
+{
+ UList::operator=(t);
+}
+
+
+template
+inline void Foam::SortableList::operator=(const UList& rhs)
{
List::operator=(rhs);
indices_.clear();
@@ -131,7 +156,7 @@ void Foam::SortableList::operator=(const UList& rhs)
template
-void Foam::SortableList::operator=(const SortableList& rhs)
+inline void Foam::SortableList::operator=(const SortableList& rhs)
{
List::operator=(rhs);
indices_ = rhs.indices();
diff --git a/src/OpenFOAM/containers/Lists/SortableList/SortableList.H b/src/OpenFOAM/containers/Lists/SortableList/SortableList.H
index 01cdce3780..40dc6dd68c 100644
--- a/src/OpenFOAM/containers/Lists/SortableList/SortableList.H
+++ b/src/OpenFOAM/containers/Lists/SortableList/SortableList.H
@@ -60,6 +60,9 @@ class SortableList
//- Original indices
labelList indices_;
+ //- Resize and sort the parameter according to the list values
+ void sortIndices(List&) const;
+
public:
// Constructors
@@ -99,9 +102,6 @@ public:
return indices_;
}
- //- Size the list. Growing can cause undefined indices (until next sort)
- void setSize(const label);
-
//- Clear the list and the indices
void clear();
@@ -112,13 +112,19 @@ public:
// also resizes the indices as required
void sort();
+ //- Reverse (stable) sort the list
+ void reverseSort();
+
// Member Operators
+ //- Assignment of all entries to the given value
+ inline void operator=(const Type&);
+
//- Assignment from UList operator. Takes linear time.
- void operator=(const UList&);
+ inline void operator=(const UList&);
//- Assignment operator. Takes linear time.
- void operator=(const SortableList&);
+ inline void operator=(const SortableList&);
};
diff --git a/src/OpenFOAM/containers/Lists/SubList/SubList.H b/src/OpenFOAM/containers/Lists/SubList/SubList.H
index dd5bcc673e..0e8336dc50 100644
--- a/src/OpenFOAM/containers/Lists/SubList/SubList.H
+++ b/src/OpenFOAM/containers/Lists/SubList/SubList.H
@@ -87,6 +87,9 @@ public:
//- Allow cast to a const List&
inline operator const Foam::List&() const;
+
+ //- Assignment of all entries to the given value
+ inline void operator=(const T&);
};
diff --git a/src/OpenFOAM/containers/Lists/SubList/SubListI.H b/src/OpenFOAM/containers/Lists/SubList/SubListI.H
index 6b142cab1f..7a112370fb 100644
--- a/src/OpenFOAM/containers/Lists/SubList/SubListI.H
+++ b/src/OpenFOAM/containers/Lists/SubList/SubListI.H
@@ -54,7 +54,7 @@ inline Foam::SubList::SubList
# ifdef FULLDEBUG
// Artificially allowing the start of a zero-sized subList to be
- // one past the end of the original list.
+ // one past the end of the original list.
if (subSize > 0)
{
list.checkStart(startIndex);
@@ -73,7 +73,7 @@ inline Foam::SubList::SubList
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template
-const Foam::SubList& Foam::SubList::null()
+inline const Foam::SubList