-
Andrew Heather authoredAndrew Heather authored
surfaceFieldValue.C 19.83 KiB
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 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 "surfaceFieldValue.H"
#include "fvMesh.H"
#include "emptyPolyPatch.H"
#include "coupledPolyPatch.H"
#include "sampledSurface.H"
#include "mergePoints.H"
#include "indirectPrimitivePatch.H"
#include "PatchTools.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace functionObjects
{
namespace fieldValues
{
defineTypeNameAndDebug(surfaceFieldValue, 0);
addToRunTimeSelectionTable(fieldValue, surfaceFieldValue, dictionary);
addToRunTimeSelectionTable(functionObject, surfaceFieldValue, dictionary);
}
}
}
template<>
const char* Foam::NamedEnum
<
Foam::functionObjects::fieldValues::surfaceFieldValue::regionTypes,
3
>::names[] =
{
"faceZone",
"patch",
"sampledSurface"
};
template<>
const char* Foam::NamedEnum
<
Foam::functionObjects::fieldValues::surfaceFieldValue::operationType,
16
>::names[] =
{
"none",
"sum",
"sumMag",
"sumDirection",
"sumDirectionBalance",
"average",
"weightedAverage",
"areaAverage",
"weightedAreaAverage",
"areaIntegrate",
"weightedAreaIntegrate",
"min",
"max",
"CoV",
"areaNormalAverage",
"areaNormalIntegrate"
};
template<>
const char* Foam::NamedEnum
<
Foam::functionObjects::fieldValues::surfaceFieldValue::postOperationType,
2
>::names[] =
{
"none",
"sqrt"
};
const Foam::NamedEnum
<
Foam::functionObjects::fieldValues::surfaceFieldValue::regionTypes,
3
> Foam::functionObjects::fieldValues::surfaceFieldValue::regionTypeNames_;
const Foam::NamedEnum
<
Foam::functionObjects::fieldValues::surfaceFieldValue::operationType,
16
> Foam::functionObjects::fieldValues::surfaceFieldValue::operationTypeNames_;
const Foam::NamedEnum
<
Foam::functionObjects::fieldValues::surfaceFieldValue::postOperationType,
2
>
Foam::functionObjects::fieldValues::surfaceFieldValue::postOperationTypeNames_;
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::functionObjects::fieldValues::surfaceFieldValue::setFaceZoneFaces()
{
label zoneId = mesh_.faceZones().findZoneID(regionName_);
if (zoneId < 0)
{
FatalErrorInFunction
<< type() << " " << name() << ": "
<< regionTypeNames_[regionType_] << "(" << regionName_ << "):" << nl
<< " Unknown face zone name: " << regionName_
<< ". Valid face zones are: " << mesh_.faceZones().names()
<< nl << exit(FatalError);
}
const faceZone& fZone = mesh_.faceZones()[zoneId];
DynamicList<label> faceIds(fZone.size());
DynamicList<label> facePatchIds(fZone.size());
DynamicList<label> faceSigns(fZone.size());
forAll(fZone, i)
{
label facei = fZone[i];
label faceId = -1;
label facePatchId = -1;
if (mesh_.isInternalFace(facei))
{
faceId = facei;
facePatchId = -1;
}
else
{
facePatchId = mesh_.boundaryMesh().whichPatch(facei);
const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
if (isA<coupledPolyPatch>(pp))
{
if (refCast<const coupledPolyPatch>(pp).owner())
{
faceId = pp.whichFace(facei);
}
else
{
faceId = -1;
}
}
else if (!isA<emptyPolyPatch>(pp))
{
faceId = facei - pp.start();
}
else
{
faceId = -1;
facePatchId = -1;
}
}
if (faceId >= 0)
{
if (fZone.flipMap()[i])
{
faceSigns.append(-1);
}
else
{
faceSigns.append(1);
}
faceIds.append(faceId);
facePatchIds.append(facePatchId);
}
}
faceId_.transfer(faceIds);
facePatchId_.transfer(facePatchIds);
faceSign_.transfer(faceSigns);
nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
if (debug)
{
Pout<< "Original face zone size = " << fZone.size()
<< ", new size = " << faceId_.size() << endl;
}
}
void Foam::functionObjects::fieldValues::surfaceFieldValue::setPatchFaces()
{
const label patchid = mesh_.boundaryMesh().findPatchID(regionName_);
if (patchid < 0)
{
FatalErrorInFunction
<< type() << " " << name() << ": "
<< regionTypeNames_[regionType_] << "(" << regionName_ << "):" << nl
<< " Unknown patch name: " << regionName_
<< ". Valid patch names are: "
<< mesh_.boundaryMesh().names() << nl
<< exit(FatalError);
}
const polyPatch& pp = mesh_.boundaryMesh()[patchid];
label nFaces = pp.size();
if (isA<emptyPolyPatch>(pp))
{
nFaces = 0;
}
faceId_.setSize(nFaces);
facePatchId_.setSize(nFaces);
faceSign_.setSize(nFaces);
nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
forAll(faceId_, facei)
{
faceId_[facei] = facei;
facePatchId_[facei] = patchid;
faceSign_[facei] = 1;
}
}
void Foam::functionObjects::fieldValues::surfaceFieldValue::sampledSurfaceFaces
(
const dictionary& dict
)
{
surfacePtr_ = sampledSurface::New
(
name(),
mesh_,
dict.subDict("sampledSurfaceDict")
);
surfacePtr_().update();
nFaces_ = returnReduce(surfacePtr_().faces().size(), sumOp<label>());
}
void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry
(
faceList& faces,
pointField& points
) const
{
List<faceList> allFaces(Pstream::nProcs());
List<pointField> allPoints(Pstream::nProcs());
labelList globalFacesIs(faceId_);
forAll(globalFacesIs, i)
{
if (facePatchId_[i] != -1)
{
label patchi = facePatchId_[i];
globalFacesIs[i] += mesh_.boundaryMesh()[patchi].start();
}
}
// Add local faces and points to the all* lists
indirectPrimitivePatch pp
(
IndirectList<face>(mesh_.faces(), globalFacesIs),
mesh_.points()
);
allFaces[Pstream::myProcNo()] = pp.localFaces();
allPoints[Pstream::myProcNo()] = pp.localPoints();
Pstream::gatherList(allFaces);
Pstream::gatherList(allPoints);
// Renumber and flatten
label nFaces = 0;
label nPoints = 0;
forAll(allFaces, proci)
{
nFaces += allFaces[proci].size();
nPoints += allPoints[proci].size();
}
faces.setSize(nFaces);
points.setSize(nPoints);
nFaces = 0;
nPoints = 0;
// My own data first
{
const faceList& fcs = allFaces[Pstream::myProcNo()];
forAll(fcs, i)
{
const face& f = fcs[i];
face& newF = faces[nFaces++];
newF.setSize(f.size());
forAll(f, fp)
{
newF[fp] = f[fp] + nPoints;
}
}
const pointField& pts = allPoints[Pstream::myProcNo()];
forAll(pts, i)
{
points[nPoints++] = pts[i];
}
}
// Other proc data follows
forAll(allFaces, proci)
{
if (proci != Pstream::myProcNo())
{
const faceList& fcs = allFaces[proci];
forAll(fcs, i)
{
const face& f = fcs[i];
face& newF = faces[nFaces++];
newF.setSize(f.size());
forAll(f, fp)
{
newF[fp] = f[fp] + nPoints;
}
}
const pointField& pts = allPoints[proci];
forAll(pts, i)
{
points[nPoints++] = pts[i];
}
}
}
// Merge
labelList oldToNew;
pointField newPoints;
bool hasMerged = mergePoints
(
points,
SMALL,
false,
oldToNew,
newPoints
);
if (hasMerged)
{
if (debug)
{
Pout<< "Merged from " << points.size()
<< " down to " << newPoints.size() << " points" << endl;
}
points.transfer(newPoints);
forAll(faces, i)
{
inplaceRenumber(oldToNew, faces[i]);
}
}
}
void Foam::functionObjects::fieldValues::surfaceFieldValue::
combineSurfaceGeometry
(
faceList& faces,
pointField& points
) const
{
if (surfacePtr_.valid())
{
const sampledSurface& s = surfacePtr_();
if (Pstream::parRun())
{
// Dimension as fraction of mesh bounding box
scalar mergeDim = 1e-10*mesh_.bounds().mag();
labelList pointsMap;
PatchTools::gatherAndMerge
(
mergeDim,
primitivePatch
(
SubList<face>(s.faces(), s.faces().size()),
s.points()
),
points,
faces,
pointsMap
);
}
else
{
faces = s.faces();
points = s.points();
}
}
}
Foam::scalar
Foam::functionObjects::fieldValues::surfaceFieldValue::totalArea() const
{
scalar totalArea;
if (surfacePtr_.valid())
{
totalArea = gSum(surfacePtr_().magSf());
}
else
{
totalArea = gSum(filterField(mesh_.magSf(), false));
}
return totalArea;
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::functionObjects::fieldValues::surfaceFieldValue::initialise
(
const dictionary& dict
)
{
dict.lookup("name") >> regionName_;
switch (regionType_)
{
case stFaceZone:
{
setFaceZoneFaces();
break;
}
case stPatch:
{
setPatchFaces();
break;
}
case stSampledSurface:
{
sampledSurfaceFaces(dict);
break;
}
default:
{
FatalErrorInFunction
<< type() << " " << name() << ": "
<< regionTypeNames_[regionType_] << "(" << regionName_ << "):"
<< nl << " Unknown region type. Valid region types are:"
<< regionTypeNames_.sortedToc() << nl << exit(FatalError);
}
}
if (nFaces_ == 0)
{
FatalErrorInFunction
<< type() << " " << name() << ": "
<< regionTypeNames_[regionType_] << "(" << regionName_ << "):" << nl
<< " Region has no faces" << exit(FatalError);
}
if (surfacePtr_.valid())
{
surfacePtr_().update();
}
totalArea_ = totalArea();
Info<< type() << " " << name() << ":" << nl
<< " total faces = " << nFaces_ << nl
<< " total area = " << totalArea_ << nl;
if (dict.readIfPresent("weightField", weightFieldName_))
{
Info<< " weight field = " << weightFieldName_ << nl;
if (regionType_ == stSampledSurface)
{
FatalIOErrorInFunction(dict)
<< "Cannot use weightField for a sampledSurface"
<< exit(FatalIOError);
}
}
if (dict.found("orientedWeightField"))
{
if (weightFieldName_ == "none")
{
dict.lookup("orientedWeightField") >> weightFieldName_;
Info<< " weight field = " << weightFieldName_ << nl;
orientWeightField_ = true;
}
else
{
FatalIOErrorInFunction(dict)
<< "Either weightField or orientedWeightField can be supplied, "
<< "but not both"
<< exit(FatalIOError);
}
}
List<word> orientedFields;
if (dict.readIfPresent("orientedFields", orientedFields))
{
orientedFieldsStart_ = fields_.size();
fields_.append(orientedFields);
}
Info<< nl << endl;
if (writeFields_)
{
const word surfaceFormat(dict.lookup("surfaceFormat"));
surfaceWriterPtr_.reset
(
surfaceWriter::New
(
surfaceFormat,
dict.subOrEmptyDict("formatOptions").
subOrEmptyDict(surfaceFormat)
).ptr()
);
}
}
void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader
(
Ostream& os
) const
{
if (operation_ != opNone)
{
writeCommented(os, "Region type : ");
os << regionTypeNames_[regionType_] << " " << regionName_ << endl;
writeHeaderValue(os, "Faces", nFaces_);
writeHeaderValue(os, "Area", totalArea_);
writeHeaderValue(os, "Scale factor", scaleFactor_);
if (weightFieldName_ != "none")
{
writeHeaderValue(os, "Weight field", weightFieldName_);
}
writeCommented(os, "Time");
if (writeArea_)
{
os << tab << "Area";
}
forAll(fields_, i)
{
os << tab << operationTypeNames_[operation_]
<< "(" << fields_[i] << ")";
}
os << endl;
}
}
template<>
Foam::scalar Foam::functionObjects::fieldValues::surfaceFieldValue::
processValues
(
const Field<scalar>& values,
const vectorField& Sf,
const scalarField& weightField
) const
{
switch (operation_)
{
case opSumDirection:
{
vector n(dict_.lookup("direction"));
return gSum(pos(values*(Sf & n))*mag(values));
}
case opSumDirectionBalance:
{
vector n(dict_.lookup("direction"));
const scalarField nv(values*(Sf & n));
return gSum(pos(nv)*mag(values) - neg(nv)*mag(values));
}
default:
{
// Fall through to other operations
return processSameTypeValues(values, Sf, weightField);
}
}
}
template<>
Foam::vector Foam::functionObjects::fieldValues::surfaceFieldValue::
processValues
(
const Field<vector>& values,
const vectorField& Sf,
const scalarField& weightField
) const
{
switch (operation_)
{
case opSumDirection:
{
vector n(dict_.lookup("direction"));
n /= mag(n) + ROOTVSMALL;
const scalarField nv(n & values);
return gSum(pos(nv)*n*(nv));
}
case opSumDirectionBalance:
{
vector n(dict_.lookup("direction"));
n /= mag(n) + ROOTVSMALL;
const scalarField nv(n & values);
return gSum(pos(nv)*n*(nv));
}
case opAreaNormalAverage:
{
scalar result = gSum(values & Sf)/gSum(mag(Sf));
return vector(result, 0.0, 0.0);
}
case opAreaNormalIntegrate:
{
scalar result = gSum(values & Sf);
return vector(result, 0.0, 0.0);
}
default:
{
// Fall through to other operations
return processSameTypeValues(values, Sf, weightField);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
(
const word& name,
const Time& runTime,
const dictionary& dict
)
:
fieldValue(name, runTime, dict, typeName),
surfaceWriterPtr_(nullptr),
regionType_(regionTypeNames_.read(dict.lookup("regionType"))),
operation_(operationTypeNames_.read(dict.lookup("operation"))),
postOperation_
(
postOperationTypeNames_
[dict.lookupOrDefault<word>("postOperation", "none")]
),
weightFieldName_("none"),
orientWeightField_(false),
orientedFieldsStart_(labelMax),
writeArea_(dict.lookupOrDefault("writeArea", false)),
nFaces_(0),
faceId_(),
facePatchId_(),
faceSign_()
{
read(dict);
writeFileHeader(file());
}
Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
(
const word& name,
const objectRegistry& obr,
const dictionary& dict
)
:
fieldValue(name, obr, dict, typeName),
surfaceWriterPtr_(nullptr),
regionType_(regionTypeNames_.read(dict.lookup("regionType"))),
operation_(operationTypeNames_.read(dict.lookup("operation"))),
postOperation_
(
postOperationTypeNames_
[dict.lookupOrDefault<word>("postOperation", "none")]
),
weightFieldName_("none"),
orientWeightField_(false),
orientedFieldsStart_(labelMax),
writeArea_(dict.lookupOrDefault("writeArea", false)),
nFaces_(0),
faceId_(),
facePatchId_(),
faceSign_()
{
read(dict);
writeFileHeader(file());
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::functionObjects::fieldValues::surfaceFieldValue::~surfaceFieldValue()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
(
const dictionary& dict
)
{
fieldValue::read(dict);
initialise(dict);
return true;
}
bool Foam::functionObjects::fieldValues::surfaceFieldValue::write()
{
if (operation_ != opNone)
{
fieldValue::write();
}
if (surfacePtr_.valid())
{
surfacePtr_().update();
}
if (operation_ != opNone && Pstream::master())
{
writeTime(file());
}
if (writeArea_)
{
totalArea_ = totalArea();
if (operation_ != opNone && Pstream::master())
{
file() << tab << totalArea_;
}
Log << " total area = " << totalArea_ << endl;
}
// Construct weight field. Note: zero size means weight = 1
scalarField weightField;
if (weightFieldName_ != "none")
{
weightField =
getFieldValues<scalar>
(
weightFieldName_,
true,
orientWeightField_
);
}
// Process the fields
forAll(fields_, i)
{
const word& fieldName = fields_[i];
bool ok = false;
bool orient = i >= orientedFieldsStart_;
ok = ok || writeValues<scalar>(fieldName, weightField, orient);
ok = ok || writeValues<vector>(fieldName, weightField, orient);
ok = ok
|| writeValues<sphericalTensor>(fieldName, weightField, orient);
ok = ok || writeValues<symmTensor>(fieldName, weightField, orient);
ok = ok || writeValues<tensor>(fieldName, weightField, orient);
if (!ok)
{
WarningInFunction
<< "Requested field " << fieldName
<< " not found in database and not processed"
<< endl;
}
}
if (operation_ != opNone && Pstream::master())
{
file() << endl;
}
Log << endl;
return true;
}
// ************************************************************************* //