diff --git a/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C b/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C
index 3e9f8cd159e1bac2283dd32447c55ec850e1de70..e0dda32b2a745a56c3ada9931dfac57fd0d7dd30 100644
--- a/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C
+++ b/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C
@@ -46,8 +46,8 @@ int main(int argc, char *argv[])
);
argList::noBanner();
- argList::noParallel();
argList::addBoolOption("times", "list available times");
+ argList::addBoolOption("latestTime", "list last time");
argList::addBoolOption
(
"keywords",
@@ -75,6 +75,15 @@ int main(int argc, char *argv[])
Info<< times[i].name() << endl;
}
}
+ else if (args.optionFound("latestTime"))
+ {
+ instantList times
+ (
+ Foam::Time::findTimes(args.rootPath()/args.caseName())
+ );
+
+ Info<< times.last().name() << endl;
+ }
if (args.optionFound("dict"))
{
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
index 21ba1cdfa1c7bbb13dc7d3cdc05710323cb8afd6..b569601734da365092cb2c6937b279fffb1d45ec 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
@@ -99,6 +99,40 @@ Usage
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+const labelIOList& procAddressing
+(
+ const PtrList<fvMesh>& procMeshList,
+ const label procI,
+ const word& name,
+ PtrList<labelIOList>& procAddressingList
+)
+{
+ const fvMesh& procMesh = procMeshList[procI];
+
+ if (!procAddressingList.set(procI))
+ {
+ procAddressingList.set
+ (
+ procI,
+ new labelIOList
+ (
+ IOobject
+ (
+ name,
+ procMesh.facesInstance(),
+ procMesh.meshSubDir,
+ procMesh,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ )
+ )
+ );
+ }
+ return procAddressingList[procI];
+}
+
+
+
int main(int argc, char *argv[])
{
argList::addNote
@@ -805,74 +839,29 @@ int main(int argc, char *argv[])
}
const fvMesh& procMesh = procMeshList[procI];
+ const labelIOList& faceProcAddressing = procAddressing
+ (
+ procMeshList,
+ procI,
+ "faceProcAddressing",
+ faceProcAddressingList
+ );
- if (!faceProcAddressingList.set(procI))
- {
- faceProcAddressingList.set
- (
- procI,
- new labelIOList
- (
- IOobject
- (
- "faceProcAddressing",
- procMesh.facesInstance(),
- procMesh.meshSubDir,
- procMesh,
- IOobject::MUST_READ,
- IOobject::NO_WRITE
- )
- )
- );
- }
- const labelIOList& faceProcAddressing =
- faceProcAddressingList[procI];
-
-
- if (!cellProcAddressingList.set(procI))
- {
- cellProcAddressingList.set
- (
- procI,
- new labelIOList
- (
- IOobject
- (
- "cellProcAddressing",
- procMesh.facesInstance(),
- procMesh.meshSubDir,
- procMesh,
- IOobject::MUST_READ,
- IOobject::NO_WRITE
- )
- )
- );
- }
- const labelIOList& cellProcAddressing =
- cellProcAddressingList[procI];
-
+ const labelIOList& cellProcAddressing = procAddressing
+ (
+ procMeshList,
+ procI,
+ "cellProcAddressing",
+ cellProcAddressingList
+ );
- if (!boundaryProcAddressingList.set(procI))
- {
- boundaryProcAddressingList.set
- (
- procI,
- new labelIOList
- (
- IOobject
- (
- "boundaryProcAddressing",
- procMesh.facesInstance(),
- procMesh.meshSubDir,
- procMesh,
- IOobject::MUST_READ,
- IOobject::NO_WRITE
- )
- )
- );
- }
- const labelIOList& boundaryProcAddressing =
- boundaryProcAddressingList[procI];
+ const labelIOList& boundaryProcAddressing = procAddressing
+ (
+ procMeshList,
+ procI,
+ "boundaryProcAddressing",
+ boundaryProcAddressingList
+ );
// FV fields
@@ -959,27 +948,13 @@ int main(int argc, char *argv[])
|| pointTensorFields.size()
)
{
- if (!pointProcAddressingList.set(procI))
- {
- pointProcAddressingList.set
- (
- procI,
- new labelIOList
- (
- IOobject
- (
- "pointProcAddressing",
- procMesh.facesInstance(),
- procMesh.meshSubDir,
- procMesh,
- IOobject::MUST_READ,
- IOobject::NO_WRITE
- )
- )
- );
- }
- const labelIOList& pointProcAddressing =
- pointProcAddressingList[procI];
+ const labelIOList& pointProcAddressing = procAddressing
+ (
+ procMeshList,
+ procI,
+ "pointProcAddressing",
+ pointProcAddressingList
+ );
const pointMesh& procPMesh = pointMesh::New(procMesh);
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
index 4eda411482090d25f72fd88f4f2c209d0f6a5fb7..923d9ae1cf402f956c32c3d9c8c5d5244b1978b4 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
@@ -38,6 +38,7 @@ License
#include "cellSet.H"
#include "faceSet.H"
#include "pointSet.H"
+#include "uniformDimensionedFields.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -195,6 +196,60 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
}
+ autoPtr<labelIOList> cellLevelPtr;
+ {
+ IOobject io
+ (
+ "cellLevel",
+ facesInstance(),
+ polyMesh::meshSubDir,
+ *this,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ );
+ if (io.headerOk())
+ {
+ Info<< "Reading hexRef8 data : " << io.name() << endl;
+ cellLevelPtr.reset(new labelIOList(io));
+ }
+ }
+ autoPtr<labelIOList> pointLevelPtr;
+ {
+ IOobject io
+ (
+ "pointLevel",
+ facesInstance(),
+ polyMesh::meshSubDir,
+ *this,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ );
+ if (io.headerOk())
+ {
+ Info<< "Reading hexRef8 data : " << io.name() << endl;
+ pointLevelPtr.reset(new labelIOList(io));
+ }
+ }
+ autoPtr<uniformDimensionedScalarField> level0EdgePtr;
+ {
+ IOobject io
+ (
+ "level0Edge",
+ facesInstance(),
+ polyMesh::meshSubDir,
+ *this,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ );
+ if (io.headerOk())
+ {
+ Info<< "Reading hexRef8 data : " << io.name() << endl;
+ level0EdgePtr.reset(new uniformDimensionedScalarField(io));
+ }
+ }
+
+
+
label maxProcCells = 0;
label totProcFaces = 0;
label maxProcPatches = 0;
@@ -767,8 +822,28 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
procMesh.write();
+ // Write points if pointsInstance differing from facesInstance
+ if (facesInstancePointsPtr_.valid())
+ {
+ pointIOField pointsInstancePoints
+ (
+ IOobject
+ (
+ "points",
+ pointsInstance(),
+ polyMesh::meshSubDir,
+ procMesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE,
+ false
+ ),
+ xferMove(procPoints)
+ );
+ pointsInstancePoints.write();
+ }
+ // Decompose any sets
if (decomposeSets)
{
forAll(cellSets, i)
@@ -813,25 +888,67 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
}
- // Write points if pointsInstance differing from facesInstance
- if (facesInstancePointsPtr_.valid())
+ // hexRef8 data
+ if (cellLevelPtr.valid())
{
- pointIOField pointsInstancePoints
+ labelIOList
(
IOobject
(
- "points",
- pointsInstance(),
+ cellLevelPtr().name(),
+ facesInstance(),
polyMesh::meshSubDir,
procMesh,
IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
+ IOobject::AUTO_WRITE
),
- xferMove(procPoints)
- );
- pointsInstancePoints.write();
+ UIndirectList<label>
+ (
+ cellLevelPtr(),
+ procCellAddressing_[procI]
+ )()
+ ).write();
}
+ if (pointLevelPtr.valid())
+ {
+ labelIOList
+ (
+ IOobject
+ (
+ pointLevelPtr().name(),
+ facesInstance(),
+ polyMesh::meshSubDir,
+ procMesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ UIndirectList<label>
+ (
+ pointLevelPtr(),
+ procPointAddressing_[procI]
+ )()
+ ).write();
+ }
+ if (level0EdgePtr.valid())
+ {
+ uniformDimensionedScalarField
+ (
+ IOobject
+ (
+ level0EdgePtr().name(),
+ facesInstance(),
+ polyMesh::meshSubDir,
+ procMesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ level0EdgePtr()
+ ).write();
+ }
+
+
+
+ // Statistics
Info<< endl
<< "Processor " << procI << nl
diff --git a/src/OpenFOAM/fields/pointPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValuePointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValuePointPatchField.C
index 01cfc4b7bd1da8091781fbe44eb1420161e02459..2590a08dcdae62668dcfb7ef84c8d736edb1fbe6 100644
--- a/src/OpenFOAM/fields/pointPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValuePointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValuePointPatchField.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -55,7 +55,9 @@ timeVaryingUniformFixedValuePointPatchField
:
fixedValuePointPatchField<Type>(ptf, p, iF, mapper),
timeSeries_(ptf.timeSeries_)
-{}
+{
+ updateCoeffs();
+}
template<class Type>
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
index fc2b484cb0a4a478eb2a1591a8f504c5afd74eea..7c94b9b1380af34a51b40c347de9306e64d2a78d 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
@@ -251,27 +251,35 @@ public:
{
public:
- template<class T>
+ template<class Type>
void operator()
(
const vectorTensorTransform& vt,
const bool forward,
- List<T>& fld
+ List<Type>& fld
) const
{
- if (forward)
- {
- transformList(vt.R(), fld);
- }
- else
+ const tensor T(forward ? vt.R() : vt.R().T());
+ transformList(T, fld);
+ }
+
+ template<class Type>
+ void operator()
+ (
+ const vectorTensorTransform& vt,
+ const bool forward,
+ List<List<Type> >& flds
+ ) const
+ {
+ forAll(flds, i)
{
- transformList(vt.R().T(), fld);
+ operator()(vt, forward, flds[i]);
}
}
//- Transform patch-based field
- template<class T>
- void operator()(const coupledPolyPatch& cpp, UList<T>& fld) const
+ template<class Type>
+ void operator()(const coupledPolyPatch& cpp, UList<Type>& fld) const
{
if (!cpp.parallel())
{
@@ -280,8 +288,8 @@ public:
}
//- Transform sparse field
- template<class T, template<class> class Container>
- void operator()(const coupledPolyPatch& cpp, Container<T>& map)
+ template<class Type, template<class> class Container>
+ void operator()(const coupledPolyPatch& cpp, Container<Type>& map)
const
{
if (!cpp.parallel())
@@ -313,6 +321,18 @@ public:
fld = vt.invTransformPosition(pfld);
}
}
+ void operator()
+ (
+ const vectorTensorTransform& vt,
+ const bool forward,
+ List<List<point> >& flds
+ ) const
+ {
+ forAll(flds, i)
+ {
+ operator()(vt, forward, flds[i]);
+ }
+ }
//- Transform patch-based field
void operator()(const coupledPolyPatch& cpp, pointField& fld) const
{
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C
index a5cb7fdbb982c4a7d8eae8d1be45bf09d378e8a0..8ac95cd8b1361ee90968007394fedeaa1a5f9034 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C
@@ -28,43 +28,6 @@ License
#include "indirectPrimitivePatch.H"
#include "globalMeshData.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
- //- Transformation
- class listTransform
- {
- public:
-
- void operator()
- (
- const vectorTensorTransform& vt,
- const bool forward,
- List<List<point> >& fld
- ) const
- {
- const tensor T
- (
- forward
- ? vt.R()
- : vt.R().T()
- );
-
- forAll(fld, i)
- {
- List<point>& elems = fld[i];
- forAll(elems, elemI)
- {
- elems[elemI] = transform(T, elems[elemI]);
- }
- }
- }
- };
-}
-
-
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template
@@ -141,7 +104,7 @@ Foam::PatchTools::pointNormals
(
transforms,
pointFaceNormals,
- listTransform()
+ mapDistribute::transform()
);
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
index b628afde9a6ed642b2d5772df9503a27dcf3e318..1fb4e2452a49890b27887bdce8b87e699a5e005c 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
@@ -219,7 +219,7 @@ Foam::label Foam::removePoints::countPointUsage
pointCanBeDeleted.setSize(mesh_.nPoints());
pointCanBeDeleted = false;
- label nDeleted = 0;
+ //label nDeleted = 0;
forAll(edge0, pointI)
{
@@ -243,14 +243,14 @@ Foam::label Foam::removePoints::countPointUsage
if ((e0Vec & e1Vec) > minCos)
{
pointCanBeDeleted[pointI] = true;
- nDeleted++;
+ //nDeleted++;
}
}
else if (edge0[pointI] == -1)
{
// point not used at all
pointCanBeDeleted[pointI] = true;
- nDeleted++;
+ //nDeleted++;
}
}
edge0.clear();
@@ -300,6 +300,15 @@ Foam::label Foam::removePoints::countPointUsage
true // null value
);
+ label nDeleted = 0;
+ forAll(pointCanBeDeleted, pointI)
+ {
+ if (pointCanBeDeleted[pointI])
+ {
+ nDeleted++;
+ }
+ }
+
return returnReduce(nDeleted, sumOp<label>());
}
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 3c56fa999b000d1cf4551c648d7b5548ee488a97..07ff697df6a010fed0e6e41bcd3eeccf3f81f971 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -33,6 +33,7 @@ Description
#include "fvMesh.H"
#include "Time.H"
#include "OFstream.H"
+#include "OBJstream.H"
#include "mapPolyMesh.H"
#include "pointEdgePoint.H"
#include "PointEdgeWave.H"
@@ -1193,12 +1194,29 @@ void Foam::autoSnapDriver::detectNearSurfaces
const pointField avgCc(avgCellCentres(mesh, pp));
- // Construct rays from localPoints to beyond cell centre
+ // Construct rays through localPoints to beyond cell centre
+ pointField start(pp.nPoints());
pointField end(pp.nPoints());
forAll(localPoints, pointI)
{
const point& pt = localPoints[pointI];
- end[pointI] = pt + 2*(avgCc[pointI]-pt);
+ const vector d = 2*(avgCc[pointI]-pt);
+ start[pointI] = pt - d;
+ end[pointI] = pt + d;
+ }
+
+
+ autoPtr<OBJstream> gapStr;
+ if (debug&meshRefinement::OBJINTERSECTIONS)
+ {
+ gapStr.reset
+ (
+ new OBJstream
+ (
+ mesh.time().path()
+ / "detectNearSurfaces_" + meshRefiner_.timeName() + ".obj"
+ )
+ );
}
@@ -1226,7 +1244,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
surfaces.findNearestIntersection
(
unzonedSurfaces,
- localPoints,
+ start,
end,
surface1,
@@ -1248,38 +1266,67 @@ void Foam::autoSnapDriver::detectNearSurfaces
bool override = false;
- if (hit1[pointI].hit())
+ //if (hit1[pointI].hit())
+ //{
+ // if
+ // (
+ // meshRefiner_.isGap
+ // (
+ // planarCos,
+ // nearestPoint[pointI],
+ // nearestNormal[pointI],
+ // hit1[pointI].hitPoint(),
+ // normal1[pointI]
+ // )
+ // )
+ // {
+ // disp[pointI] = hit1[pointI].hitPoint()-pt;
+ // override = true;
+ // }
+ //}
+ //if (hit2[pointI].hit())
+ //{
+ // if
+ // (
+ // meshRefiner_.isGap
+ // (
+ // planarCos,
+ // nearestPoint[pointI],
+ // nearestNormal[pointI],
+ // hit2[pointI].hitPoint(),
+ // normal2[pointI]
+ // )
+ // )
+ // {
+ // disp[pointI] = hit2[pointI].hitPoint()-pt;
+ // override = true;
+ // }
+ //}
+
+ if (hit1[pointI].hit() && hit2[pointI].hit())
{
if
(
meshRefiner_.isGap
(
planarCos,
- nearestPoint[pointI],
- nearestNormal[pointI],
hit1[pointI].hitPoint(),
- normal1[pointI]
- )
- )
- {
- disp[pointI] = hit1[pointI].hitPoint()-pt;
- override = true;
- }
- }
- if (hit2[pointI].hit())
- {
- if
- (
- meshRefiner_.isGap
- (
- planarCos,
- nearestPoint[pointI],
- nearestNormal[pointI],
+ normal1[pointI],
hit2[pointI].hitPoint(),
normal2[pointI]
)
)
{
+ // TBD: check if the attraction (to nearest) would attract
+ // good enough and not override attraction
+
+ if (gapStr.valid())
+ {
+ const point& intPt = hit2[pointI].hitPoint();
+ gapStr().write(linePointRef(pt, intPt));
+ }
+
+ // Choose hit2 : nearest to end point (so inside the domain)
disp[pointI] = hit2[pointI].hitPoint()-pt;
override = true;
}
@@ -1337,7 +1384,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
surfaces.findNearestIntersection
(
surfacesToTest,
- pointField(localPoints, zonePointIndices),
+ pointField(start, zonePointIndices),
pointField(end, zonePointIndices),
surface1,
@@ -1363,38 +1410,63 @@ void Foam::autoSnapDriver::detectNearSurfaces
bool override = false;
- if (hit1[i].hit())
+ //if (hit1[i].hit())
+ //{
+ // if
+ // (
+ // meshRefiner_.isGap
+ // (
+ // planarCos,
+ // nearestPoint[pointI],
+ // nearestNormal[pointI],
+ // hit1[i].hitPoint(),
+ // normal1[i]
+ // )
+ // )
+ // {
+ // disp[pointI] = hit1[i].hitPoint()-pt;
+ // override = true;
+ // }
+ //}
+ //if (hit2[i].hit())
+ //{
+ // if
+ // (
+ // meshRefiner_.isGap
+ // (
+ // planarCos,
+ // nearestPoint[pointI],
+ // nearestNormal[pointI],
+ // hit2[i].hitPoint(),
+ // normal2[i]
+ // )
+ // )
+ // {
+ // disp[pointI] = hit2[i].hitPoint()-pt;
+ // override = true;
+ // }
+ //}
+
+ if (hit1[i].hit() && hit2[i].hit())
{
if
(
meshRefiner_.isGap
(
planarCos,
- nearestPoint[pointI],
- nearestNormal[pointI],
hit1[i].hitPoint(),
- normal1[i]
- )
- )
- {
- disp[pointI] = hit1[i].hitPoint()-pt;
- override = true;
- }
- }
- if (hit2[i].hit())
- {
- if
- (
- meshRefiner_.isGap
- (
- planarCos,
- nearestPoint[pointI],
- nearestNormal[pointI],
+ normal1[i],
hit2[i].hitPoint(),
normal2[i]
)
)
{
+ if (gapStr.valid())
+ {
+ const point& intPt = hit2[i].hitPoint();
+ gapStr().write(linePointRef(pt, intPt));
+ }
+
disp[pointI] = hit2[i].hitPoint()-pt;
override = true;
}
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 86e129a3e92cd3ed90426216517f7a4b198a7472..2406b8431021fed46c8b0a747966ef32bf1e0f50 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -623,6 +623,149 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
}
+// Split faces
+Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitFaces
+(
+ const labelList& splitFaces,
+ const labelPairList& splits
+)
+{
+ polyTopoChange meshMod(mesh_);
+
+ forAll(splitFaces, i)
+ {
+ label faceI = splitFaces[i];
+ const face& f = mesh_.faces()[faceI];
+
+
+ // Split as start and end index in face
+ const labelPair& split = splits[i];
+
+ label nVerts = split[1]-split[0];
+ if (nVerts < 0)
+ {
+ nVerts += f.size();
+ }
+ nVerts += 1;
+
+
+ // Split into f0, f1
+ face f0(nVerts);
+
+ label fp = split[0];
+ forAll(f0, i)
+ {
+ f0[i] = f[fp];
+ fp = f.fcIndex(fp);
+ }
+
+ face f1(f.size()-f0.size()+2);
+ fp = split[1];
+ forAll(f1, i)
+ {
+ f1[i] = f[fp];
+ fp = f.fcIndex(fp);
+ }
+
+
+ // Determine face properties
+ label own = mesh_.faceOwner()[faceI];
+ label nei = -1;
+ label patchI = -1;
+ if (faceI >= mesh_.nInternalFaces())
+ {
+ patchI = mesh_.boundaryMesh().whichPatch(faceI);
+ }
+ else
+ {
+ nei = mesh_.faceNeighbour()[faceI];
+ }
+
+ label zoneI = mesh_.faceZones().whichZone(faceI);
+ bool zoneFlip = false;
+ if (zoneI != -1)
+ {
+ const faceZone& fz = mesh_.faceZones()[zoneI];
+ zoneFlip = fz.flipMap()[fz.whichFace(faceI)];
+ }
+
+
+Pout<< "face:" << faceI << " verts:" << f
+ << " split into f0:" << f0
+ << " f1:" << f1 << endl;
+
+ // Change/add faces
+ meshMod.modifyFace
+ (
+ f0, // modified face
+ faceI, // label of face
+ own, // owner
+ nei, // neighbour
+ false, // face flip
+ patchI, // patch for face
+ zoneI, // zone for face
+ zoneFlip // face flip in zone
+ );
+ meshMod.addFace
+ (
+ f1, // modified face
+ own, // owner
+ nei, // neighbour
+ -1, // master point
+ -1, // master edge
+ faceI, // master face
+ false, // face flip
+ patchI, // patch for face
+ zoneI, // zone for face
+ zoneFlip // face flip in zone
+ );
+ }
+
+
+
+ // Change the mesh (no inflation)
+ autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
+
+ // Update fields
+ mesh_.updateMesh(map);
+
+ // Move mesh (since morphing might not do this)
+ if (map().hasMotionPoints())
+ {
+ mesh_.movePoints(map().preMotionPoints());
+ }
+ else
+ {
+ // Delete mesh volumes. No other way to do this?
+ mesh_.clearOut();
+ }
+
+ // Reset the instance for if in overwrite mode
+ mesh_.setInstance(timeName());
+ setInstance(mesh_.facesInstance());
+
+ // Update local mesh data
+ const labelList& oldToNew = map().reverseFaceMap();
+ labelList newSplitFaces(renumber(oldToNew, splitFaces));
+ // Add added faces (every splitFaces becomes two faces)
+ label sz = newSplitFaces.size();
+ newSplitFaces.setSize(2*sz);
+ forAll(map().faceMap(), faceI)
+ {
+ label oldFaceI = map().faceMap()[faceI];
+ if (oldToNew[oldFaceI] != faceI)
+ {
+ // Added face
+ newSplitFaces[sz++] = faceI;
+ }
+ }
+
+ updateMesh(map, newSplitFaces);
+
+ return map;
+}
+
+
//// Determine for multi-processor regions the lowest numbered cell on the
//// lowest numbered processor.
//void Foam::meshRefinement::getCoupledRegionMaster
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index c149378729ace764bbf924d2c166e55239a7342c..c37e9c6f833827502e393000d1595b8a050bde0a 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -866,6 +866,13 @@ public:
const point& keepPoint
);
+ //- Split faces into two
+ autoPtr<mapPolyMesh> splitFaces
+ (
+ const labelList& splitFaces,
+ const labelPairList& splits
+ );
+
//- Update local numbering for mesh redistribution
void distribute(const mapDistributePolyMesh&);
@@ -1004,6 +1011,26 @@ public:
//- Do any one of above IO functions. flag is combination of
// writeFlag values.
void write(const label flag, const fileName&) const;
+
+ //- Helper: calculate average
+ template<class T>
+ static T gAverage
+ (
+ const polyMesh& mesh,
+ const PackedBoolList& isMasterElem,
+ const UList<T>& values
+ );
+
+ //- Helper: calculate average over selected elements
+ template<class T>
+ static T gAverage
+ (
+ const polyMesh& mesh,
+ const PackedBoolList& isMasterElem,
+ const labelList& meshPoints,
+ const UList<T>& values
+ );
+
};
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index db515c0bdcb9bdd2f87d49358621d9e98e806c84..7444493e2982bfe088d20f804db3facc52b5c3b9 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -267,7 +267,12 @@ void Foam::meshRefinement::markFeatureCellLevel
// Find all start cells of features. Is done by tracking from keepPoint.
- Cloud<trackedParticle> cloud(mesh_, IDLList<trackedParticle>());
+ Cloud<trackedParticle> startPointCloud
+ (
+ mesh_,
+ "startPointCloud",
+ IDLList<trackedParticle>()
+ );
// Features are identical on all processors. Number them so we know
@@ -315,7 +320,7 @@ void Foam::meshRefinement::markFeatureCellLevel
}
// Non-manifold point. Create particle.
- cloud.addParticle
+ startPointCloud.addParticle
(
new trackedParticle
(
@@ -359,7 +364,7 @@ void Foam::meshRefinement::markFeatureCellLevel
}
// Non-manifold point. Create particle.
- cloud.addParticle
+ startPointCloud.addParticle
(
new trackedParticle
(
@@ -384,12 +389,13 @@ void Foam::meshRefinement::markFeatureCellLevel
maxFeatureLevel = labelList(mesh_.nCells(), -1);
// Database to pass into trackedParticle::move
- trackedParticle::trackingData td(cloud, maxFeatureLevel);
+ trackedParticle::trackingData td(startPointCloud, maxFeatureLevel);
// Track all particles to their end position (= starting feature point)
// Note that the particle might have started on a different processor
// so this will transfer across nicely until we can start tracking proper.
- cloud.move(td, GREAT);
+ startPointCloud.move(td, GREAT);
+
// Reset level
maxFeatureLevel = -1;
@@ -403,8 +409,64 @@ void Foam::meshRefinement::markFeatureCellLevel
featureEdgeVisited[featI] = 0u;
}
+
+ Cloud<trackedParticle> cloud
+ (
+ mesh_,
+ "featureCloud",
+ IDLList<trackedParticle>()
+ );
+
+ forAllIter(Cloud<trackedParticle>, startPointCloud, iter)
+ {
+ const trackedParticle& startTp = iter();
+
+ label featI = startTp.i();
+ label pointI = startTp.j();
+
+ const featureEdgeMesh& featureMesh = features_[featI];
+ const labelList& pEdges = featureMesh.pointEdges()[pointI];
+
+ // Now shoot particles down all pEdges.
+ forAll(pEdges, pEdgeI)
+ {
+ label edgeI = pEdges[pEdgeI];
+
+ if (featureEdgeVisited[featI].set(edgeI, 1u))
+ {
+ // Unvisited edge. Make the particle go to the other point
+ // on the edge.
+
+ const edge& e = featureMesh.edges()[edgeI];
+ label otherPointI = e.otherVertex(pointI);
+
+ trackedParticle* tp(new trackedParticle(startTp));
+ tp->end() = featureMesh.points()[otherPointI];
+ tp->j() = otherPointI;
+
+ if (debug&meshRefinement::FEATURESEEDS)
+ {
+ Pout<< "Adding particle for point:" << pointI
+ << " coord:" << tp->position()
+ << " feature:" << featI
+ << " to track to:" << tp->end()
+ << endl;
+ }
+
+ cloud.addParticle(tp);
+ }
+ }
+ }
+
+ startPointCloud.clear();
+
+
while (true)
{
+ // Track all particles to their end position.
+ cloud.move(td, GREAT);
+
+
label nParticles = 0;
// Make particle follow edge.
@@ -460,10 +522,23 @@ void Foam::meshRefinement::markFeatureCellLevel
{
break;
}
-
- // Track all particles to their end position.
- cloud.move(td, GREAT);
}
+
+
+
+ //if (debug&meshRefinement::FEATURESEEDS)
+ //{
+ // forAll(maxFeatureLevel, cellI)
+ // {
+ // if (maxFeatureLevel[cellI] != -1)
+ // {
+ // Pout<< "Feature went through cell:" << cellI
+ // << " coord:" << mesh_.cellCentres()[cellI]
+ // << " leve:" << maxFeatureLevel[cellI]
+ // << endl;
+ // }
+ // }
+ //}
}
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
index beeb47ca3cbfb0e8eb97ac18454f0b86050aa86d..812f313abe402b2b8e080d7a2b6c985b0f8fe559 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
@@ -57,6 +57,108 @@ template<class T> void meshRefinement::updateList
}
+template<class T>
+T meshRefinement::gAverage
+(
+ const polyMesh& mesh,
+ const PackedBoolList& isMasterElem,
+ const UList<T>& values
+)
+{
+ if (values.size() != isMasterElem.size())
+ {
+ FatalErrorIn
+ (
+ "meshRefinement::gAverage\n"
+ "(\n"
+ " const polyMesh&,\n"
+ " const PackedBoolList& isMasterElem,\n"
+ " const UList<T>& values\n"
+ ")\n"
+ ) << "Number of elements in list " << values.size()
+ << " does not correspond to number of elements in isMasterElem "
+ << isMasterElem.size()
+ << exit(FatalError);
+ }
+
+ T sum = pTraits<T>::zero;
+ label n = 0;
+
+ forAll(values, i)
+ {
+ if (isMasterElem[i])
+ {
+ sum += values[i];
+ n++;
+ }
+ }
+
+ reduce(sum, sumOp<T>());
+ reduce(n, sumOp<label>());
+
+ if (n > 0)
+ {
+ return sum/n;
+ }
+ else
+ {
+ return pTraits<T>::max;
+ }
+}
+
+
+template<class T>
+T meshRefinement::gAverage
+(
+ const polyMesh& mesh,
+ const PackedBoolList& isMasterElem,
+ const labelList& meshPoints,
+ const UList<T>& values
+)
+{
+ if (values.size() != meshPoints.size())
+ {
+ FatalErrorIn
+ (
+ "meshRefinement::gAverage\n"
+ "(\n"
+ " const polyMesh&,\n"
+ " const labelList&,\n"
+ " const PackedBoolList& isMasterElem,\n"
+ " const UList<T>& values\n"
+ ")\n"
+ ) << "Number of elements in list " << values.size()
+ << " does not correspond to number of elements in meshPoints "
+ << meshPoints.size()
+ << exit(FatalError);
+ }
+
+ T sum = pTraits<T>::zero;
+ label n = 0;
+
+ forAll(values, i)
+ {
+ if (isMasterElem[meshPoints[i]])
+ {
+ sum += values[i];
+ n++;
+ }
+ }
+
+ reduce(sum, sumOp<T>());
+ reduce(n, sumOp<label>());
+
+ if (n > 0)
+ {
+ return sum/n;
+ }
+ else
+ {
+ return pTraits<T>::max;
+ }
+}
+
+
// Compare two lists over all boundary faces
template<class T>
void meshRefinement::testSyncBoundaryFaceList
diff --git a/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C b/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C
index 2b90da54b72887fb97b2d01f74059d6a84d41215..b861bdb037221a652ed4575f991f2a8a60c0bf67 100644
--- a/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C
+++ b/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C
@@ -169,7 +169,7 @@ void Foam::trackedParticle::hitProcessorPatch
trackingData& td
)
{
- // Remove particle
+ // Move to different processor
td.switchProcessor = true;
}
diff --git a/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.H b/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.H
index de92709cf8c171943fe4a693678dc9fe3a7c046d..30acbc2fcc3bfe521f61deafca5806b93b766c47 100644
--- a/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.H
+++ b/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -163,12 +163,24 @@ public:
return end_;
}
+ //- transported label
+ label i() const
+ {
+ return i_;
+ }
+
//- transported label
label& i()
{
return i_;
}
+ //- transported label
+ label j() const
+ {
+ return j_;
+ }
+
//- transported label
label& j()
{
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 1d7770f050e8c322797f0166ef2185bc4f085b93..20b57d80b38a5ab7cb491961ba7dc6e05d64c2cd 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -135,145 +135,12 @@ Foam::sampledTriSurfaceMesh::nonCoupledboundaryTree() const
}
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
-(
- const word& name,
- const polyMesh& mesh,
- const word& surfaceName,
- const samplingSource sampleSource
-)
-:
- sampledSurface(name, mesh),
- surface_
- (
- IOobject
- (
- surfaceName,
- mesh.time().constant(), // instance
- "triSurface", // local
- mesh, // registry
- IOobject::MUST_READ,
- IOobject::NO_WRITE,
- false
- )
- ),
- sampleSource_(sampleSource),
- needsUpdate_(true),
- sampleElements_(0),
- samplePoints_(0)
-{}
-
-
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
-(
- const word& name,
- const polyMesh& mesh,
- const dictionary& dict
-)
-:
- sampledSurface(name, mesh, dict),
- surface_
- (
- IOobject
- (
- dict.lookup("surface"),
- mesh.time().constant(), // instance
- "triSurface", // local
- mesh, // registry
- IOobject::MUST_READ,
- IOobject::NO_WRITE,
- false
- )
- ),
- sampleSource_(samplingSourceNames_[dict.lookup("source")]),
- needsUpdate_(true),
- sampleElements_(0),
- samplePoints_(0)
-{}
-
-
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
-(
- const word& name,
- const polyMesh& mesh,
- const triSurface& surface,
- const word& sampleSourceName
-)
-:
- sampledSurface(name, mesh),
- surface_
- (
- IOobject
- (
- name,
- mesh.time().constant(), // instance
- "triSurface", // local
- mesh, // registry
- IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
- ),
- surface
- ),
- sampleSource_(samplingSourceNames_[sampleSourceName]),
- needsUpdate_(true),
- sampleElements_(0),
- samplePoints_(0)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-Foam::sampledTriSurfaceMesh::~sampledTriSurfaceMesh()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-bool Foam::sampledTriSurfaceMesh::needsUpdate() const
-{
- return needsUpdate_;
-}
-
-
-bool Foam::sampledTriSurfaceMesh::expire()
-{
- // already marked as expired
- if (needsUpdate_)
- {
- return false;
- }
-
- sampledSurface::clearGeom();
- MeshStorage::clear();
-
- boundaryTreePtr_.clear();
- sampleElements_.clear();
- samplePoints_.clear();
-
- needsUpdate_ = true;
- return true;
-}
-
-
-bool Foam::sampledTriSurfaceMesh::update()
+bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
{
- if (!needsUpdate_)
- {
- return false;
- }
-
-
// Find the cells the triangles of the surface are in.
// Does approximation by looking at the face centres only
const pointField& fc = surface_.faceCentres();
- // Mesh search engine, no triangulation of faces.
- meshSearch meshSearcher(mesh(), polyMesh::FACEPLANES);
-
-
List<nearInfo> nearest(fc.size());
// Global numbering for cells/faces - only used to uniquely identify local
@@ -611,6 +478,157 @@ bool Foam::sampledTriSurfaceMesh::update()
}
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+(
+ const word& name,
+ const polyMesh& mesh,
+ const word& surfaceName,
+ const samplingSource sampleSource
+)
+:
+ sampledSurface(name, mesh),
+ surface_
+ (
+ IOobject
+ (
+ surfaceName,
+ mesh.time().constant(), // instance
+ "triSurface", // local
+ mesh, // registry
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE,
+ false
+ )
+ ),
+ sampleSource_(sampleSource),
+ needsUpdate_(true),
+ sampleElements_(0),
+ samplePoints_(0)
+{}
+
+
+Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+(
+ const word& name,
+ const polyMesh& mesh,
+ const dictionary& dict
+)
+:
+ sampledSurface(name, mesh, dict),
+ surface_
+ (
+ IOobject
+ (
+ dict.lookup("surface"),
+ mesh.time().constant(), // instance
+ "triSurface", // local
+ mesh, // registry
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE,
+ false
+ )
+ ),
+ sampleSource_(samplingSourceNames_[dict.lookup("source")]),
+ needsUpdate_(true),
+ sampleElements_(0),
+ samplePoints_(0)
+{}
+
+
+Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+(
+ const word& name,
+ const polyMesh& mesh,
+ const triSurface& surface,
+ const word& sampleSourceName
+)
+:
+ sampledSurface(name, mesh),
+ surface_
+ (
+ IOobject
+ (
+ name,
+ mesh.time().constant(), // instance
+ "triSurface", // local
+ mesh, // registry
+ IOobject::NO_READ,
+ IOobject::NO_WRITE,
+ false
+ ),
+ surface
+ ),
+ sampleSource_(samplingSourceNames_[sampleSourceName]),
+ needsUpdate_(true),
+ sampleElements_(0),
+ samplePoints_(0)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::sampledTriSurfaceMesh::~sampledTriSurfaceMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+bool Foam::sampledTriSurfaceMesh::needsUpdate() const
+{
+ return needsUpdate_;
+}
+
+
+bool Foam::sampledTriSurfaceMesh::expire()
+{
+ // already marked as expired
+ if (needsUpdate_)
+ {
+ return false;
+ }
+
+ sampledSurface::clearGeom();
+ MeshStorage::clear();
+
+ boundaryTreePtr_.clear();
+ sampleElements_.clear();
+ samplePoints_.clear();
+
+ needsUpdate_ = true;
+ return true;
+}
+
+
+bool Foam::sampledTriSurfaceMesh::update()
+{
+ if (!needsUpdate_)
+ {
+ return false;
+ }
+
+ // Mesh search engine, no triangulation of faces.
+ meshSearch meshSearcher(mesh(), polyMesh::FACEPLANES);
+
+ return update(meshSearcher);
+}
+
+
+bool Foam::sampledTriSurfaceMesh::update(const treeBoundBox& bb)
+{
+ if (!needsUpdate_)
+ {
+ return false;
+ }
+
+ // Mesh search engine on subset, no triangulation of faces.
+ meshSearch meshSearcher(mesh(), bb, polyMesh::FACEPLANES);
+
+ return update(meshSearcher);
+}
+
+
Foam::tmp<Foam::scalarField> Foam::sampledTriSurfaceMesh::sample
(
const volScalarField& vField
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
index 1ea10b01e9b5f195ec7bce5e757a392440ad396b..7ef68570c43b0a67ac48c7b13005de6ca249bc0a 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
@@ -72,6 +72,7 @@ namespace Foam
{
class treeDataFace;
+class meshSearch;
/*---------------------------------------------------------------------------*\
Class sampledTriSurfaceMesh Declaration
@@ -136,6 +137,8 @@ private:
tmp<Field<Type> >
interpolateField(const interpolation<Type>&) const;
+ bool update(const meshSearch& meshSearcher);
+
public:
//- Runtime type information
@@ -189,6 +192,10 @@ public:
// Do nothing (and return false) if no update was needed
virtual bool update();
+ //- Update the surface using a bound box to limit the searching.
+ // For direct use, i.e. not through sample.
+ // Do nothing (and return false) if no update was needed
+ bool update(const treeBoundBox&);
//- Points of surface
virtual const pointField& points() const
diff --git a/src/thermophysicalModels/laminarFlameSpeed/Make/files b/src/thermophysicalModels/laminarFlameSpeed/Make/files
index 468c086eea8169d5892cf82ca95a81c22cf57292..63e0f3c1cae30a86857e36be92e252546e1b35d2 100644
--- a/src/thermophysicalModels/laminarFlameSpeed/Make/files
+++ b/src/thermophysicalModels/laminarFlameSpeed/Make/files
@@ -3,5 +3,6 @@ laminarFlameSpeed/laminarFlameSpeedNew.C
constant/constant.C
Gulders/Gulders.C
GuldersEGR/GuldersEGR.C
+RaviPetersen/RaviPetersen.C
LIB = $(FOAM_LIBBIN)/liblaminarFlameSpeedModels
diff --git a/src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.C b/src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.C
new file mode 100644
index 0000000000000000000000000000000000000000..320b779539eb35d9b87d854ae2da5f0c50f12acf
--- /dev/null
+++ b/src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.C
@@ -0,0 +1,365 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "RaviPetersen.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarFlameSpeedModels
+{
+ defineTypeNameAndDebug(RaviPetersen, 0);
+
+ addToRunTimeSelectionTable
+ (
+ laminarFlameSpeed,
+ RaviPetersen,
+ dictionary
+ );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::laminarFlameSpeedModels::RaviPetersen::RaviPetersen
+(
+ const dictionary& dict,
+ const psiuReactionThermo& ct
+)
+:
+ laminarFlameSpeed(dict, ct),
+ coeffsDict_(dict.subDict(typeName + "Coeffs").subDict(fuel_)),
+ pPoints_(coeffsDict_.lookup("pPoints")),
+ EqRPoints_(coeffsDict_.lookup("EqRPoints")),
+ alpha_(coeffsDict_.lookup("alpha")),
+ beta_(coeffsDict_.lookup("beta")),
+ TRef_(readScalar(coeffsDict_.lookup("TRef")))
+{
+ checkPointsMonotonicity("equivalenceRatio", EqRPoints_);
+ checkPointsMonotonicity("pressure", pPoints_);
+ checkCoefficientArrayShape("alpha", alpha_);
+ checkCoefficientArrayShape("beta", beta_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::laminarFlameSpeedModels::RaviPetersen::~RaviPetersen()
+{}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::laminarFlameSpeedModels::RaviPetersen::checkPointsMonotonicity
+(
+ const word& name,
+ const List<scalar>& x
+) const
+{
+ for (label i = 1; i < x.size(); i ++)
+ {
+ if (x[i] <= x[i-1])
+ {
+ FatalIOErrorIn
+ (
+ "laminarFlameSpeedModels::RaviPetersen::checkPointsMonotonicity"
+ "("
+ "const word&, "
+ "const List<scalar>&"
+ ") const",
+ coeffsDict_
+ ) << "Data points for the " << name
+ << " do not increase monotonically" << endl
+ << exit(FatalIOError);
+ }
+ }
+}
+
+
+void Foam::laminarFlameSpeedModels::RaviPetersen::checkCoefficientArrayShape
+(
+ const word& name,
+ const List<List<List<scalar> > >& x
+) const
+{
+ bool ok = true;
+
+ ok &= x.size() == EqRPoints_.size() - 1;
+
+ forAll(x, i)
+ {
+ ok &= x[i].size() == pPoints_.size();
+
+ forAll(x[i], j)
+ {
+ ok &= x[i][j].size() == x[i][0].size();
+ }
+ }
+
+ if (!ok)
+ {
+ FatalIOErrorIn
+ (
+ "laminarFlameSpeedModels::RaviPetersen::checkCoefficientArrayShape"
+ "("
+ "const word&, "
+ "const List<List<List<scalar> > >&"
+ ") const",
+ coeffsDict_
+ ) << "Inconsistent size of " << name << " coefficients array" << endl
+ << exit(FatalIOError);
+ }
+}
+
+
+inline bool Foam::laminarFlameSpeedModels::RaviPetersen::interval
+(
+ const List<scalar>& xPoints,
+ const scalar x,
+ label& xIndex,
+ scalar& xXi,
+ scalar& xLim
+) const
+{
+ if (x < xPoints.first())
+ {
+ xIndex = 0;
+ xXi = 0.0;
+ xLim = xPoints.first();
+ return false;
+ }
+
+ else if (x > xPoints.last())
+ {
+ xIndex = xPoints.size() - 2;
+ xXi = 1.0;
+ xLim = xPoints.last();
+ return false;
+ }
+
+ for (xIndex = 0; x > xPoints[xIndex+1]; xIndex ++)
+ {
+ // increment xIndex until xPoints[xIndex] < x < xPoints[xIndex+1]
+ }
+
+ xXi = (x - xPoints[xIndex])/(xPoints[xIndex+1] - xPoints[xIndex]);
+ xLim = x;
+
+ return true;
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::polynomial
+(
+ const List<scalar>& coeffs,
+ const scalar x
+) const
+{
+ scalar xPow = 1.0;
+ scalar y = 0.0;
+ forAll(coeffs, i)
+ {
+ y += coeffs[i]*xPow;
+ xPow *= x;
+ }
+ return y;
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::dPolynomial
+(
+ const List<scalar>& coeffs,
+ const scalar x
+) const
+{
+ scalar xPow = 1.0;
+ scalar y = 0.0;
+ for (label i = 1; i < coeffs.size(); i++)
+ {
+ y += i*coeffs[i]*xPow;
+ xPow *= x;
+ }
+ return y;
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::THatPowB
+(
+ const label EqRIndex,
+ const label pIndex,
+ const scalar EqR,
+ const scalar Tu
+) const
+{
+ return pow
+ (
+ Tu/TRef_,
+ polynomial(beta_[EqRIndex][pIndex],EqR)
+ );
+}
+
+
+inline Foam::scalar
+Foam::laminarFlameSpeedModels::RaviPetersen::correlationInRange
+(
+ const label EqRIndex,
+ const label pIndex,
+ const scalar EqR,
+ const scalar Tu
+) const
+{
+ // standard correlation
+ return
+ polynomial(alpha_[EqRIndex][pIndex],EqR)
+ *THatPowB(EqRIndex, pIndex, EqR, Tu);
+}
+
+
+inline Foam::scalar
+Foam::laminarFlameSpeedModels::RaviPetersen::correlationOutOfRange
+(
+ const label EqRIndex,
+ const label pIndex,
+ const scalar EqR,
+ const scalar EqRLim,
+ const scalar Tu
+) const
+{
+ scalar A = polynomial(alpha_[EqRIndex][pIndex], EqRLim);
+ scalar dA = dPolynomial(alpha_[EqRIndex][pIndex], EqRLim);
+ scalar dB = dPolynomial(beta_[EqRIndex][pIndex], EqRLim);
+ scalar TB = THatPowB(EqRIndex, pIndex, EqRLim, Tu);
+
+ // linear extrapolation from the bounds of the correlation
+ return max(TB*(A + (dA + A*log(Tu/TRef_)*dB)*(EqR - EqRLim)), 0.0);
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::speed
+(
+ const scalar EqR,
+ const scalar p,
+ const scalar Tu
+) const
+{
+ scalar Su = 0, s;
+
+ label EqRIndex, pIndex;
+ scalar EqRXi, pXi;
+ scalar EqRLim, pLim;
+ bool EqRInRange;
+
+ EqRInRange = interval(EqRPoints_, EqR, EqRIndex, EqRXi, EqRLim);
+
+ interval(pPoints_, p, pIndex, pXi, pLim);
+
+ for (label pI = 0; pI < 2; pI ++)
+ {
+ if (EqRInRange)
+ {
+ s = correlationInRange(EqRIndex, pIndex + pI, EqR, Tu);
+ }
+ else
+ {
+ s = correlationOutOfRange(EqRIndex, pIndex + pI, EqR, EqRLim, Tu);
+ }
+
+ Su += (1 - pXi)*s;
+ pXi = 1 - pXi;
+ }
+
+ return Su;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::laminarFlameSpeedModels::RaviPetersen::operator()() const
+{
+ const volScalarField& p = psiuReactionThermo_.p();
+ const volScalarField& Tu = psiuReactionThermo_.Tu();
+
+ volScalarField EqR
+ (
+ IOobject
+ (
+ "EqR",
+ p.time().timeName(),
+ p.db(),
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ p.mesh(),
+ dimensionedScalar("EqR", dimless, 0.0)
+ );
+
+ if (psiuReactionThermo_.composition().contains("ft"))
+ {
+ const volScalarField& ft = psiuReactionThermo_.composition().Y("ft");
+
+ EqR =
+ dimensionedScalar
+ (
+ psiuReactionThermo_.lookup("stoichiometricAirFuelMassRatio")
+ )*ft/max(1 - ft, SMALL);
+ }
+ else
+ {
+ EqR = equivalenceRatio_;
+ }
+
+ tmp<volScalarField> tSu0
+ (
+ new volScalarField
+ (
+ IOobject
+ (
+ "Su0",
+ p.time().timeName(),
+ p.db(),
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ p.mesh(),
+ dimensionedScalar("Su0", dimVelocity, 0.0)
+ )
+ );
+
+ volScalarField& Su0 = tSu0();
+
+ forAll(Su0, cellI)
+ {
+ Su0[cellI] = speed(EqR[cellI], p[cellI], Tu[cellI]);
+ }
+
+ return tSu0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.H b/src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.H
new file mode 100644
index 0000000000000000000000000000000000000000..964fee4032b00586a70ff80709aeb09f9bfe21dd
--- /dev/null
+++ b/src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.H
@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::laminarFlameSpeedModels::RaviPetersen
+
+Description
+ Laminar flame speed obtained from Ravi and Petersen's correlation.
+
+ The correlation for the laminar flame speed \f$Su\f$ is of the following
+ form:
+ \f[
+ Su = \left( \sum \alpha_i \phi^i \right)
+ \left( \frac{T}{T_{ref}} \right)^{\left( \sum \beta_j \phi^j \right)}
+ \f]
+
+ Where \f$\phi\f$ is the equivalence ratio, and \f$\alpha\f$ and \f$\beta\f$
+ are polynomial coefficients given for a number of pressure and equivalence
+ ratio points.
+
+SourceFiles
+ RaviPetersen.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef RaviPetersen_H
+#define RaviPetersen_H
+
+#include "laminarFlameSpeed.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarFlameSpeedModels
+{
+
+/*---------------------------------------------------------------------------*\
+ Class RaviPetersen Declaration
+\*---------------------------------------------------------------------------*/
+
+class RaviPetersen
+:
+ public laminarFlameSpeed
+{
+ // Private Data
+
+ dictionary coeffsDict_;
+
+ //- Correlation pressure values
+ List<scalar> pPoints_;
+
+ //- Correlation equivalence ratios
+ List<scalar> EqRPoints_;
+
+ //- Correlation alpha coefficients
+ List<List<List<scalar> > > alpha_;
+
+ //- Correlation beta coefficients
+ List<List<List<scalar> > > beta_;
+
+ //- Reference temperature
+ scalar TRef_;
+
+
+ // Private Member Functions
+
+ //- Check that input points are ordered
+ void checkPointsMonotonicity
+ (
+ const word& name,
+ const List<scalar>& x
+ ) const;
+
+ //- Check that the coefficient arrays are of the correct shape
+ void checkCoefficientArrayShape
+ (
+ const word& name,
+ const List<List<List<scalar> > >& x
+ ) const;
+
+ //- Find and interpolate a value in the data point arrays
+ inline bool interval
+ (
+ const List<scalar>& xPoints,
+ const scalar x,
+ label& xIndex,
+ scalar& xXi,
+ scalar& xLim
+ ) const;
+
+ //- Evaluate a polynomial
+ inline scalar polynomial
+ (
+ const List<scalar>& coeffs,
+ const scalar x
+ ) const;
+
+ //- Evaluate a polynomial differential
+ inline scalar dPolynomial
+ (
+ const List<scalar>& coeffs,
+ const scalar x
+ ) const;
+
+ //- Calculate normalised temperature to the power of the B polynomial
+ inline scalar THatPowB
+ (
+ const label EqRIndex,
+ const label pIndex,
+ const scalar EqR,
+ const scalar Tu
+ ) const;
+
+ //- Return the flame speed within the correlation range
+ inline scalar correlationInRange
+ (
+ const label EqRIndex,
+ const label pIndex,
+ const scalar EqR,
+ const scalar Tu
+ ) const;
+
+ //- Extrapolate the flame speed correlation outside its range
+ inline scalar correlationOutOfRange
+ (
+ const label EqRIndex,
+ const label pIndex,
+ const scalar EqR,
+ const scalar EqRLim,
+ const scalar Tu
+ ) const;
+
+ //- Return the laminar flame speed [m/s]
+ inline scalar speed
+ (
+ const scalar EqR,
+ const scalar p,
+ const scalar Tu
+ ) const;
+
+
+ //- Construct as copy (not implemented)
+ RaviPetersen(const RaviPetersen&);
+ void operator=(const RaviPetersen&);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("RaviPetersen");
+
+ // Constructors
+
+ //- Construct from dictionary and psiuReactionThermo
+ RaviPetersen
+ (
+ const dictionary&,
+ const psiuReactionThermo&
+ );
+
+
+ //- Destructor
+ virtual ~RaviPetersen();
+
+
+ // Member functions
+
+ //- Return the laminar flame speed [m/s]
+ tmp<volScalarField> operator()() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End laminarFlameSpeedModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/XiFoam/ras/Allclean b/tutorials/combustion/XiFoam/ras/Allclean
index 9768a9af4697ea684e8ac7c83ff95160a70b5c8f..b35a5b7730885770ad908e72e8d32cb9af8df8fe 100755
--- a/tutorials/combustion/XiFoam/ras/Allclean
+++ b/tutorials/combustion/XiFoam/ras/Allclean
@@ -5,7 +5,7 @@ cd ${0%/*} || exit 1 # run from this directory
. $WM_PROJECT_DIR/bin/tools/CleanFunctions
keepCases="moriyoshiHomogeneous"
-loseCases="moriyoshiHomogeneousPart2"
+loseCases="moriyoshiHomogeneousPart2 moriyoshiHomogeneousHydrogen"
for caseName in $keepCases
do
diff --git a/tutorials/combustion/XiFoam/ras/Allrun b/tutorials/combustion/XiFoam/ras/Allrun
index 352a424e97fc57620dfdbefe6006b5b697d900af..e533fb24c1489406634f5d01854ba9c6ab7e8d17 100755
--- a/tutorials/combustion/XiFoam/ras/Allrun
+++ b/tutorials/combustion/XiFoam/ras/Allrun
@@ -6,21 +6,25 @@ cd ${0%/*} || exit 1 # run from this directory
setControlDict()
{
- controlDict="system/controlDict"
-
sed \
- -e s/"\(deltaT[ \t]*\) 5e-06;"/"\1 1e-05;"/g \
- -e s/"\(endTime[ \t]*\) 0.005;"/"\1 0.015;"/g \
- -e s/"\(writeInterval[ \t]*\) 10;"/"\1 50;"/g \
- $controlDict > temp.$$
- mv temp.$$ $controlDict
+ -e "s/\(deltaT[ \t]*\) 5e-06;/\1 1e-05;/g" \
+ -e "s/\(endTime[ \t]*\) 0.005;/\1 0.015;/g" \
+ -e "s/\(writeInterval[ \t]*\) 10;/\1 50;/g" \
+ -i system/controlDict
}
+setCombustionProperties()
+{
+ sed \
+ -e "s/\(laminarFlameSpeedCorrelation[ \t]*\) Gulders;/\1 RaviPetersen;/g" \
+ -e "s/\(fuel[ \t]*\) Propane;/\1 HydrogenInAir;/g" \
+ -i constant/combustionProperties
+}
# Do moriyoshiHomogeneous
( cd moriyoshiHomogeneous &&