Commit 8e4035f2 authored by andy's avatar andy
Browse files

Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

parents 28030899 cb5796df
......@@ -44,6 +44,12 @@
fvc::interpolate(U) & mesh.Sf()
);
if (args.optionFound("initialiseUBCs"))
{
U.correctBoundaryConditions();
phi = fvc::interpolate(U) & mesh.Sf();
}
label pRefCell = 0;
scalar pRefValue = 0.0;
......
......@@ -38,6 +38,11 @@ Description
int main(int argc, char *argv[])
{
argList::addBoolOption("writep", "write the final pressure field");
argList::addBoolOption
(
"initialiseUBCs",
"initialise U boundary conditions"
);
#include "setRootCase.H"
#include "createTime.H"
......
......@@ -5,8 +5,8 @@
- fvm::Sp(fvc::div(phi), h)
- fvm::laplacian(turb.alphaEff(), h)
==
fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
- (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
fvc::div(phi/fvc::interpolate(rho), p, "div(U,p)")
- p*fvc::div(phi/fvc::interpolate(rho))
+ rad.Sh(thermo)
);
......
......@@ -844,16 +844,16 @@ int main(int argc, char *argv[])
List<faceList> patchFaceVerts;
labelList nrFaceCells(boundaryFaces.size(),0);
HashTable<label,label> faceToCell[2];
labelList own(boundaryFaces.size(), -1);
labelList nei(boundaryFaces.size(), -1);
HashTable<label, label> faceToCell[2];
{
HashTable<label, face, Hash<face> > faceToFaceID(boundaryFaces.size());
forAll(boundaryFaces, faceI)
{
SortableList<label> foo(boundaryFaces[faceI]);
face theFace(foo);
faceToFaceID.insert(theFace,faceI);
SortableList<label> sortedVerts(boundaryFaces[faceI]);
faceToFaceID.insert(face(sortedVerts), faceI);
}
forAll(cellVerts, cellI)
......@@ -861,31 +861,57 @@ int main(int argc, char *argv[])
faceList faces = cellVerts[cellI].faces();
forAll(faces, i)
{
SortableList<label> foo(faces[i]);
face theFace(foo);
if (faceToFaceID.found(theFace))
SortableList<label> sortedVerts(faces[i]);
HashTable<label, face, Hash<face> >::const_iterator fnd =
faceToFaceID.find(face(sortedVerts));
if (fnd != faceToFaceID.end())
{
label faceI = faceToFaceID[theFace];
if (nrFaceCells[faceI] < 2)
label faceI = fnd();
int stat = face::compare(faces[i], boundaryFaces[faceI]);
if (stat == 1)
{
// Same orientation. Cell is owner.
own[faceI] = cellI;
}
else if (stat == -1)
{
faceToCell[nrFaceCells[faceI]].insert(faceI,cellI);
// Opposite orientation. Cell is neighbour.
nei[faceI] = cellI;
}
nrFaceCells[faceI]++;
}
}
}
label nReverse = 0;
forAll(own, faceI)
{
if (own[faceI] == -1 && nei[faceI] != -1)
{
// Boundary face with incorrect orientation
boundaryFaces[faceI] = boundaryFaces[faceI].reverseFace();
Swap(own[faceI], nei[faceI]);
nReverse++;
}
}
if (nReverse > 0)
{
Info << "Found " << nReverse << " reversed boundary faces out of "
<< boundaryFaces.size() << endl;
}
label cnt = 0;
forAll(nrFaceCells, faceI)
forAll(own, faceI)
{
assert(nrFaceCells[faceI] == 1 || nrFaceCells[faceI] == 2);
if (nrFaceCells[faceI]>1)
if (own[faceI] != -1 && nei[faceI] != -1)
{
cnt++;
}
}
if (cnt>0)
if (cnt > 0)
{
Info << "Of " << boundaryFaces.size() << " so-called"
<< " boundary faces " << cnt << " belong to two cells "
......@@ -994,7 +1020,8 @@ int main(int argc, char *argv[])
if (boundaryFaceToIndex.found(faceIndices[i]))
{
label bFaceI = boundaryFaceToIndex[faceIndices[i]];
if (nrFaceCells[bFaceI] == 1)
if (own[bFaceI] != -1 && nei[bFaceI] == -1)
{
patchFaces[cnt] = boundaryFaces[bFaceI];
cnt++;
......
......@@ -1294,11 +1294,12 @@ Foam::labelList Foam::backgroundMeshDecomposition::processorPosition
"Foam::labelList"
"Foam::backgroundMeshDecomposition::processorPosition"
"("
"const List<point>& pts"
"const List<point>&"
") const"
)
<< "The position " << pts[pI]
<< " is not in any part of the background mesh. "
<< " is not in any part of the background mesh "
<< globalBackgroundBounds_ << endl
<< "A background mesh with a wider margin around "
<< "the geometry may help."
<< exit(FatalError);
......@@ -1311,11 +1312,11 @@ Foam::labelList Foam::backgroundMeshDecomposition::processorPosition
"Foam::labelList"
"Foam::backgroundMeshDecomposition::processorPosition"
"("
"const List<point>& pts"
"const List<point>&"
") const"
)
<< "The position " << pts[pI]
<< " was not found in the background mesh, finding nearest."
) << "The position " << pts[pI]
<< " was not found in the background mesh "
<< globalBackgroundBounds_ << ", finding nearest."
<< endl;
}
......
......@@ -30,6 +30,16 @@ License
#include "backgroundMeshDecomposition.H"
#include "meshSearch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(conformalVoronoiMesh, 0);
}
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
......
......@@ -722,6 +722,13 @@ private:
const Delaunay::Finite_facets_iterator& fit
) const;
//- Determines if the edge constructed from the face is on
// a processor patch
inline bool isParallelDualEdge
(
const Delaunay::Finite_facets_iterator& fit
) const;
//- Merge vertices that are very close together
void mergeCloseDualVertices
(
......@@ -771,7 +778,7 @@ private:
) const;
//- Collapse a face to an edge, updating the point and point
// map. Returns the collapse mode that was applied.
// map. Returns the collapse mode that was applied.
faceCollapseMode collapseFace
(
const face& f,
......@@ -783,7 +790,7 @@ private:
label maxFC
) const;
// Identify the index of the longest edge on the face
//- Identify the index of the longest edge on the face
label longestEdge(const face& f, const pointField& pts) const;
//- Identify the face labels of the deferred collapse faces
......@@ -886,6 +893,10 @@ private:
public:
//- Runtime type information
ClassName("conformalVoronoiMesh");
// Constructors
//- Construct from Time and cvMeshDict
......
......@@ -543,8 +543,21 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
{
label nPtsMerged = 0;
label nIdentical = 0;
label nProcEdge = 0;
// Relative distance for points to be merged
scalar closenessTolerance = cvMeshControls().mergeClosenessCoeff();
// Absolute distance for points to be considered coincident. Bit adhoc
// but points were seen with distSqr ~ 1E-30 which is SMALL^2. Add a few
// digits to account for truncation errors.
scalar coincidentDistanceSqr = sqr
(
SMALL*1E2*geometryToConformTo_.globalBounds().mag()
);
for
(
Delaunay::Finite_facets_iterator fit = finite_facets_begin();
......@@ -567,14 +580,14 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
continue;
}
if (!c1->farCell() && !c2->farCell() && (c1I != c2I))
if ((c1I != c2I) && !c1->farCell() && !c2->farCell())
{
if
(
magSqr(pts[c1I] - pts[c2I])
< sqr(averageAnyCellSize(fit)*closenessTolerance)
)
scalar distSqr = magSqr(pts[c1I] - pts[c2I]);
if (pts[c1I] == pts[c2I] || distSqr < coincidentDistanceSqr)
{
nIdentical++;
if (boundaryPts[c2I] == true)
{
// If c2I is a boundary point, then it is kept.
......@@ -583,18 +596,59 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
dualPtIndexMap.insert(c1I, c2I);
dualPtIndexMap.insert(c2I, c2I);
nPtsMerged++;
}
else
{
dualPtIndexMap.insert(c1I, c1I);
dualPtIndexMap.insert(c2I, c1I);
nPtsMerged++;
}
nPtsMerged++;
}
else if (distSqr < sqr(averageAnyCellSize(fit)*closenessTolerance))
{
if (c1->parallelDualVertex() || c2->parallelDualVertex())
//if (isParallelDualEdge(fit))
{
// Skip if face uses any edge that becomes a processor
// dual face.
// Note: the real test should be whether the Delaunay edge
// will form a processor patch.
nProcEdge++;
}
else if (boundaryPts[c2I] == true)
{
// If c2I is a boundary point, then it is kept.
// If both are boundary points then c2I is chosen
// arbitrarily to be kept.
dualPtIndexMap.insert(c1I, c2I);
dualPtIndexMap.insert(c2I, c2I);
nPtsMerged++;
}
else
{
dualPtIndexMap.insert(c1I, c1I);
dualPtIndexMap.insert(c2I, c1I);
nPtsMerged++;
}
}
}
}
if (debug)
{
Info<< "mergeCloseDualVertices:"
<< " coincident distance:" << coincidentDistanceSqr
<< " closenessTolerance:" << closenessTolerance << endl
<< " identical points : "
<< returnReduce(nIdentical, sumOp<label>()) << endl
<< " processor edges : "
<< returnReduce(nProcEdge, sumOp<label>()) << endl
<< endl;
}
return nPtsMerged;
}
......
......@@ -384,6 +384,23 @@ inline Foam::List<Foam::label> Foam::conformalVoronoiMesh::processorsAttached
}
inline bool Foam::conformalVoronoiMesh::isParallelDualEdge
(
const Delaunay::Finite_facets_iterator& fit
) const
{
const Cell_handle c1(fit->first);
const int oppositeVertex = fit->second;
return
(
c1->vertex(vertex_triple_index(oppositeVertex, 0))->referred()
|| c1->vertex(vertex_triple_index(oppositeVertex, 1))->referred()
|| c1->vertex(vertex_triple_index(oppositeVertex, 2))->referred()
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
#ifdef CGAL_INEXACT
......
......@@ -208,7 +208,6 @@ public:
);
}
//- Does the Dual vertex form part of a processor patch
inline Foam::label dualVertexMasterProc() const
{
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -55,7 +55,6 @@ Description
this if you don't mind having disconnected domains in a single region.
This option requires all cells to be in one (and one only) cellZone.
- Should work in parallel.
cellZones can differ on either side of processor boundaries in which case
the faces get moved from processor patch to directMapped patch. Not
......@@ -66,6 +65,13 @@ Description
region that covers more than 50% of the zone. It has to be a subset
so cannot have any cells in any other zone.
- If explicitly a single region has been selected (-largestOnly or
-insidePoint) its region name will be either
- name of a cellZone it matches to or
- "largestOnly" respectively "insidePoint" or
- polyMesh::defaultRegion if additionally -overwrite
(so it will overwrite the input mesh!)
- writes maps like decomposePar back to original mesh:
- pointRegionAddressing : for every point in this region the point in
the original mesh
......@@ -1866,6 +1872,29 @@ int main(int argc, char *argv[])
regionNames,
zoneToRegion
);
// Override any default region names if single region selected
if (largestOnly || insidePoint)
{
forAll(regionToZone, regionI)
{
if (regionToZone[regionI] == -1)
{
if (overwrite)
{
regionNames[regionI] = polyMesh::defaultRegion;
}
else if (insidePoint)
{
regionNames[regionI] = "insidePoint";
}
else if (largestOnly)
{
regionNames[regionI] = "largestOnly";
}
}
}
}
}
Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
......@@ -2173,7 +2202,8 @@ int main(int argc, char *argv[])
Info<< nl
<< "Subsetting region " << regionI
<< " of size " << regionSizes[regionI] << endl;
<< " of size " << regionSizes[regionI]
<< " as named region " << regionNames[regionI] << endl;
createAndWriteRegion
(
......
......@@ -541,29 +541,51 @@ int main(int argc, char *argv[])
IOobjectList objects(mesh, runTime.timeName());
HashSet<word> selectedFields;
args.optionReadIfPresent("fields", selectedFields);
bool specifiedFields = args.optionReadIfPresent
(
"fields",
selectedFields
);
// Construct the vol fields (on the original mesh if subsetted)
PtrList<volScalarField> vsf;
readFields(vMesh, vMesh.baseMesh(), objects, selectedFields, vsf);
print(" volScalarFields :", Info, vsf);
PtrList<volVectorField> vvf;
readFields(vMesh, vMesh.baseMesh(), objects, selectedFields, vvf);
print(" volVectorFields :", Info, vvf);
PtrList<volSphericalTensorField> vSpheretf;
readFields(vMesh, vMesh.baseMesh(), objects, selectedFields, vSpheretf);
print(" volSphericalTensorFields :", Info, vSpheretf);
PtrList<volSymmTensorField> vSymmtf;
readFields(vMesh, vMesh.baseMesh(), objects, selectedFields, vSymmtf);
print(" volSymmTensorFields :", Info, vSymmtf);
PtrList<volTensorField> vtf;
readFields(vMesh, vMesh.baseMesh(), objects, selectedFields, vtf);
print(" volTensorFields :", Info, vtf);
if (!specifiedFields || selectedFields.size())
{
readFields(vMesh, vMesh.baseMesh(), objects, selectedFields, vsf);
print(" volScalarFields :", Info, vsf);
readFields(vMesh, vMesh.baseMesh(), objects, selectedFields, vvf);
print(" volVectorFields :", Info, vvf);
readFields
(
vMesh,
vMesh.baseMesh(),
objects,
selectedFields,
vSpheretf
);
print(" volSphericalTensorFields :", Info, vSpheretf);
readFields
(
vMesh,
vMesh.baseMesh(),
objects,
selectedFields,
vSymmtf
);
print(" volSymmTensorFields :", Info, vSymmtf);
readFields(vMesh, vMesh.baseMesh(), objects, selectedFields, vtf);
print(" volTensorFields :", Info, vtf);
}
label nVolFields =
vsf.size()
......@@ -589,7 +611,7 @@ int main(int argc, char *argv[])
PtrList<pointSymmTensorField> pSymmtf;
PtrList<pointTensorField> ptf;
if (!noPointValues)
if (!noPointValues && !(specifiedFields && selectedFields.empty()))
{
readFields
(
......
......@@ -263,6 +263,26 @@ checkLineLengthNonDirective()
}
#
# check that OpenCFD copyright is currents
#
checkCopyright()
{
year=$(date +%Y)
echo "$hookName: check copyright ..." 1>&2
for f in $fileList
do
sYear=`grep "Copyright.*OpenCFD" $f | sed 's/[^0-9]//g' | cut -c 5-9`
if [ "$year" != "" ] && [ "$year" != "$sYear" ]; then
echo "Updated copyright for: $f"
sed -i "s/$sYear OpenCFD/$year OpenCFD/g" $f
fi
done
}
#------------------------------------------------------------------------------
# Main code : do all checks
#
......@@ -276,6 +296,7 @@ checkIllegalCode
# ensure code conforms to 80 columns max
checkLineLengthNonDirective
checkCopyright
exit 0
#------------------------------------------------------------------------------
......@@ -66,8 +66,13 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
// Master reads headerclassname from file. Make sure this gets
// transfered as well as contents.
Pstream::scatter(comms, const_cast<word&>(headerClassName()));
Pstream::scatter(comms, note());
Pstream::scatter
(
comms,
const_cast<word&>(headerClassName()),
Pstream::msgType()
);
Pstream::scatter(comms, note(), Pstream::msgType());
// Get my communication order
const Pstream::commsStruct& myComm = comms[Pstream::myProcNo()];
......@@ -88,6 +93,7 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
Pstream::scheduled,
myComm.above(),
0,
Pstream::msgType(),
IOstream::ASCII
);
IOdictionary::readData(fromAbove);
......