Merge branch 'updated-surface-handling' into 'develop'

Updated surface handling

See merge request Development/openfoam!394
This commit is contained in:
Andrew Heather
2020-11-30 17:13:39 +00:00
36 changed files with 3055 additions and 440 deletions

View File

@ -47,12 +47,19 @@ Usage
where the entries mean:
\table
Property | Description | Type | Req'd | Dflt
writePrecision | Number of decimal points | label | no | \<system dflt\>
writeToFile | Flag to produce text file output | bool | no | true
useUserTime | Flag to use user time, e.g. degrees | bool | no | true
Property | Description | Type | Reqd | Dflt
writePrecision | Number of decimal points | int | no | \<system dflt\>
writeToFile | Produce text file output? | bool | no | true
useUserTime | Use user time (e.g. degrees)? | bool | no | true
updateHeader | Update header on mesh changes? | bool | no | true
\endtable
Note
The file header is normally updated whenver the mesh points or
topology changes. In some cases, the function object is actually
unaffected by these changes.
Use the \c updateHeader flag to override the default behaviour.
See also
- Foam::functionObject
- Foam::functionObjects::logFiles

View File

@ -30,8 +30,9 @@ License
#include "ListOps.H"
#include "stringOps.H"
#include "UIListStream.H"
#undef Foam_readAbaqusSurface
#include "cellModel.H"
#include <algorithm>
#include <cctype>
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -319,11 +320,10 @@ Foam::fileFormats::ABAQUSCore::readHelper::readPoints
const label initialCount = points_.size();
char sep; // Comma separator (dummy)
string line;
label id;
point p;
string line;
// Read nodes (points) until next "*Section"
while (is.peek() != '*' && is.peek() != EOF)
{
@ -368,10 +368,10 @@ Foam::fileFormats::ABAQUSCore::readHelper::readElements
const label initialCount = elemTypes_.size();
char sep; // Comma separator (dummy)
label id;
labelList elemNodes(nNodes, Zero);
string line;
label id;
labelList elemNodes(nNodes, Zero);
// Read element connectivity until next "*Section"
@ -403,6 +403,130 @@ Foam::fileFormats::ABAQUSCore::readHelper::readElements
}
Foam::label
Foam::fileFormats::ABAQUSCore::readHelper::readSurfaceElements
(
ISstream& is,
const label setId
)
{
// Info<< "*Surface" << nl;
// Models for supported solids (need to face mapping)
const cellModel& tet = cellModel::ref(cellModel::TET);
const cellModel& prism = cellModel::ref(cellModel::PRISM);
const cellModel& hex = cellModel::ref(cellModel::HEX);
// Face mapping from Abaqus cellModel to OpenFOAM cellModel
const auto& abqToFoamFaceMap = abaqusToFoamFaceAddr();
const label initialCount = elemTypes_.size();
char sep; // Comma separator (dummy)
string line;
label id;
// Read until next "*Section"
// Parse for elemId, sideId.
// Eg, "1235, S1"
while (is.peek() != '*' && is.peek() != EOF)
{
is >> id >> sep;
is.getLine(line);
const word sideName(word::validate(stringOps::upper(line)));
if
(
sideName.size() != 2
|| sideName[0] != 'S'
|| !std::isdigit(sideName[1])
)
{
Info<< "Abaqus reader: unsupported surface element side "
<< id << ", " << sideName << nl;
continue;
}
const label index = elemIds_.find(id);
if (id <= 0 || index < 0)
{
Info<< "Abaqus reader: unsupported surface element "
<< id << nl;
continue;
}
const auto faceIdIter = abqToFoamFaceMap.cfind(elemTypes_[index]);
if (!faceIdIter.found())
{
Info<< "Abaqus reader: reject non-solid shape: " << nl;
}
// The abaqus element side number (1-based)
const label sideNum = (sideName[1] - '0');
const label foamFaceNum = (*faceIdIter)[sideNum - 1];
const labelList& connect = connectivity_[index];
// Nodes for the derived shell element
labelList elemNodes;
switch (elemTypes_[index])
{
case shapeType::abaqusTet:
{
elemNodes = labelList(connect, tet.modelFaces()[foamFaceNum]);
break;
}
case shapeType::abaqusPrism:
{
elemNodes = labelList(connect, prism.modelFaces()[foamFaceNum]);
break;
}
case shapeType::abaqusHex:
{
elemNodes = labelList(connect, hex.modelFaces()[foamFaceNum]);
break;
}
default:
break;
}
enum shapeType shape = shapeType::abaqusUnknownShape;
if (elemNodes.size() == 3)
{
shape = shapeType::abaqusTria;
}
else if (elemNodes.size() == 4)
{
shape = shapeType::abaqusQuad;
}
else
{
// Cannot happen
FatalErrorInFunction
<< "Could not map face side for "
<< id << ", " << sideName << nl
<< exit(FatalError);
}
// Synthesize face Id from solid element Id and side Id
const label newElemId = ABAQUSCore::encodeSolidId(id, sideNum);
// Further checks?
connectivity_.append(std::move(elemNodes));
elemTypes_.append(shape);
elemIds_.append(newElemId);
elsetIds_.append(setId);
}
return (elemTypes_.size() - initialCount);
}
void Foam::fileFormats::ABAQUSCore::readHelper::read
(
ISstream& is
@ -427,9 +551,11 @@ void Foam::fileFormats::ABAQUSCore::readHelper::read
// Some abaqus files use upper-case or mixed-case for section names,
// convert all to upper-case for ease.
string upperLine(stringOps::upper(line));
const string upperLine(stringOps::upper(line));
//
// "*Nodes" section
//
if (upperLine.starts_with("*NODE"))
{
// Ignore "NSET=...", we cannot do anything useful with it
@ -446,14 +572,16 @@ void Foam::fileFormats::ABAQUSCore::readHelper::read
continue;
}
//
// "*Element" section
//
if (upperLine.starts_with("*ELEMENT,"))
{
// Must have "TYPE=..."
auto elemTypeName = getIdentifier("TYPE", line);
const string elemTypeName(getIdentifier("TYPE", line));
// May have "ELSET=..." on the same line
string elsetName(getIdentifier("ELSET", line));
const string elsetName(getIdentifier("ELSET", line));
const shapeType shape(getElementType(elemTypeName));
@ -485,18 +613,42 @@ void Foam::fileFormats::ABAQUSCore::readHelper::read
continue;
}
//
// "*Surface" section
//
if (upperLine.starts_with("*SURFACE,"))
{
#ifdef Foam_readAbaqusSurface
// Require "NAME=..." on the same line
const string elsetName(getIdentifier("NAME", line));
// May have "TYPE=..." on the same line.
// If missing, default is ELEMENT.
const string surfTypeName(getIdentifier("TYPE", line));
if
(
!surfTypeName.empty()
&& stringOps::upper(surfTypeName) != "ELEMENT"
)
{
Info<< "Reading abaqus surface type "
<< surfTypeName << " is not implemented" << nl;
continue;
}
// Treat like an element set
const label elsetId = addNewElset(elsetName);
skipComments(is);
#else
Info<< "Reading of abaqus surfaces not implemented" << nl;
#endif
nread = readSurfaceElements(is, elsetId);
if (verbose_)
{
InfoErr
<< "Read " << nread << " *SURFACE entries for "
<< elsetName << nl;
}
continue;
}
}
@ -623,6 +775,15 @@ void Foam::fileFormats::ABAQUSCore::readHelper::compact_nodes()
}
void Foam::fileFormats::ABAQUSCore::readHelper::renumber_elements_1to0()
{
for (label& elemId : elemIds_)
{
renumber0_elemId(elemId);
}
}
void Foam::fileFormats::ABAQUSCore::writePoints
(
Ostream& os,

View File

@ -90,7 +90,6 @@ SourceFiles
namespace Foam
{
namespace fileFormats
{
@ -148,11 +147,62 @@ public:
return (tag & 0x80);
}
//- Is a combined (synthetic) face Id?
inline static bool isEncodedSolidId(const label combinedId)
{
return (combinedId < 0);
}
//- Combine solid element Id and side Id into synthetic face Id
inline static label encodeSolidId(const label id, const label side)
{
return -(10 * id + side);
}
//- Entangle solid element id from synthetic face Id
inline static label decodeSolidElementId(const label combinedId)
{
return
(
isEncodedSolidId(combinedId)
? (-combinedId / 10)
: combinedId
);
}
//- Entangle solid side id from synthetic face Id
//- Synthesize faceId from solid element Id and sideId
inline static label decodeSolidSideNum(const label combinedId)
{
return
(
isEncodedSolidId(combinedId)
? (-combinedId % 10)
: 0
);
}
protected:
// Protected Member Functions
//- From 1-based to 0-based
inline static void renumber0_elemId(label& elemId)
{
if (isEncodedSolidId(elemId))
{
// Eg,
// 1-based (solid 4, side 2) is -42
// 0-based (solid 3, side 2) is -32
elemId += 10;
}
else
{
--elemId;
}
}
//- Face addressing from ABAQUS faces to OpenFOAM faces.
// For hex, prism, tet primitive shapes.
static const Map<labelList>& abaqusToFoamFaceAddr();
@ -249,12 +299,26 @@ protected:
const label setId = 0
);
//- Read elements within an "*Surface" section.
// If the shape is known/supported, appends to
// connectivity, elemType, elemIds lists.
//
// \return the number of elements read
label readSurfaceElements
(
ISstream& is,
const label setId = 0
);
//- Remove non-shell elements and compact the points
void purge_solids();
//- Compact unused points and relabel connectivity
void compact_nodes();
//- Renumber elements from 1-based to 0-based
void renumber_elements_1to0();
};

View File

@ -32,7 +32,7 @@ License
#include "coupledPolyPatch.H"
#include "sampledSurface.H"
#include "mergePoints.H"
#include "indirectPrimitivePatch.H"
#include "uindirectPrimitivePatch.H"
#include "PatchTools.H"
#include "addToRunTimeSelectionTable.H"
@ -130,115 +130,225 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::obr() const
void Foam::functionObjects::fieldValues::surfaceFieldValue::setFaceZoneFaces()
{
const label zoneId = mesh_.faceZones().findZoneID(regionName_);
// Indices for all matches, already sorted
const labelList zoneIds
(
mesh_.faceZones().indices(selectionNames_)
);
if (zoneId < 0)
// Total number of faces selected
label numFaces = 0;
for (const label zoneId : zoneIds)
{
FatalErrorInFunction
<< type() << " " << name() << ": "
<< regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
<< " Unknown face zone name: " << regionName_
<< ". Valid face zones are: " << mesh_.faceZones().names()
<< nl << exit(FatalError);
numFaces += mesh_.faceZones()[zoneId].size();
}
const faceZone& fZone = mesh_.faceZones()[zoneId];
DynamicList<label> faceIds(fZone.size());
DynamicList<label> facePatchIds(fZone.size());
DynamicList<bool> faceFlip(fZone.size());
forAll(fZone, i)
if (zoneIds.empty())
{
const label facei = fZone[i];
FatalErrorInFunction
<< type() << ' ' << name() << ": "
<< regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
<< " No matching face zone(s): "
<< flatOutput(selectionNames_) << nl
<< " Known face zones: "
<< flatOutput(mesh_.faceZones().names()) << nl
<< exit(FatalError);
}
label faceId = -1;
label facePatchId = -1;
if (mesh_.isInternalFace(facei))
// Could also check this
#if 0
if (!returnReduce(bool(numFaces), orOp<bool>()))
{
WarningInFunction
<< type() << ' ' << name() << ": "
<< regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
<< " The faceZone specification: "
<< flatOutput(selectionNames_) << nl
<< " resulted in 0 faces" << nl
<< exit(FatalError);
}
#endif
faceId_.resize(numFaces);
facePatchId_.resize(numFaces);
faceFlip_.resize(numFaces);
numFaces = 0;
for (const label zoneId : zoneIds)
{
const faceZone& fZone = mesh_.faceZones()[zoneId];
forAll(fZone, i)
{
faceId = facei;
facePatchId = -1;
}
else
{
facePatchId = mesh_.boundaryMesh().whichPatch(facei);
const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
if (isA<coupledPolyPatch>(pp))
const label meshFacei = fZone[i];
const bool isFlip = fZone.flipMap()[i];
// Internal faces
label faceId = meshFacei;
label facePatchId = -1;
// Boundary faces
if (!mesh_.isInternalFace(meshFacei))
{
if (refCast<const coupledPolyPatch>(pp).owner())
facePatchId = mesh_.boundaryMesh().whichPatch(meshFacei);
const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
if (isA<coupledPolyPatch>(pp))
{
faceId = pp.whichFace(facei);
if (refCast<const coupledPolyPatch>(pp).owner())
{
faceId = pp.whichFace(meshFacei);
}
else
{
faceId = -1;
}
}
else if (!isA<emptyPolyPatch>(pp))
{
faceId = meshFacei - pp.start();
}
else
{
faceId = -1;
facePatchId = -1;
}
}
else if (!isA<emptyPolyPatch>(pp))
{
faceId = facei - pp.start();
}
else
{
faceId = -1;
facePatchId = -1;
}
}
if (faceId >= 0)
{
faceIds.append(faceId);
facePatchIds.append(facePatchId);
faceFlip.append(fZone.flipMap()[i] ? true : false);
if (faceId >= 0)
{
faceId_[numFaces] = faceId;
facePatchId_[numFaces] = facePatchId;
faceFlip_[numFaces] = isFlip;
++numFaces;
}
}
}
faceId_.transfer(faceIds);
facePatchId_.transfer(facePatchIds);
faceFlip_.transfer(faceFlip);
// Shrink to size used
faceId_.resize(numFaces);
facePatchId_.resize(numFaces);
faceFlip_.resize(numFaces);
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_);
// Patch indices for all matches
labelList patchIds;
if (patchid < 0)
// Total number of faces selected
label numFaces = 0;
labelList selected
(
mesh_.boundaryMesh().patchSet
(
selectionNames_,
false // warnNotFound - we do that ourselves
).sortedToc()
);
DynamicList<label> bad;
for (const label patchi : selected)
{
const polyPatch& pp = mesh_.boundaryMesh()[patchi];
if (isA<emptyPolyPatch>(pp))
{
bad.append(patchi);
}
else
{
numFaces += pp.size();
}
}
if (bad.size())
{
label nGood = (selected.size() - bad.size());
auto& os = (nGood > 0 ? WarningInFunction : FatalErrorInFunction);
os << "Cannot sample an empty patch" << nl;
for (const label patchi : bad)
{
os << " "
<< mesh_.boundaryMesh()[patchi].name() << nl;
}
if (nGood)
{
os << "No non-empty patches selected" << endl
<< exit(FatalError);
}
else
{
os << "Selected " << nGood << " non-empty patches" << nl;
}
patchIds.resize(nGood);
nGood = 0;
for (const label patchi : selected)
{
if (!bad.found(patchi))
{
patchIds[nGood] = patchi;
++nGood;
}
}
}
else
{
patchIds = std::move(selected);
}
if (patchIds.empty())
{
FatalErrorInFunction
<< type() << " " << name() << ": "
<< type() << ' ' << name() << ": "
<< regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
<< " Unknown patch name: " << regionName_
<< ". Valid patch names are: "
<< " No matching patch name(s): "
<< flatOutput(selectionNames_) << nl
<< " Known patch names:" << nl
<< mesh_.boundaryMesh().names() << nl
<< exit(FatalError);
}
const polyPatch& pp = mesh_.boundaryMesh()[patchid];
label nFaces = pp.size();
if (isA<emptyPolyPatch>(pp))
// Could also check this
#if 0
if (!returnReduce(bool(numFaces), orOp<bool>()))
{
nFaces = 0;
WarningInFunction
<< type() << ' ' << name() << ": "
<< regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
<< " The patch specification: "
<< flatOutput(selectionNames_) << nl
<< " resulted in 0 faces" << nl
<< exit(FatalError);
}
#endif
faceId_.setSize(nFaces);
facePatchId_.setSize(nFaces);
faceFlip_.setSize(nFaces);
faceId_.resize(numFaces);
facePatchId_.resize(numFaces);
faceFlip_.resize(numFaces);
nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
forAll(faceId_, facei)
numFaces = 0;
for (const label patchi : patchIds)
{
faceId_[facei] = facei;
facePatchId_[facei] = patchid;
faceFlip_[facei] = false;
const polyPatch& pp = mesh_.boundaryMesh()[patchi];
const label len = pp.size();
SubList<label>(faceId_, len, numFaces) = identity(len);
SubList<label>(facePatchId_, len, numFaces) = patchi;
SubList<bool>(faceFlip_, len, numFaces) = false;
numFaces += len;
}
}
@ -252,24 +362,25 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry
List<faceList> allFaces(Pstream::nProcs());
List<pointField> allPoints(Pstream::nProcs());
labelList globalFacesIs(faceId_);
forAll(globalFacesIs, i)
{
if (facePatchId_[i] != -1)
IndirectList<face> selectedFaces(mesh_.faces(), labelList(faceId_));
labelList& meshFaceIds = selectedFaces.addressing();
forAll(meshFaceIds, i)
{
const label patchi = facePatchId_[i];
globalFacesIs[i] += mesh_.boundaryMesh()[patchi].start();
if (patchi != -1)
{
meshFaceIds[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();
// Add local faces and points to the all* lists
uindirectPrimitivePatch pp(selectedFaces, mesh_.points());
allFaces[Pstream::myProcNo()] = pp.localFaces();
allPoints[Pstream::myProcNo()] = pp.localPoints();
}
Pstream::gatherList(allFaces);
Pstream::gatherList(allPoints);
@ -283,53 +394,41 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry
nPoints += allPoints[proci].size();
}
faces.setSize(nFaces);
points.setSize(nPoints);
faces.resize(nFaces);
points.resize(nPoints);
nFaces = 0;
nPoints = 0;
// My own data first
// My data first
{
const faceList& fcs = allFaces[Pstream::myProcNo()];
for (const face& f : fcs)
for (const face& f : allFaces[Pstream::myProcNo()])
{
face& newF = faces[nFaces++];
newF.setSize(f.size());
forAll(f, fp)
{
newF[fp] = f[fp] + nPoints;
}
faces[nFaces++] = offsetOp<face>()(f, nPoints);
}
const pointField& pts = allPoints[Pstream::myProcNo()];
for (const point& pt : pts)
for (const point& p : allPoints[Pstream::myProcNo()])
{
points[nPoints++] = pt;
points[nPoints++] = p;
}
}
// Other proc data follows
forAll(allFaces, proci)
{
if (proci != Pstream::myProcNo())
if (proci == Pstream::myProcNo())
{
const faceList& fcs = allFaces[proci];
for (const face& f : fcs)
{
face& newF = faces[nFaces++];
newF.setSize(f.size());
forAll(f, fp)
{
newF[fp] = f[fp] + nPoints;
}
}
continue;
}
const pointField& pts = allPoints[proci];
for (const point& pt : pts)
{
points[nPoints++] = pt;
}
for (const face& f : allFaces[proci])
{
faces[nFaces++] = offsetOp<face>()(f, nPoints);
}
for (const point& p : allPoints[proci])
{
points[nPoints++] = p;
}
}
@ -383,11 +482,7 @@ combineSurfaceGeometry
PatchTools::gatherAndMerge
(
mergeDim,
primitivePatch
(
SubList<face>(s.faces(), s.faces().size()),
s.points()
),
primitivePatch(SubList<face>(s.faces()), s.points()),
points,
faces,
pointsMap
@ -413,11 +508,7 @@ combineSurfaceGeometry
PatchTools::gatherAndMerge
(
mergeDim,
primitivePatch
(
SubList<face>(s.faces(), s.faces().size()),
s.points()
),
primitivePatch(SubList<face>(s.faces()), s.points()),
points,
faces,
pointsMap
@ -522,7 +613,7 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::update()
if (nFaces_ == 0)
{
FatalErrorInFunction
<< type() << " " << name() << ": "
<< type() << ' ' << name() << ": "
<< regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
<< " Region has no faces" << exit(FatalError);
}
@ -548,15 +639,20 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader
if (canWriteHeader() && (operation_ != opNone))
{
writeCommented(os, "Region type : ");
os << regionTypeNames_[regionType_] << " " << regionName_ << endl;
os << regionTypeNames_[regionType_] << ' ' << regionName_ << nl;
writeHeaderValue(os, "Faces", nFaces_);
writeHeaderValue(os, "Area", totalArea_);
writeHeaderValue(os, "Scale factor", scaleFactor_);
if (weightFieldName_ != "none")
if (weightFieldNames_.size())
{
writeHeaderValue(os, "Weight field", weightFieldName_);
writeHeaderValue
(
os,
"Weight field",
flatOutput(weightFieldNames_, FlatOutput::BareComma{})
);
}
writeCommented(os, "Time");
@ -570,7 +666,7 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader
for (const word& fieldName : fields_)
{
os << tab << operationTypeNames_[operation_]
<< "(" << fieldName << ")";
<< '(' << fieldName << ')';
}
os << endl;
@ -821,9 +917,10 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
true // Failsafe behaviour
)
),
weightFieldName_("none"),
needsUpdate_(true),
writeArea_(false),
selectionNames_(),
weightFieldNames_(),
totalArea_(0),
nFaces_(0),
faceId_(),
@ -854,9 +951,10 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
true // Failsafe behaviour
)
),
weightFieldName_("none"),
needsUpdate_(true),
writeArea_(false),
selectionNames_(),
weightFieldNames_(),
totalArea_(0),
nFaces_(0),
faceId_(),
@ -876,18 +974,55 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
{
fieldValue::read(dict);
weightFieldName_ = "none";
needsUpdate_ = true;
writeArea_ = dict.getOrDefault("writeArea", false);
weightFieldNames_.clear();
totalArea_ = 0;
nFaces_ = 0;
faceId_.clear();
facePatchId_.clear();
faceFlip_.clear();
sampledPtr_.clear();
surfaceWriterPtr_.clear();
sampledPtr_.reset(nullptr);
surfaceWriterPtr_.reset(nullptr);
// Can have "name" (word) and/or "names" (wordRes)
//
// If "names" exists AND contains a literal (non-regex) that can be used
// as a suitable value for "name", the "name" entry becomes optional.
regionName_.clear();
selectionNames_.clear();
{
dict.readIfPresent("names", selectionNames_);
for (const auto& item : selectionNames_)
{
if (item.isLiteral())
{
regionName_ = item;
break;
}
}
// Mandatory if we didn't pick up a value from selectionNames_
dict.readEntry
(
"name",
regionName_,
keyType::LITERAL,
regionName_.empty()
);
// Ensure there is always content for selectionNames_
if (selectionNames_.empty())
{
selectionNames_.resize(1);
selectionNames_.first() = regionName_;
}
}
dict.readEntry("name", regionName_);
// Create sampled surface, but leave 'expired' (ie, no update) since it
// may depend on fields or data that do not yet exist
@ -901,7 +1036,7 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
);
}
Info<< type() << " " << name() << ":" << nl
Info<< type() << ' ' << name() << ':' << nl
<< " operation = ";
if (postOperation_ != postOpNone)
@ -925,11 +1060,29 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
<< exit(FatalIOError);
}
if (dict.readIfPresent("weightField", weightFieldName_))
// Can have "weightFields" or "weightField"
bool missing = true;
if (dict.readIfPresent("weightFields", weightFieldNames_))
{
Info<< " weight field = " << weightFieldName_ << nl;
missing = false;
}
else
{
weightFieldNames_.resize(1);
if (dict.readIfPresent("weightField", weightFieldNames_.first()))
{
missing = false;
if ("none" == weightFieldNames_.first())
{
// "none" == no weighting
weightFieldNames_.clear();
}
}
}
if (missing)
{
// Suggest possible alternative unweighted operation?
FatalIOErrorInFunction(dict)
@ -940,6 +1093,16 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
<< "or use a different operation."
<< exit(FatalIOError);
}
Info<< " weight field = ";
if (weightFieldNames_.empty())
{
Info<< "none" << nl;
}
else
{
Info<< flatOutput(weightFieldNames_) << nl;
}
}
// Backwards compatibility for v1612 and older
@ -1048,46 +1211,79 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::write()
}
}
// Only a few weight types (scalar, vector)
if (weightFieldName_ != "none")
// Check availability and type of weight field
// Only support a few weight types:
// scalar: 0-N fields
// vector: 0-1 fields
// Default is a zero-size scalar weight field (ie, weight = 1)
scalarField scalarWeights;
vectorField vectorWeights;
for (const word& weightName : weightFieldNames_)
{
if (validField<scalar>(weightFieldName_))
if (validField<scalar>(weightName))
{
scalarField weightField
(
getFieldValues<scalar>(weightFieldName_, true)
);
tmp<scalarField> tfld = getFieldValues<scalar>(weightName, true);
// Process the fields
writeAll(Sf, weightField, points, faces);
if (scalarWeights.empty())
{
scalarWeights = tfld;
}
else
{
scalarWeights *= tfld;
}
}
else if (validField<vector>(weightFieldName_))
else if (validField<vector>(weightName))
{
vectorField weightField
(
getFieldValues<vector>(weightFieldName_, true)
);
tmp<vectorField> tfld = getFieldValues<vector>(weightName, true);
// Process the fields
writeAll(Sf, weightField, points, faces);
if (vectorWeights.empty())
{
vectorWeights = tfld;
}
else
{
FatalErrorInFunction
<< "weightField " << weightName
<< " - only one vector weight field allowed. " << nl
<< "weights: " << flatOutput(weightFieldNames_) << nl
<< abort(FatalError);
}
}
else
else if (weightName != "none")
{
// Silently ignore "none", flag everything else as an error
// TBD: treat missing "rho" like incompressible with rho = 1
// and/or provided rhoRef value
FatalErrorInFunction
<< "weightField " << weightFieldName_
<< " not found or an unsupported type"
<< "weightField " << weightName
<< " not found or an unsupported type" << nl
<< abort(FatalError);
}
}
// Process the fields
if (vectorWeights.size())
{
if (scalarWeights.size())
{
vectorWeights *= scalarWeights;
}
writeAll(Sf, vectorWeights, points, faces);
}
else
{
// Default is a zero-size scalar weight field (ie, weight = 1)
scalarField weightField;
// Process the fields
writeAll(Sf, weightField, points, faces);
writeAll(Sf, scalarWeights, points, faces);
}
if (operation_ != opNone)
{
file() << endl;

View File

@ -31,17 +31,17 @@ Group
grpFieldFunctionObjects
Description
Provides a 'face regionType' variant of the \c fieldValues function object.
A \c face regionType variant of the \c fieldValues function object.
Given a list of user-specified fields and a selection of mesh (or general
surface) faces, a number of operations can be performed, such as sums,
averages and integrations.
For example, to calculate the volumetric or mass flux across a patch,
apply the 'sum' operator to the flux field (typically \c phi).
apply the \c sum operator to the flux field (typically \c phi).
Usage
Minimal example by using \c system/controlDict.functions:
A minimal example:
\verbatim
surfaceFieldValuePatch1
{
@ -56,6 +56,8 @@ Usage
name <patch>;
// Optional entries (runtime modifiable)
names (<patch-name> <patch-regex>);
postOperation none;
weightField alpha1;
scaleFactor 1.0;
@ -79,8 +81,10 @@ Usage
name <faceZone>;
// Optional entries (runtime modifiable)
names (<zone-name> <zone-regex>);
postOperation none;
weightField alpha1;
weightFields (rho U);
scaleFactor 1.0;
writeArea false;
surfaceFormat none;
@ -92,15 +96,17 @@ Usage
where the entries mean:
\table
Property | Description | Type | Req'd | Dflt
Property | Description | Type | Reqd | Dflt
type | Type name: surfaceFieldValue | word | yes | -
libs | Library name: fieldFunctionObjects | word | yes | -
fields | Names of operand fields | wordList | yes | -
libs | Libraries: fieldFunctionObjects | wordList | yes | -
regionType | Face regionType: see below | word | yes | -
name | Name for regionType | word | yes | -
fields | Names of operand fields | wordList | yes | -
name | Name of the regionType | word | yes | -
names | Extended selection | word/regex list | no | -
operation | Operation type: see below | word | yes | -
postOperation | Post-operation type: see below | word | no | none
weightField | Name of field to apply weighting | word | no | none
weightField | Name of field to apply weighting | word | maybe |
weightFields | Names of fields to apply weighting | wordList | maybe |
scaleFactor | Output value scaling factor | scalar | no | 1.0
writeArea | Write the surface area | bool | no | false
surfaceFormat | Output value format | word <!--
@ -112,9 +118,9 @@ Usage
Options for the \c regionType entry:
\plaintable
faceZone | The \b name entry to specify the faceZone
patch | The \b name entry to specify the patch
functionObjectSurface | The \b name entry to specify a polySurface
faceZone | The \b name entry specifies a faceZone. Supports \b names
patch | The \b name entry specifies a patch. Supports \b names
functionObjectSurface | The \b name entry specifies a polySurface
sampledSurface | A \b sampledSurfaceDict sub-dictionary and \b name
\endplaintable
@ -156,12 +162,16 @@ Usage
Usage by the \c postProcess utility is not available.
Note
- Some types (eg, patch or faceZone) support the selection of multiple
entries, which can be specified as list of words or regular expressions
by \b names entry.<br>
If the \b names enty exists \em and contains a literal that can be used
as a suitable value for \b name, the \b name entry becomes optional
instead of being mandatory.
- The values reported by the \c areaNormalAverage and \c areaNormalIntegrate
operations are written as the first component of a field with the same
rank as the input field.
- Faces on empty patches get ignored.
- If the field is a volField the \c faceZone
can only consist of boundary faces.
- Faces on empty patches are ignored.
- Using \c functionObjectSurface:
- The keyword %subRegion should not be used to select surfaces.
Instead specify the regionType 'functionObjectSurface' and provide
@ -242,8 +252,8 @@ public:
//- Region type enumeration
enum regionTypes
{
stFaceZone = 0x01, //!< Calculate on a faceZone
stPatch = 0x02, //!< Calculate on a patch
stFaceZone = 0x01, //!< Calculate with faceZone(s)
stPatch = 0x02, //!< Calculate with patch(es)
stObject = 0x11, //!< Calculate with function object surface
stSampled = 0x12 //!< Sample onto surface and calculate
};
@ -385,15 +395,18 @@ protected:
//- Optional post-evaluation operation
postOperationType postOperation_;
//- Weight field name - optional
word weightFieldName_;
//- Track if the surface needs an update
bool needsUpdate_;
//- Optionally write the area of the surfaceFieldValue
bool writeArea_;
//- Extended selections
wordRes selectionNames_;
//- Weight field name(s) - optional
wordList weightFieldNames_;
//- Total area of the surfaceFieldValue
scalar totalArea_;
@ -401,17 +414,20 @@ protected:
label nFaces_;
// If operating on mesh faces (faceZone, patch)
// If operating on mesh faces (faceZone, patch)
//- Local list of face IDs
labelList faceId_;
//- Local list of face IDs
labelList faceId_;
//- Local list of patch ID per face
labelList facePatchId_;
//- Local list of patch ID per face
labelList facePatchId_;
//- List representing the face flip map
// (false: use as-is, true: negate)
boolList faceFlip_;
//- List representing the face flip map
// (false: use as-is, true: negate)
boolList faceFlip_;
// Demand-driven
//- The sampledSurface (when operating on sampledSurface)
autoPtr<sampledSurface> sampledPtr_;

View File

@ -118,7 +118,7 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::getFieldValues
if (mandatory)
{
FatalErrorInFunction
<< "Field " << fieldName << " not found in database"
<< "Field " << fieldName << " not found in database" << nl
<< abort(FatalError);
}

View File

@ -135,9 +135,15 @@ void Foam::functionObjects::fieldValues::volFieldValue::writeFileHeader
) const
{
volRegion::writeFileHeader(*this, os);
if (weightFieldName_ != "none")
if (weightFieldNames_.size())
{
writeHeaderValue(os, "Weight field", weightFieldName_);
writeHeaderValue
(
os,
"Weight field",
flatOutput(weightFieldNames_, FlatOutput::BareComma{})
);
}
writeCommented(os, "Time");
@ -210,7 +216,7 @@ Foam::functionObjects::fieldValues::volFieldValue::volFieldValue
true // Failsafe behaviour
)
),
weightFieldName_("none")
weightFieldNames_()
{
read(dict);
writeFileHeader(file());
@ -237,7 +243,7 @@ Foam::functionObjects::fieldValues::volFieldValue::volFieldValue
true // Failsafe behaviour
)
),
weightFieldName_("none")
weightFieldNames_()
{
read(dict);
}
@ -252,15 +258,33 @@ bool Foam::functionObjects::fieldValues::volFieldValue::read
{
fieldValue::read(dict);
weightFieldName_ = "none";
weightFieldNames_.clear();
if (usesWeight())
{
if (dict.readIfPresent("weightField", weightFieldName_))
// Can have "weightFields" or "weightField"
bool missing = true;
if (dict.readIfPresent("weightFields", weightFieldNames_))
{
Info<< " weight field = " << weightFieldName_;
missing = false;
}
else
{
weightFieldNames_.resize(1);
if (dict.readIfPresent("weightField", weightFieldNames_.first()))
{
missing = false;
if ("none" == weightFieldNames_.first())
{
// "none" == no weighting
weightFieldNames_.clear();
}
}
}
if (missing)
{
// Suggest possible alternative unweighted operation?
FatalIOErrorInFunction(dict)
@ -271,6 +295,16 @@ bool Foam::functionObjects::fieldValues::volFieldValue::read
<< "or use a different operation."
<< exit(FatalIOError);
}
Info<< " weight field = ";
if (weightFieldNames_.empty())
{
Info<< "none" << nl;
}
else
{
Info<< flatOutput(weightFieldNames_) << nl;
}
}
Info<< nl << endl;
@ -297,14 +331,45 @@ bool Foam::functionObjects::fieldValues::volFieldValue::write()
V = filterField(fieldValue::mesh_.V());
}
// Weight field - zero-size means weight = 1
scalarField weightField;
if (weightFieldName_ != "none")
// Check availability and type of weight field
// Only support a few weight types:
// scalar: 0-N fields
// Default is a zero-size scalar weight field (ie, weight = 1)
scalarField scalarWeights;
for (const word& weightName : weightFieldNames_)
{
weightField = getFieldValues<scalar>(weightFieldName_, true);
if (validField<scalar>(weightName))
{
tmp<scalarField> tfld = getFieldValues<scalar>(weightName, true);
if (scalarWeights.empty())
{
scalarWeights = tfld;
}
else
{
scalarWeights *= tfld;
}
}
else if (weightName != "none")
{
// Silently ignore "none", flag everything else as an error
// TBD: treat missing "rho" like incompressible with rho = 1
// and/or provided rhoRef value
FatalErrorInFunction
<< "weightField " << weightName
<< " not found or an unsupported type" << nl
<< abort(FatalError);
}
}
writeAll(V, weightField);
// Process the fields
writeAll(V, scalarWeights);
if (Pstream::master())
{

View File

@ -69,7 +69,8 @@ Usage
name | Name for regionType | word | yes | -
operation | Operation type: see below | word | yes | -
postOperation | Post-operation type: see below | word | no | none
weightField | Name of field to apply weighting | word | no | none
weightField | Name of field to apply weighting | word | maybe |
weightFields | Names of fields to apply weighting | wordList | maybe |
\endtable
The inherited entries are elaborated in:
@ -214,8 +215,8 @@ protected:
//- Optional post-evaluation operation
postOperationType postOperation_;
//- Weight field name - only used for weighted modes
word weightFieldName_;
//- Weight field name(s) - optional
wordList weightFieldNames_;
// Protected Member Functions
@ -233,7 +234,7 @@ protected:
// Checks for availability on any processor.
inline bool canWeight(const scalarField& weightField) const;
//- True if the field name is valid (exists, and a supported type)
//- Return true if the field name is valid
template<class Type>
bool validField(const word& fieldName) const;

View File

@ -71,7 +71,7 @@ Foam::functionObjects::fieldValues::volFieldValue::getFieldValues
if (mandatory)
{
FatalErrorInFunction
<< "Field " << fieldName << " not found in database"
<< "Field " << fieldName << " not found in database" << nl
<< abort(FatalError);
}

View File

@ -35,6 +35,7 @@ surface/isoSurface/isoSurfaceTopo.C
surface/thresholdCellFaces/thresholdCellFaces.C
sampledSurface/sampledNone/sampledNone.C
sampledSurface/sampledFaceZone/sampledFaceZone.C
sampledSurface/sampledPatch/sampledPatch.C
sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C
sampledSurface/sampledPlane/sampledPlane.C

View File

@ -0,0 +1,402 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "sampledFaceZone.H"
#include "dictionary.H"
#include "polyMesh.H"
#include "polyPatch.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "volPointInterpolation.H"
#include "uindirectPrimitivePatch.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(sampledFaceZone, 0);
addNamedToRunTimeSelectionTable
(
sampledSurface,
sampledFaceZone,
word,
faceZone
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sampledFaceZone::sampledFaceZone
(
const word& name,
const polyMesh& mesh,
const UList<wordRe>& zoneNames,
const bool triangulate
)
:
sampledSurface(name, mesh),
selectionNames_(zoneNames),
triangulate_(triangulate),
needsUpdate_(true)
{}
Foam::sampledFaceZone::sampledFaceZone
(
const word& name,
const polyMesh& mesh,
const dictionary& dict
)
:
sampledSurface(name, mesh, dict),
selectionNames_(dict.get<wordRes>("zones")),
triangulate_(dict.getOrDefault("triangulate", false)),
needsUpdate_(true)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::labelList& Foam::sampledFaceZone::zoneIDs() const
{
if (zoneIds_.empty())
{
// Zone indices for all matches, already sorted
zoneIds_ = mesh().faceZones().indices(selectionNames_);
}
return zoneIds_;
}
bool Foam::sampledFaceZone::needsUpdate() const
{
return needsUpdate_;
}
bool Foam::sampledFaceZone::expire()
{
// already marked as expired
if (needsUpdate_)
{
return false;
}
sampledSurface::clearGeom();
Mesh::clear();
zoneIds_.clear();
faceId_.clear();
facePatchId_.clear();
needsUpdate_ = true;
return true;
}
bool Foam::sampledFaceZone::update()
{
if (!needsUpdate_)
{
return false;
}
// Total number of faces selected
label numFaces = 0;
for (const label zonei : zoneIDs())
{
numFaces += mesh().faceZones()[zonei].size();
}
if (zoneIDs().empty())
{
WarningInFunction
<< type() << ' ' << name() << ": "
<< " No matching face zone(s): "
<< flatOutput(selectionNames_) << nl
<< " Known face zones: "
<< flatOutput(mesh().faceZones().names()) << nl;
}
// Could also check numFaces
// The mesh face or local patch face and the patch id
faceId_.resize(numFaces);
facePatchId_.resize(numFaces);
IndirectList<face> selectedFaces(mesh().faces(), labelList());
labelList& meshFaceIds = selectedFaces.addressing();
meshFaceIds.resize(numFaces);
numFaces = 0;
forAll(zoneIDs(), idx)
{
const label zonei = zoneIDs()[idx];
const faceZone& fZone = mesh().faceZones()[zonei];
for (const label meshFacei : fZone)
{
// Internal faces
label faceId = meshFacei;
label facePatchId = -1;
// Boundary faces
if (!mesh().isInternalFace(meshFacei))
{
facePatchId = mesh().boundaryMesh().whichPatch(meshFacei);
const polyPatch& pp = mesh().boundaryMesh()[facePatchId];
if (isA<emptyPolyPatch>(pp))
{
// Do not sample an empty patch
continue;
}
const auto* procPatch = isA<processorPolyPatch>(pp);
if (procPatch && !procPatch->owner())
{
// Do not sample neighbour-side, retain owner-side only
continue;
}
const auto* cpp = isA<coupledPolyPatch>(pp);
if (cpp)
{
faceId = (cpp->owner() ? pp.whichFace(meshFacei) : -1);
}
else
{
faceId = meshFacei - pp.start();
}
}
if (faceId >= 0)
{
faceId_[numFaces] = faceId;
facePatchId_[numFaces] = facePatchId;
meshFaceIds[numFaces] = meshFacei;
++numFaces;
}
}
}
// Shrink to size used
faceId_.resize(numFaces);
facePatchId_.resize(numFaces);
meshFaceIds.resize(numFaces);
uindirectPrimitivePatch zoneFaces(selectedFaces, mesh().points());
this->storedPoints() = zoneFaces.localPoints();
this->storedFaces() = zoneFaces.localFaces();
// triangulate - uses remapFaces()
if (triangulate_)
{
Mesh::triangulate();
}
needsUpdate_ = false;
return true;
}
// remap action on triangulation
void Foam::sampledFaceZone::remapFaces(const labelUList& faceMap)
{
if (!faceMap.empty())
{
Mesh::remapFaces(faceMap);
faceId_ = labelList
(
labelUIndList(faceId_, faceMap)
);
facePatchId_ = labelList
(
labelUIndList(facePatchId_, faceMap)
);
}
}
Foam::tmp<Foam::scalarField> Foam::sampledFaceZone::sample
(
const interpolation<scalar>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::vectorField> Foam::sampledFaceZone::sample
(
const interpolation<vector>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::sphericalTensorField> Foam::sampledFaceZone::sample
(
const interpolation<sphericalTensor>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::symmTensorField> Foam::sampledFaceZone::sample
(
const interpolation<symmTensor>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::tensorField> Foam::sampledFaceZone::sample
(
const interpolation<tensor>& sampler
) const
{
return sampleOnFaces(sampler);
}
bool Foam::sampledFaceZone::withSurfaceFields() const
{
return true;
}
Foam::tmp<Foam::scalarField> Foam::sampledFaceZone::sample
(
const surfaceScalarField& sField
) const
{
return sampleOnFaces(sField);
}
Foam::tmp<Foam::vectorField> Foam::sampledFaceZone::sample
(
const surfaceVectorField& sField
) const
{
return sampleOnFaces(sField);
}
Foam::tmp<Foam::sphericalTensorField> Foam::sampledFaceZone::sample
(
const surfaceSphericalTensorField& sField
) const
{
return sampleOnFaces(sField);
}
Foam::tmp<Foam::symmTensorField> Foam::sampledFaceZone::sample
(
const surfaceSymmTensorField& sField
) const
{
return sampleOnFaces(sField);
}
Foam::tmp<Foam::tensorField> Foam::sampledFaceZone::sample
(
const surfaceTensorField& sField
) const
{
return sampleOnFaces(sField);
}
Foam::tmp<Foam::scalarField> Foam::sampledFaceZone::interpolate
(
const interpolation<scalar>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::vectorField> Foam::sampledFaceZone::interpolate
(
const interpolation<vector>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::sphericalTensorField>
Foam::sampledFaceZone::interpolate
(
const interpolation<sphericalTensor>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::symmTensorField> Foam::sampledFaceZone::interpolate
(
const interpolation<symmTensor>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::tensorField> Foam::sampledFaceZone::interpolate
(
const interpolation<tensor>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
void Foam::sampledFaceZone::print(Ostream& os) const
{
os << "faceZone: " << name() << " :"
<< " zones: " << flatOutput(selectionNames_)
<< " faces:" << faces().size()
<< " points:" << points().size();
}
// ************************************************************************* //

View File

@ -0,0 +1,345 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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::sampledFaceZone
Description
A sampledSurface defined by the cell faces corresponding to a threshold
value.
This is often embedded as part of a sampled surfaces function object.
Usage
Example of function object partial specification:
\verbatim
surfaces
(
surface1
{
type faceZones;
zones (zone1 "sides.*");
}
);
\endverbatim
Where the sub-entries comprise:
\table
Property | Description | Required | Default
type | faceZones | yes |
zones | zone selection as word/regex list | yes |
triangulate | triangulate faces | no | false
\endtable
SourceFiles
sampledFaceZone.C
sampledFaceZoneTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef sampledFaceZone_H
#define sampledFaceZone_H
#include "sampledSurface.H"
#include "MeshedSurfaces.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class sampledFaceZone Declaration
\*---------------------------------------------------------------------------*/
class sampledFaceZone
:
public meshedSurface,
public sampledSurface
{
//- Mesh storage type
typedef meshedSurface Mesh;
// Private Data
//- Selection (word/regex) of face zones
const wordRes selectionNames_;
//- IDs for selected face zones (sorted)
mutable labelList zoneIds_;
//- Triangulated faces or keep faces as is
bool triangulate_;
//- Track if the surface needs an update
mutable bool needsUpdate_;
//- Local list of face IDs
labelList faceId_;
//- Local list of patch ID per face. Is -1 for internal face
labelList facePatchId_;
// Private Member Functions
//- Sample volume/boundary field onto surface faces
template<class Type>
tmp<Field<Type>> sampleOnFaces
(
const interpolation<Type>& sampler
) const;
//- Sample surface field onto surface faces
template<class Type>
tmp<Field<Type>> sampleOnFaces
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& sField
) const;
//- Interpolate volume/boundary field onto surface points
template<class Type>
tmp<Field<Type>> sampleOnPoints
(
const interpolation<Type>& interpolator
) const;
//- Re-map action on triangulation or cleanup
virtual void remapFaces(const labelUList& faceMap);
protected:
//- The selected face zones (sorted)
const labelList& zoneIDs() const;
public:
//- Runtime type information
TypeName("faceZone");
// Constructors
//- Construct from components
sampledFaceZone
(
const word& name,
const polyMesh& mesh,
const UList<wordRe>& zoneNames,
const bool triangulate = false
);
//- Construct from dictionary
sampledFaceZone
(
const word& name,
const polyMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~sampledFaceZone() = default;
// Member Functions
//- Does the surface need an update?
virtual bool needsUpdate() const;
//- Mark the surface as needing an update.
// May also free up unneeded data.
// Return false if surface was already marked as expired.
virtual bool expire();
//- Update the surface as required.
// Do nothing (and return false) if no update was needed
virtual bool update();
//- Points of surface
virtual const pointField& points() const
{
return Mesh::points();
}
//- Faces of surface
virtual const faceList& faces() const
{
return Mesh::surfFaces();
}
//- Per-face zone/region information
virtual const labelList& zoneIds() const
{
return labelList::null();
}
//- Face area vectors (normals)
virtual const vectorField& Sf() const
{
return Mesh::Sf();
}
//- Face area magnitudes
virtual const scalarField& magSf() const
{
return Mesh::magSf();
}
//- Face centres
virtual const vectorField& Cf() const
{
return Mesh::Cf();
}
// Sample
//- Sample volume field onto surface faces
virtual tmp<scalarField> sample
(
const interpolation<scalar>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<vectorField> sample
(
const interpolation<vector>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<sphericalTensorField> sample
(
const interpolation<sphericalTensor>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<symmTensorField> sample
(
const interpolation<symmTensor>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<tensorField> sample
(
const interpolation<tensor>& sampler
) const;
//- Can it sample surface-fields?
virtual bool withSurfaceFields() const;
//- Sample surface field on face zone
virtual tmp<scalarField> sample
(
const surfaceScalarField&
) const;
//- Sample surface field on face zone
virtual tmp<vectorField> sample
(
const surfaceVectorField&
) const;
//- Sample surface field on face zone
virtual tmp<sphericalTensorField> sample
(
const surfaceSphericalTensorField&
) const;
//- Sample surface field on face zone
virtual tmp<symmTensorField> sample
(
const surfaceSymmTensorField&
) const;
//- Sample surface field on face zone
virtual tmp<tensorField> sample
(
const surfaceTensorField&
) const;
// Interpolate
//- Interpolate volume field onto surface points
virtual tmp<scalarField> interpolate
(
const interpolation<scalar>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<vectorField> interpolate
(
const interpolation<vector>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<sphericalTensorField> interpolate
(
const interpolation<sphericalTensor>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<symmTensorField> interpolate
(
const interpolation<symmTensor>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<tensorField> interpolate
(
const interpolation<tensor>& interpolator
) const;
// Output
//- Print information
virtual void print(Ostream& os) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "sampledFaceZoneTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,153 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "sampledFaceZone.H"
#include "volFieldsFwd.H"
#include "pointFields.H"
#include "volPointInterpolation.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::sampledFaceZone::sampleOnFaces
(
const interpolation<Type>& sampler
) const
{
const auto& vField = sampler.psi();
const labelList& own = mesh().faceOwner();
const labelList& nei = mesh().faceNeighbour();
// One value per face
auto tvalues = tmp<Field<Type>>::New(faceId_.size());
auto& values = tvalues.ref();
forAll(faceId_, i)
{
const label facei = faceId_[i];
const label patchi = facePatchId_[i];
if (patchi != -1)
{
// Boundary face - face id is the patch-local face id
values[i] = vField.boundaryField()[patchi][facei];
}
else
{
// Internal face
values[i] = 0.5*(vField[own[facei]] + vField[nei[facei]]);
}
}
return tvalues;
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::sampledFaceZone::sampleOnFaces
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& sField
) const
{
// One value per face
auto tvalues = tmp<Field<Type>>::New(faceId_.size());
auto& values = tvalues.ref();
forAll(faceId_, i)
{
const label facei = faceId_[i];
const label patchi = facePatchId_[i];
if (patchi != -1)
{
// Boundary face - face id is the patch-local face id
values[i] = sField.boundaryField()[patchi][facei];
}
else
{
// Internal face
values[i] = sField[facei];
}
}
return tvalues;
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::sampledFaceZone::sampleOnPoints
(
const interpolation<Type>& interpolator
) const
{
// One value per vertex
auto tvalues = tmp<Field<Type>>::New(points().size(), Zero);
auto& values = tvalues.ref();
const labelList& own = mesh().faceOwner();
bitSet pointDone(points().size());
forAll(faces(), i)
{
const face& f = faces()[i];
label facei = faceId_[i];
const label patchi = facePatchId_[i];
if (patchi != -1)
{
// Boundary face. patch-local face id -> mesh face id
const polyPatch& pp = mesh().boundaryMesh()[patchi];
facei += pp.start();
}
const label celli = own[facei];
for (const label pointi : f)
{
if (pointDone.set(pointi))
{
values[pointi] = interpolator.interpolate
(
points()[pointi],
celli,
facei
);
}
}
}
return tvalues;
}
// ************************************************************************* //

View File

@ -170,9 +170,9 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
// Does approximation by looking at the face centres only
const pointField& fc = surface_.faceCentres();
List<nearInfo> nearest(fc.size(), nearInfo(GREAT, labelMax));
List<nearInfo> nearest(fc.size(), nearInfo(Foam::sqr(GREAT), labelMax));
if (sampleSource_ == cells)
if (sampleSource_ == samplingSource::cells)
{
// Search for nearest cell
@ -181,17 +181,18 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
forAll(fc, facei)
{
const point& pt = fc[facei];
auto& near = nearest[facei];
pointIndexHit info = cellTree.findNearest(pt, sqr(GREAT));
pointIndexHit info = cellTree.findNearest(pt, Foam::sqr(GREAT));
if (info.hit())
{
nearest[facei].first() = magSqr(info.hitPoint()-pt);
nearest[facei].second() = globalCells.toGlobal(info.index());
near.first() = magSqr(info.hitPoint()-pt);
near.second() = globalCells.toGlobal(info.index());
}
}
}
else if (sampleSource_ == insideCells)
else if (sampleSource_ == samplingSource::insideCells)
{
// Search for cell containing point
@ -200,35 +201,37 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
forAll(fc, facei)
{
const point& pt = fc[facei];
auto& near = nearest[facei];
if (cellTree.bb().contains(pt))
{
const label index = cellTree.findInside(pt);
if (index != -1)
{
nearest[facei].first() = 0;
nearest[facei].second() = globalCells.toGlobal(index);
near.first() = 0;
near.second() = globalCells.toGlobal(index);
}
}
}
}
else
else // samplingSource::boundaryFaces
{
// Search for nearest boundaryFace
// Search for nearest boundary face
// on all non-coupled boundary faces
//- Search on all non-coupled boundary faces
const auto& bndTree = meshSearcher.nonCoupledBoundaryTree();
forAll(fc, facei)
{
const point& pt = fc[facei];
auto& near = nearest[facei];
pointIndexHit info = bndTree.findNearest(pt, sqr(GREAT));
if (info.hit())
{
nearest[facei].first() = magSqr(info.hitPoint()-pt);
nearest[facei].second() =
near.first() = magSqr(info.hitPoint()-pt);
near.second() =
globalCells.toGlobal
(
bndTree.shapes().faceLabels()[info.index()]
@ -250,16 +253,26 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
forAll(nearest, facei)
{
const label index = nearest[facei].second();
const auto& near = nearest[facei];
const label index = near.second();
if (index == labelMax)
{
// Not found on any processor. How to map?
// Not found on any processor, or simply too far.
// How to map?
}
else if (globalCells.isLocal(index))
{
cellOrFaceLabels[facei] = globalCells.toLocal(index);
facesToSubset.set(facei);
// Mark as special if too far away
cellOrFaceLabels[facei] =
(
(near.first() < maxDistanceSqr_)
? globalCells.toLocal(index)
: -1
);
}
}
@ -295,7 +308,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
{
// With point interpolation
samplePoints_.resize(pointMap.size());
samplePoints_ = points();
sampleElements_.resize(pointMap.size(), -1);
// Store any face per point (without using pointFaces())
@ -312,34 +325,27 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
}
if (sampleSource_ == cells)
if (sampleSource_ == samplingSource::cells)
{
// samplePoints_ : per surface point a location inside the cell
// sampleElements_ : per surface point the cell
// sampleElements_ : per surface point, the cell
// samplePoints_ : per surface point, a location inside the cell
forAll(points(), pointi)
forAll(samplePoints_, pointi)
{
const point& pt = points()[pointi];
// Use _copy_ of point
const point pt = samplePoints_[pointi];
const label celli = cellOrFaceLabels[pointToFace[pointi]];
sampleElements_[pointi] = celli;
// Check if point inside cell
if
(
mesh().pointInCell
(
pt,
sampleElements_[pointi],
meshSearcher.decompMode()
)
celli >= 0
&& !mesh().pointInCell(pt, celli, meshSearcher.decompMode())
)
{
samplePoints_[pointi] = pt;
}
else
{
// Point not inside cell
// Find nearest point on faces of cell
scalar minDistSqr = VGREAT;
@ -348,7 +354,12 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
{
const face& f = mesh().faces()[facei];
pointHit info = f.nearestPoint(pt, mesh().points());
pointHit info =
f.nearestPoint
(
pt,
mesh().points()
);
if (info.distance() < minDistSqr)
{
@ -359,53 +370,56 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
}
}
}
else if (sampleSource_ == insideCells)
else if (sampleSource_ == samplingSource::insideCells)
{
// samplePoints_ : per surface point a location inside the cell
// sampleElements_ : per surface point the cell
// samplePoints_ : per surface point a location inside the cell
forAll(points(), pointi)
forAll(samplePoints_, pointi)
{
const point& pt = points()[pointi];
const label celli = cellOrFaceLabels[pointToFace[pointi]];
sampleElements_[pointi] = celli;
samplePoints_[pointi] = pt;
}
}
else
else // samplingSource::boundaryFaces
{
// samplePoints_ : per surface point a location on the boundary
// sampleElements_ : per surface point the boundary face containing
// sampleElements_ : per surface point, the boundary face containing
// the location
// samplePoints_ : per surface point, a location on the boundary
forAll(points(), pointi)
forAll(samplePoints_, pointi)
{
const point& pt = points()[pointi];
const point& pt = samplePoints_[pointi];
const label facei = cellOrFaceLabels[pointToFace[pointi]];
sampleElements_[pointi] = facei;
samplePoints_[pointi] = mesh().faces()[facei].nearestPoint
(
pt,
mesh().points()
).rawPoint();
if (facei >= 0)
{
samplePoints_[pointi] =
mesh().faces()[facei].nearestPoint
(
pt,
mesh().points()
).rawPoint();
}
}
}
}
else
{
// if sampleSource_ == cells:
// sampleElements_ : per surface triangle the cell
// sampleElements_ : per surface face, the cell
// samplePoints_ : n/a
// if sampleSource_ == insideCells:
// sampleElements_ : -1 or per surface triangle the cell
// sampleElements_ : -1 or per surface face, the cell
// samplePoints_ : n/a
// else:
// sampleElements_ : per surface triangle the boundary face
// if sampleSource_ == boundaryFaces:
// sampleElements_ : per surface face, the boundary face
// samplePoints_ : n/a
sampleElements_.transfer(cellOrFaceLabels);
samplePoints_.clear();
}
@ -431,14 +445,14 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
forAll(samplePoints_, pointi)
{
meshTools::writeOBJ(str, points()[pointi]);
vertI++;
++vertI;
meshTools::writeOBJ(str, samplePoints_[pointi]);
vertI++;
++vertI;
label elemi = sampleElements_[pointi];
meshTools::writeOBJ(str, centres[elemi]);
vertI++;
++vertI;
str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI << nl;
}
}
@ -448,11 +462,11 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
forAll(sampleElements_, triI)
{
meshTools::writeOBJ(str, faceCentres()[triI]);
vertI++;
++vertI;
label elemi = sampleElements_[triI];
meshTools::writeOBJ(str, centres[elemi]);
vertI++;
++vertI;
str << "l " << vertI-1 << ' ' << vertI << nl;
}
}
@ -474,7 +488,7 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
)
:
sampledSurface(name, mesh),
MeshStorage(),
meshedSurface(),
surfaceName_(surfaceName),
surface_
(
@ -486,7 +500,9 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
keepIds_(true),
zoneIds_(),
sampleElements_(),
samplePoints_()
samplePoints_(),
maxDistanceSqr_(Foam::sqr(GREAT)),
defaultValues_()
{}
@ -498,7 +514,7 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
)
:
sampledSurface(name, mesh, dict),
MeshStorage(),
meshedSurface(),
surfaceName_
(
meshedSurface::findFile
@ -517,8 +533,25 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
keepIds_(dict.getOrDefault("keepIds", true)),
zoneIds_(),
sampleElements_(),
samplePoints_()
samplePoints_(),
maxDistanceSqr_(Foam::sqr(GREAT)),
defaultValues_(dict.subOrEmptyDict("defaultValue"))
{
if (dict.readIfPresent("maxDistance", maxDistanceSqr_))
{
// Info<< "Limit to maxDistance " << maxDistanceSqr_ << nl;
// if (defaultValues_.empty())
// {
// Info<< "defaultValues = zero" << nl;
// }
// else
// {
// defaultValues_.writeEntry(Info);
// }
maxDistanceSqr_ = Foam::sqr(maxDistanceSqr_);
}
wordRes includePatches;
dict.readIfPresent("patches", includePatches);
includePatches.uniq();
@ -596,7 +629,7 @@ bool Foam::sampledMeshedSurface::expire()
}
sampledSurface::clearGeom();
MeshStorage::clear();
meshedSurface::clear();
zoneIds_.clear();
sampleElements_.clear();

View File

@ -35,33 +35,35 @@ Description
- 6 different modes:
- source=cells, interpolate=false:
finds per triangle centre the nearest cell centre and uses its value
- source=cells, interpolate=true
finds per triangle centre the nearest cell centre.
finds per surface face centre the \em nearest cell centre
and uses its value
- source=cells, interpolate=true:
finds per surface face centre the \em nearest cell centre.
Per surface point checks if this nearest cell is the one containing
point; otherwise projects the point onto the nearest point on
the boundary of the cell (to make sure interpolateCellPoint
gets a valid location)
- source=insideCells, interpolate=false:
finds per triangle centre the cell containing it and uses its value.
Trims triangles outside mesh.
- source=insideCells, interpolate=true
finds per surface face centre the cell containing it
and uses its value.
Trims surface faces outside of the mesh.
- source=insideCells, interpolate=true:
Per surface point interpolate cell containing it.
- source=boundaryFaces, interpolate=false:
finds per triangle centre the nearest point on the boundary
finds per surface face centre the \em nearest point on the boundary
(uncoupled faces only) and uses the value (or 0 if the nearest
is on an empty boundary)
- source=boundaryFaces, interpolate=true:
finds per triangle centre the nearest point on the boundary
finds per surface face centre the \em nearest point on the boundary
(uncoupled faces only).
Per surface point projects the point onto this boundary face
(to make sure interpolateCellPoint gets a valid location)
- since it finds a nearest per triangle each triangle is guaranteed
to be on one processor only. So after stitching (by sampledSurfaces)
the original surface should be complete.
- since it finds the nearest per surface face, each surface face
is guaranteed to be on one processor only.
So after stitching the original surface should be complete.
This is often embedded as part of a sampled surfaces function object.
@ -75,6 +77,16 @@ Usage
type meshedSurface;
surface something.obj;
source cells;
//- Max sampling distance
maxDistance 0.005;
//- Fallback for missed sampling in 'cells' mode
defaultValue
{
"p.*" 1e5;
T 273.15;
}
}
}
\endverbatim
@ -90,8 +102,11 @@ Usage
file | Alternative file name | no |
fileType | The surface format | no | (extension)
scale | Surface scaling factor | no | 0
maxDistance | Max search distance | no | GREAT
defaultValue | Value beyond max distance (dictionary) | no | empty
\endtable
SourceFiles
sampledMeshedSurface.C
sampledMeshedSurfaceTemplates.C
@ -102,8 +117,7 @@ SourceFiles
#define sampledMeshedSurface_H
#include "sampledSurface.H"
#include "MeshedSurface.H"
#include "MeshedSurfacesFwd.H"
#include "MeshedSurfaces.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -123,18 +137,19 @@ class sampledMeshedSurface
public meshedSurface
{
public:
//- Types of sampling regions
enum samplingSource
enum class samplingSource : unsigned char
{
cells,
insideCells,
boundaryFaces
cells, //!< Use nearest cell value
insideCells, //!< Surface face within a cell, or trim
boundaryFaces //!< Use nearest boundary face values
};
private:
//- Private typedefs for convenience
typedef meshedSurface MeshStorage;
//- The concrete storage type for faces and points
typedef meshedSurface Mesh;
// Private Data
@ -165,6 +180,12 @@ private:
//- Local points to sample per point
pointField samplePoints_;
//- Max search distance squared
scalar maxDistanceSqr_;
//- Out-of-range fallback values (when beyond the max search distance)
dictionary defaultValues_;
// Private Member Functions
@ -239,13 +260,13 @@ public:
//- Points of surface
virtual const pointField& points() const
{
return MeshStorage::points();
return Mesh::points();
}
//- Faces of surface
virtual const faceList& faces() const
{
return MeshStorage::surfFaces();
return Mesh::surfFaces();
}
//- Per-face zone/region information
@ -257,37 +278,37 @@ public:
//- Per-face identifier (eg, element Id)
virtual const labelList& faceIds() const
{
return MeshStorage::faceIds();
return Mesh::faceIds();
}
//- Face area vectors
virtual const vectorField& Sf() const
{
return MeshStorage::Sf();
return Mesh::Sf();
}
//- Face area magnitudes
virtual const scalarField& magSf() const
{
return MeshStorage::magSf();
return Mesh::magSf();
}
//- Face centres
virtual const vectorField& Cf() const
{
return MeshStorage::Cf();
return Mesh::Cf();
}
//- If element ids/order of the original surface are kept
virtual bool hasFaceIds() const
{
return keepIds_ && !MeshStorage::faceIds().empty();
return keepIds_ && !Mesh::faceIds().empty();
}
//- Sampling boundary values instead of cell values
bool onBoundary() const
{
return sampleSource_ == boundaryFaces;
return sampleSource_ == samplingSource::boundaryFaces;
}

View File

@ -37,8 +37,14 @@ Foam::sampledMeshedSurface::sampleOnFaces
const interpolation<Type>& sampler
) const
{
const Type deflt
(
defaultValues_.getOrDefault<Type>(sampler.psi().name(), Zero)
);
const labelList& elements = sampleElements_;
if (!onBoundary())
{
// Sample cells
@ -48,7 +54,8 @@ Foam::sampledMeshedSurface::sampleOnFaces
sampler,
elements,
faces(),
points()
points(),
deflt
);
}
@ -60,33 +67,33 @@ Foam::sampledMeshedSurface::sampleOnFaces
auto tvalues = tmp<Field<Type>>::New(elements.size());
auto& values = tvalues.ref();
const polyBoundaryMesh& pbm = mesh().boundaryMesh();
const label nBnd = mesh().nBoundaryFaces();
// Create flat boundary field
Field<Type> bVals(nBnd, Zero);
const polyBoundaryMesh& pbm = mesh().boundaryMesh();
Field<Type> bVals(mesh().nBoundaryFaces(), Zero);
const auto& bField = sampler.psi().boundaryField();
forAll(bField, patchi)
{
const label bFacei = (pbm[patchi].start() - mesh().nInternalFaces());
SubList<Type>
(
bVals,
bField[patchi].size(),
bFacei
) = bField[patchi];
SubList<Type>(bVals, pbm[patchi].range()) = bField[patchi];
}
// Sample in flat boundary field
// Sample within the flat boundary field
forAll(elements, i)
{
const label bFacei = (elements[i] - mesh().nInternalFaces());
values[i] = bVals[bFacei];
if (bFacei < 0)
{
values[i] = deflt;
}
else
{
values[i] = bVals[bFacei];
}
}
return tvalues;
@ -100,34 +107,60 @@ Foam::sampledMeshedSurface::sampleOnPoints
const interpolation<Type>& interpolator
) const
{
// One value per vertex
auto tvalues = tmp<Field<Type>>::New(sampleElements_.size());
const Type deflt
(
defaultValues_.getOrDefault<Type>(interpolator.psi().name(), Zero)
);
const labelList& elements = sampleElements_;
// One value per vertex.
// - sampleElements_ and samplePoints_ have the identical size
auto tvalues = tmp<Field<Type>>::New(elements.size());
auto& values = tvalues.ref();
if (!onBoundary())
{
// Sample cells
forAll(sampleElements_, pointi)
forAll(elements, i)
{
values[pointi] = interpolator.interpolate
(
samplePoints_[pointi],
sampleElements_[pointi]
);
const label celli = elements[i];
if (celli < 0)
{
values[i] = deflt;
}
else
{
values[i] = interpolator.interpolate
(
samplePoints_[i],
celli
);
}
}
}
else
//
// Sample boundary faces
//
forAll(elements, i)
{
// Sample boundary faces
const label facei = elements[i];
forAll(samplePoints_, pointi)
if (facei < 0)
{
const label facei = sampleElements_[pointi];
values[pointi] = interpolator.interpolate
values[i] = deflt;
}
else
{
values[i] = interpolator.interpolate
(
samplePoints_[pointi],
samplePoints_[i],
mesh().faceOwner()[facei],
facei
);

View File

@ -32,6 +32,7 @@ License
#include "polyPatch.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "uindirectPrimitivePatch.H"
#include "addToRunTimeSelectionTable.H"
@ -55,7 +56,7 @@ Foam::sampledPatch::sampledPatch
)
:
sampledSurface(name, mesh),
patchNames_(patchNames),
selectionNames_(patchNames),
triangulate_(triangulate),
needsUpdate_(true)
{}
@ -69,7 +70,7 @@ Foam::sampledPatch::sampledPatch
)
:
sampledSurface(name, mesh, dict),
patchNames_(dict.get<wordRes>("patches")),
selectionNames_(dict.get<wordRes>("patches")),
triangulate_(dict.getOrDefault("triangulate", false)),
needsUpdate_(true)
{}
@ -81,8 +82,63 @@ const Foam::labelList& Foam::sampledPatch::patchIDs() const
{
if (patchIDs_.empty())
{
patchIDs_ = mesh().boundaryMesh().patchSet(patchNames_).sortedToc();
labelList selected
(
mesh().boundaryMesh().patchSet(selectionNames_).sortedToc()
);
DynamicList<label> bad;
for (const label patchi : selected)
{
const polyPatch& pp = mesh().boundaryMesh()[patchi];
if (isA<emptyPolyPatch>(pp))
{
bad.append(patchi);
}
}
if (bad.size())
{
label nGood = (selected.size() - bad.size());
auto& os = nGood > 0 ? WarningInFunction : FatalErrorInFunction;
os << "Cannot sample an empty patch" << nl;
for (const label patchi : bad)
{
os << " "
<< mesh().boundaryMesh()[patchi].name() << nl;
}
if (nGood)
{
os << "No non-empty patches selected" << endl
<< exit(FatalError);
}
else
{
os << "Selected " << nGood << " non-empty patches" << nl;
}
patchIDs_.resize(nGood);
nGood = 0;
for (const label patchi : selected)
{
if (!bad.found(patchi))
{
patchIDs_[nGood] = patchi;
++nGood;
}
}
}
else
{
patchIDs_ = std::move(selected);
}
}
return patchIDs_;
}
@ -102,11 +158,13 @@ bool Foam::sampledPatch::expire()
}
sampledSurface::clearGeom();
MeshStorage::clear();
Mesh::clear();
patchIDs_.clear();
patchStart_.clear();
patchIndex_.clear();
patchFaceLabels_.clear();
patchStart_.clear();
needsUpdate_ = true;
return true;
@ -120,52 +178,45 @@ bool Foam::sampledPatch::update()
return false;
}
label sz = 0;
// Total number of faces selected
label numFaces = 0;
for (const label patchi : patchIDs())
{
const polyPatch& pp = mesh().boundaryMesh()[patchi];
if (isA<emptyPolyPatch>(pp))
{
FatalErrorInFunction
<< "Cannot sample an empty patch. Patch " << pp.name()
<< exit(FatalError);
}
sz += pp.size();
numFaces += pp.size();
}
// For every face (or triangle) the originating patch and local face in the
// patch.
patchIndex_.setSize(sz);
patchFaceLabels_.setSize(sz);
patchStart_.setSize(patchIDs().size());
labelList meshFaceLabels(sz);
patchStart_.resize(patchIDs().size());
sz = 0;
// The originating patch and local face in the patch.
patchIndex_.resize(numFaces);
patchFaceLabels_.resize(numFaces);
forAll(patchIDs(), i)
IndirectList<face> selectedFaces(mesh().faces(), labelList());
labelList& meshFaceIds = selectedFaces.addressing();
meshFaceIds.resize(numFaces);
numFaces = 0;
forAll(patchIDs(), idx)
{
const label patchi = patchIDs()[i];
patchStart_[i] = sz;
const label patchi = patchIDs()[idx];
const polyPatch& pp = mesh().boundaryMesh()[patchi];
const label len = pp.size();
forAll(pp, j)
{
patchIndex_[sz] = i;
patchFaceLabels_[sz] = j;
meshFaceLabels[sz] = pp.start()+j;
++sz;
}
patchStart_[idx] = numFaces;
SubList<label>(patchIndex_, len, numFaces) = idx;
SubList<label>(patchFaceLabels_, len, numFaces) = identity(len);
SubList<label>(meshFaceIds, len, numFaces) = identity(len, pp.start());
numFaces += len;
}
indirectPrimitivePatch allPatches
(
IndirectList<face>(mesh().faces(), meshFaceLabels),
mesh().points()
);
uindirectPrimitivePatch allPatches(selectedFaces, mesh().points());
this->storedPoints() = allPatches.localPoints();
this->storedFaces() = allPatches.localFaces();
@ -177,7 +228,7 @@ bool Foam::sampledPatch::update()
// too often anyhow.
if (triangulate_)
{
MeshStorage::triangulate();
Mesh::triangulate();
}
if (debug)
@ -194,10 +245,9 @@ bool Foam::sampledPatch::update()
// remap action on triangulation
void Foam::sampledPatch::remapFaces(const labelUList& faceMap)
{
// Recalculate the cells cut
if (!faceMap.empty())
{
MeshStorage::remapFaces(faceMap);
Mesh::remapFaces(faceMap);
patchFaceLabels_ = labelList
(
labelUIndList(patchFaceLabels_, faceMap)
@ -207,11 +257,11 @@ void Foam::sampledPatch::remapFaces(const labelUList& faceMap)
labelUIndList(patchIndex_, faceMap)
);
// Redo patchStart.
if (patchIndex_.size() > 0)
// Update patchStart
if (patchIndex_.size())
{
patchStart_[patchIndex_[0]] = 0;
for (label i = 1; i < patchIndex_.size(); i++)
for (label i = 1; i < patchIndex_.size(); ++i)
{
if (patchIndex_[i] != patchIndex_[i-1])
{
@ -367,7 +417,7 @@ Foam::tmp<Foam::tensorField> Foam::sampledPatch::interpolate
void Foam::sampledPatch::print(Ostream& os) const
{
os << "sampledPatch: " << name() << " :"
<< " patches:" << patchNames()
<< " patches: " << flatOutput(selectionNames_)
<< " faces:" << faces().size()
<< " points:" << points().size();
}

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2016-2019 OpenCFD Ltd.
Copyright (C) 2016-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -62,8 +62,7 @@ SourceFiles
#define sampledPatch_H
#include "sampledSurface.H"
#include "MeshedSurface.H"
#include "wordRes.H"
#include "MeshedSurfaces.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -76,19 +75,19 @@ namespace Foam
class sampledPatch
:
public MeshedSurface<face>,
public meshedSurface,
public sampledSurface
{
//- Private typedefs for convenience
typedef MeshedSurface<face> MeshStorage;
//- Mesh storage type
typedef meshedSurface Mesh;
// Private data
// Private Data
//- Name of patches
const wordRes patchNames_;
//- Selection (word/regex) of patches
const wordRes selectionNames_;
//- Corresponding patchIDs
//- The IDs for selected patches. Sorted and no emptyPolyPatch
mutable labelList patchIDs_;
//- Triangulated faces or keep faces as is
@ -97,15 +96,15 @@ class sampledPatch
//- Track if the surface needs an update
mutable bool needsUpdate_;
//- Start indices (in patchFaceLabels_) of patches
labelList patchStart_;
//- For every face (or triangle) the originating patch
labelList patchIndex_;
//- For every face (or triangle) the index in the originating patch
labelList patchFaceLabels_;
//- Start indices (in patchFaceLabels_) of patches
labelList patchStart_;
// Private Member Functions
@ -137,23 +136,28 @@ class sampledPatch
protected:
//- The selection (word/regex) of patches
const wordRes& patchNames() const
{
return patchNames_;
return selectionNames_;
}
//- The patches selected
const labelList& patchIDs() const;
//- The offset into patchIndex, patchFaceLabels
const labelList& patchStart() const
{
return patchStart_;
}
//- For each face, the patch ID
const labelList& patchIndex() const
{
return patchIndex_;
}
//- For each face, the patch local face ID
const labelList& patchFaceLabels() const
{
return patchFaceLabels_;
@ -208,13 +212,13 @@ public:
//- Points of surface
virtual const pointField& points() const
{
return MeshStorage::points();
return Mesh::points();
}
//- Faces of surface
virtual const faceList& faces() const
{
return MeshStorage::surfFaces();
return Mesh::surfFaces();
}
//- Per-face zone/region information
@ -226,19 +230,19 @@ public:
//- Face area vectors
virtual const vectorField& Sf() const
{
return MeshStorage::Sf();
return Mesh::Sf();
}
//- Face area magnitudes
virtual const scalarField& magSf() const
{
return MeshStorage::magSf();
return Mesh::magSf();
}
//- Face centres
virtual const vectorField& Cf() const
{
return MeshStorage::Cf();
return Mesh::Cf();
}

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2018 OpenCFD Ltd.
Copyright (C) 2018-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -93,20 +93,21 @@ Foam::sampledPatch::sampleOnPoints
bitSet pointDone(points().size());
forAll(faces(), cutFacei)
forAll(faces(), i)
{
const label patchi = patchIDs_[patchIndex_[cutFacei]];
const face& f = faces()[i];
const label patchi = patchIDs_[patchIndex_[i]];
const label patchFacei = patchFaceLabels_[i];
const polyPatch& pp = mesh().boundaryMesh()[patchi];
const label patchFacei = patchFaceLabels()[cutFacei];
const face& f = faces()[cutFacei];
const label facei = patchFacei + pp.start();
const label celli = own[facei];
for (const label pointi : f)
{
if (pointDone.set(pointi))
{
const label facei = patchFacei + pp.start();
const label celli = own[facei];
values[pointi] = interpolator.interpolate
(
points()[pointi],

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2016-2019 OpenCFD Ltd.
Copyright (C) 2016-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -132,14 +132,16 @@ protected:
// Protected Member Functions
//- General loop for sampling elements to faces
//- General loop for sampling elements to faces.
// The defaultValue is used for invalid (negative) elements
template<class Type>
static tmp<Field<Type>> sampleOnFaces
(
const interpolation<Type>& sampler,
const labelUList& elements,
const faceList& fcs,
const pointField& pts
const pointField& pts,
const Type& defaultValue = Type(Zero)
);

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2018-2019 OpenCFD Ltd.
Copyright (C) 2018-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -37,7 +37,8 @@ Foam::sampledSurface::sampleOnFaces
const interpolation<Type>& sampler,
const labelUList& elements,
const faceList& fcs,
const pointField& pts
const pointField& pts,
const Type& defaultValue
)
{
const label len = elements.size();
@ -57,9 +58,16 @@ Foam::sampledSurface::sampleOnFaces
for (label i=0; i < len; ++i)
{
const label celli = elements[i];
const point pt = fcs[i].centre(pts);
if (celli < 0)
{
values[i] = defaultValue;
}
else
{
const point pt = fcs[i].centre(pts);
values[i] = sampler.interpolate(pt, celli);
values[i] = sampler.interpolate(pt, celli);
}
}
return tvalues;

View File

@ -63,6 +63,7 @@ writers = writers
$(writers)/surfaceWriter.C
$(writers)/caching/surfaceWriterCaching.C
$(writers)/abaqus/abaqusSurfaceWriter.C
$(writers)/boundaryData/boundaryDataSurfaceWriter.C
$(writers)/ensight/ensightSurfaceWriter.C
$(writers)/foam/foamSurfaceWriter.C

View File

@ -154,6 +154,7 @@ bool Foam::fileFormats::ABAQUSsurfaceFormat<Face>::read
// No solids
reader.purge_solids();
reader.compact_nodes();
reader.renumber_elements_1to0(); // Convert 1-based -> 0-based
// Convert connectivity to faces

View File

@ -0,0 +1,355 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "abaqusSurfaceWriter.H"
#include "ABAQUSCore.H"
#include "IOmanip.H"
#include "ListOps.H"
#include "OSspecific.H"
#include "surfaceWriterMethods.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace surfaceWriters
{
defineTypeName(abaqusWriter);
addToRunTimeSelectionTable(surfaceWriter, abaqusWriter, word);
addToRunTimeSelectionTable(surfaceWriter, abaqusWriter, wordDict);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Field writing implementation
#include "abaqusSurfaceWriterImpl.C"
// Field writing methods
defineSurfaceWriterWriteFields(Foam::surfaceWriters::abaqusWriter);
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
// Write connectivity as CSV list
inline static void writeConnectivity
(
Ostream& os,
const label elemId,
const labelUList& elem
)
{
os << " " << elemId;
for (const label vert : elem)
{
os << ", " << (vert + 1);
}
os << nl;
}
} // End namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::surfaceWriters::abaqusWriter::writeFace
(
Ostream& os,
const labelUList& f,
const label elemId,
const label propId,
bool header
) const
{
// Only called with 3 or 4 points!
if (header)
{
os << "*ELEMENT, TYPE=S" << f.size();
if (propId >= 0)
{
os << ", ELSET=_" << propId;
}
os << nl;
}
writeConnectivity(os, elemId, f);
}
void Foam::surfaceWriters::abaqusWriter::writeGeometry
(
Ostream& os,
const meshedSurf& surf,
labelList& decompOffsets,
DynamicList<face>& decompFaces
) const
{
const pointField& points = surf.points();
const faceList& faces = surf.faces();
const labelList& zones = surf.zoneIds();
const labelList& elemIds = surf.faceIds();
// Possible to use faceIds?
bool useOrigFaceIds =
(
elemIds.size() == faces.size()
&& !ListOps::found(elemIds, lessOp1<label>(0))
);
if (useOrigFaceIds)
{
// Not possible with on-the-fly face decomposition
for (const auto& f : faces)
{
if (f.size() > 4)
{
useOrigFaceIds = false;
break;
}
}
}
os << "** Geometry" << nl;
os << nl
<< "**" << nl
<< "** Points" << nl
<< "**" << nl;
fileFormats::ABAQUSCore::writePoints(os, points, geometryScale_);
// Write faces, with on-the-fly decomposition (triangulation)
decompOffsets.resize(faces.size()+1);
decompFaces.clear();
decompOffsets[0] = 0; // The first offset is always zero
os << "**" << nl
<< "** Faces" << nl
<< "**" << nl;
// Simple tracking for change of element type/set
labelPair prevOutput(-1, -1);
label elemId = 0; // The element-id
forAll(faces, facei)
{
const face& f = faces[facei];
if (useOrigFaceIds)
{
// When available and not decomposed
elemId = elemIds[facei];
}
// 1-offset for PID
const label propId = 1 + (facei < zones.size() ? zones[facei] : 0);
const label n = f.size();
bool header =
(prevOutput.first() != n || prevOutput.second() != propId);
if (header)
{
// Update values
prevOutput.first() = n;
prevOutput.second() = propId;
}
if (n == 3 || n == 4)
{
writeFace(os, f, ++elemId, propId, header);
}
else
{
// Decompose into tris
prevOutput.first() = 3;
f.triangles(points, decompFaces);
for
(
label decompi = decompOffsets[facei];
decompi < decompFaces.size();
++decompi
)
{
writeFace
(
os,
decompFaces[decompi],
++elemId,
propId,
header
);
header = false;
}
}
// The end offset, which is the next begin offset
decompOffsets[facei+1] = decompFaces.size();
}
os << "**" << nl
<< "**" << nl;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::surfaceWriters::abaqusWriter::abaqusWriter()
:
surfaceWriter(),
geometryScale_(1),
fieldScale_(),
noGeometry_(false),
outputLayout_(outputLayoutType::BY_FIELD)
{}
Foam::surfaceWriters::abaqusWriter::abaqusWriter
(
const dictionary& options
)
:
surfaceWriter(options),
geometryScale_(options.getOrDefault<scalar>("scale", 1)),
fieldScale_(options.subOrEmptyDict("fieldScale")),
noGeometry_(options.getOrDefault("noGeometry", false)),
outputLayout_(outputLayoutType::BY_FIELD)
{}
Foam::surfaceWriters::abaqusWriter::abaqusWriter
(
const meshedSurf& surf,
const fileName& outputPath,
bool parallel,
const dictionary& options
)
:
abaqusWriter(options)
{
open(surf, outputPath, parallel);
}
Foam::surfaceWriters::abaqusWriter::abaqusWriter
(
const pointField& points,
const faceList& faces,
const fileName& outputPath,
bool parallel,
const dictionary& options
)
:
abaqusWriter(options)
{
open(points, faces, outputPath, parallel);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::fileName Foam::surfaceWriters::abaqusWriter::write()
{
checkOpen();
// Geometry:
// 1) rootdir/<TIME>/surfaceName.abq
// 2) rootdir/geometry/surfaceName_<TIME>.abq
fileName outputFile;
switch (outputLayout_)
{
case outputLayoutType::BY_TIME:
{
outputFile = outputPath_;
if (useTimeDir() && !timeName().empty())
{
// Splice in time-directory
outputFile =
outputPath_.path() / timeName() / outputPath_.name();
}
break;
}
case outputLayoutType::BY_FIELD:
{
outputFile = outputPath_ / "geometry" / outputPath_.name();
if (!timeName().empty())
{
// Append time information to file name
outputFile += '_' + timeName();
}
break;
}
}
outputFile.ext("abq");
if (verbose_)
{
Info<< "Writing abaqus geometry to " << outputFile << endl;
}
const meshedSurf& surf = surface();
if (Pstream::master() || !parallel_)
{
if (!isDir(outputFile.path()))
{
mkDir(outputFile.path());
}
OFstream os(outputFile);
labelList decompOffsets;
DynamicList<face> decompFaces;
writeGeometry(os, surf, decompOffsets, decompFaces);
}
wroteGeom_ = true;
return outputFile;
}
// ************************************************************************* //

View File

@ -0,0 +1,239 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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::surfaceWriters::abaqusWriter
Description
A surface writer for the ABAQUS file format - both surface mesh and fields
The formatOptions for abaqus:
\table
Property | Description | Required | Default
scale | output geometry scaling | no | 1
fieldScale | output field scaling (dictionary) | no | empty
noGeometry | Suppress geometry output (beta feature) | no | false
\endtable
For example,
\verbatim
formatOptions
{
abaqus
{
scale 1000; // [m] -> [mm]
fieldScale
{
"p.*" 0.01; // [Pa] -> [mbar]
}
}
}
\endverbatim
\heading Output file locations
The \c rootdir normally corresponds to something like
\c postProcessing/\<name\>
\subheading Geometry
\verbatim
rootdir
`-- surfaceName0
`-- geometry
`-- surfaceName0_<time>.abq
\endverbatim
\subheading Fields
\verbatim
rootdir
`-- surfaceName0
`-- field0
| `-- surfaceName0_<time>.{inp}
`-- field1
`-- surfaceName0_<time>.{inp}
\endverbatim
Note
Output variable scaling does not apply to integer types such as Ids.
SourceFiles
abaqusSurfaceWriter.C
abaqusSurfaceWriterImpl.C
\*---------------------------------------------------------------------------*/
#ifndef abaqusSurfaceWriter_H
#define abaqusSurfaceWriter_H
#include "surfaceWriter.H"
#include "ABAQUSCore.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace surfaceWriters
{
/*---------------------------------------------------------------------------*\
Class abaqusWriter Declaration
\*---------------------------------------------------------------------------*/
class abaqusWriter
:
public surfaceWriter
{
// Private Data
//- Output geometry scaling
const scalar geometryScale_;
//- Output field scaling
const dictionary fieldScale_;
//- BETA feature
bool noGeometry_;
//- Output directory layouts
enum class outputLayoutType
{
BY_TIME = 0,
BY_FIELD,
};
//- Output directory layout
outputLayoutType outputLayout_;
// Private Member Functions
//- Write face element (tri/quad)
void writeFace
(
Ostream& os,
const labelUList& f,
const label elemId, //!< 1-based Element Id
const label propId, //!< 1-based Property Id
const bool header = true
) const;
//- Write the surface mesh geometry, tracking face decomposition
//
// \param decompOffsets begin/end offsets (size+1) into decompFaces
// \param decompFaces Non tri/quad decomposed into triangles
void writeGeometry
(
Ostream& os,
const meshedSurf& surf,
labelList& decompOffsets,
DynamicList<face>& decompFaces
) const;
//- Write a face-based value
template<class Type>
Ostream& writeFaceValue
(
Ostream& os,
const Type& value,
const label elemId //!< 0-based Element Id
) const;
//- Templated write operation
template<class Type>
fileName writeTemplate
(
const word& fieldName, //!< Name of field
const Field<Type>& localValues //!< Local field values to write
);
public:
//- Declare type-name, virtual type (with debug switch)
TypeNameNoDebug("abaqus");
// Constructors
//- Default construct
abaqusWriter();
//- Construct with some output options
explicit abaqusWriter(const dictionary& options);
//- Construct from components
abaqusWriter
(
const meshedSurf& surf,
const fileName& outputPath,
bool parallel = Pstream::parRun(),
const dictionary& options = dictionary()
);
//- Construct from components
abaqusWriter
(
const pointField& points,
const faceList& faces,
const fileName& outputPath,
bool parallel = Pstream::parRun(),
const dictionary& options = dictionary()
);
//- Destructor
virtual ~abaqusWriter() = default;
// Member Functions
//- Format uses faceIds as part of its output
virtual bool usesFaceIds() const // override
{
return true;
}
//- Write surface geometry to file.
virtual fileName write(); // override
declareSurfaceWriterWriteMethod(label);
declareSurfaceWriterWriteMethod(scalar);
declareSurfaceWriterWriteMethod(vector);
declareSurfaceWriterWriteMethod(sphericalTensor);
declareSurfaceWriterWriteMethod(symmTensor);
declareSurfaceWriterWriteMethod(tensor);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace surfaceWriters
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,308 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "OFstream.H"
#include "IOmanip.H"
#include "OSspecific.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type>
Foam::Ostream& Foam::surfaceWriters::abaqusWriter::writeFaceValue
(
Ostream& os,
const Type& value,
const label elemId
) const
{
if (fileFormats::ABAQUSCore::isEncodedSolidId(elemId))
{
const label solidId =
fileFormats::ABAQUSCore::decodeSolidElementId(elemId);
const label sideNum =
fileFormats::ABAQUSCore::decodeSolidSideNum(elemId);
// Element-id: 0-based to 1-based
// Side-id: always 1-based
os << (solidId + 1) << ", P" << sideNum;
}
else
{
// Element-id: 0-based to 1-based
os << (elemId + 1) << ", P";
}
os << ", ";
if (pTraits<Type>::nComponents == 1)
{
os << value;
}
else
{
os << mag(value);
}
os << nl;
return os;
}
template<class Type>
Foam::fileName Foam::surfaceWriters::abaqusWriter::writeTemplate
(
const word& fieldName,
const Field<Type>& localValues
)
{
checkOpen();
// Field:
// 1) rootdir/<TIME>/<field>_surfaceName.inp
// 2) rootdir/<field>/surfaceName_<TIME>.inp
fileName outputFile;
switch (outputLayout_)
{
case outputLayoutType::BY_TIME:
{
outputFile = outputPath_;
if (useTimeDir() && !timeName().empty())
{
// Splice in time-directory
outputFile =
outputPath_.path() / timeName() / outputPath_.name();
}
// Append <field>_surfaceName
outputFile /= fieldName + '_' + outputPath_.name();
break;
}
case outputLayoutType::BY_FIELD:
{
outputFile = outputPath_ / fieldName / outputPath_.name();
if (!timeName().empty())
{
// Append time information to file name
outputFile += '_' + timeName();
}
break;
}
}
outputFile.ext("inp");
// Output scaling for the variable, but not for integer types.
// could also solve with clever templating
const scalar varScale =
(
std::is_integral<Type>::value
? scalar(1)
: fieldScale_.getOrDefault<scalar>(fieldName, 1)
);
if (verbose_)
{
Info<< "Writing field " << fieldName;
if (!equal(varScale, 1))
{
Info<< " (scaling " << varScale << ')';
}
Info<< " to " << outputFile << endl;
}
// Implicit geometry merge()
tmp<Field<Type>> tfield = mergeField(localValues) * varScale;
const meshedSurf& surf = surface();
if (Pstream::master() || !parallel_)
{
const auto& values = tfield();
if (!isDir(outputFile.path()))
{
mkDir(outputFile.path());
}
// const scalar timeValue(0);
// Additional bookkeeping for decomposing non tri/quad
labelList decompOffsets;
DynamicList<face> decompFaces;
OFstream os(outputFile);
if (noGeometry_ || wroteGeom_)
{
// Geometry already written (or suppressed)
// - still need decomposition information
fileFormats::ABAQUSCore::faceDecomposition
(
surf.points(),
surf.faces(),
decompOffsets,
decompFaces
);
}
else
{
// Write geometry (separate file)
OFstream osGeom(outputFile.lessExt().ext("abq"));
writeGeometry(osGeom, surf, decompOffsets, decompFaces);
}
// Write field
// *DLOAD
// Element-based: // elemId, P, 1000
// Surface-based: // elemId, P4, 1000
os << "**" << nl
<< "** field = " << fieldName << nl
<< "** type = " << pTraits<Type>::typeName << nl;
if (useTimeDir() && !timeName().empty())
{
os << "** time = " << timeName() << nl;
}
os << "**" << nl
<< "*DLOAD" << nl;
// Regular (undecomposed) faces
const faceList& faces = surf.faces();
const labelList& elemIds = surf.faceIds();
// Possible to use faceIds?
const bool useOrigFaceIds =
(
elemIds.size() == faces.size()
&& decompFaces.empty()
);
label elemId = 0;
if (this->isPointData())
{
forAll(faces, facei)
{
if (useOrigFaceIds)
{
// When available and not decomposed
elemId = elemIds[facei];
}
const label beginElemId = elemId;
// Any face decomposition
for
(
label decompi = decompOffsets[facei];
decompi < decompOffsets[facei+1];
++decompi
)
{
const face& f = decompFaces[decompi];
Type v = Zero;
for (const label verti : f)
{
v += values[verti];
}
v /= f.size();
writeFaceValue(os, v, elemId); // 0-based
++elemId;
}
// Face not decomposed
if (beginElemId == elemId)
{
const face& f = faces[facei];
Type v = Zero;
for (const label verti : f)
{
v += values[verti];
}
v /= f.size();
writeFaceValue(os, v, elemId); // 0-based
++elemId;
}
}
}
else
{
auto valIter = values.cbegin();
forAll(faces, facei)
{
if (useOrigFaceIds)
{
// When available and not decomposed
elemId = elemIds[facei];
}
const Type v(*valIter);
++valIter;
label nValues =
max
(
label(1),
(decompOffsets[facei+1] - decompOffsets[facei])
);
while (nValues--)
{
writeFaceValue(os, v, elemId); // 0-based
++elemId;
}
}
}
os << "**" << nl
<< "**" << nl;
}
wroteGeom_ = true;
return outputFile;
}
// ************************************************************************* //

View File

@ -121,6 +121,44 @@ sampled
}
// Sample and write for meshed surfaces
sampleMesh
{
type surfaces;
libs (sampling);
log true;
writeControl writeTime;
writeInterval 1;
surfaceFormat vtk;
fields (p rho U T);
surfaces
{
// Oversized sampling - for general testing
meshed_sample
{
type meshedSurface;
surface oversized_sample.obj;
source cells;
maxDistance 0.0025;
defaultValue
{
T 273;
}
}
meshed_interpolate
{
$meshed_sample;
interpolate true;
}
}
}
// * * * * * * * * * * * * * * * Calculations * * * * * * * * * * * * * * * //
massflow
@ -144,6 +182,7 @@ areaAverage
operation weightedAreaAverage;
weightField rhoU;
weightFields ( rho U none ); // 2012 and later
fields ( p pTotal );
}
@ -156,6 +195,7 @@ areaIntegrate
operation weightedAreaIntegrate;
weightField rhoU;
weightFields ( rho U ); // 2012 and later
fields ( T );
}

View File

@ -16,8 +16,9 @@ runApplication subsetMesh c0 -patch walls -overwrite
runApplication splitMeshRegions -cellZones -overwrite
# create register face and cell zones
rm log.topoSet
runApplication topoSet -region cabin -dict system/topoSetDictRegister
rm -f log.topoSet.register
runApplication -s register \
topoSet -region cabin -dict system/topoSetDictRegister
# set the initial fields
restore0Dir

View File

@ -16,8 +16,9 @@ runApplication subsetMesh c0 -patch walls -overwrite
runApplication splitMeshRegions -cellZones -overwrite
# create register face and cell zones
rm log.topoSet
runApplication topoSet -region cabin -dict system/topoSetDictRegister
rm -f log.topoSet.register
runApplication -o -s register \
topoSet -region cabin -dict system/topoSetDictRegister
# set the initial fields
restore0Dir

View File

@ -83,6 +83,27 @@ functions
( 60 "<system>/solverControls.60")
);
}
// Verify values for sampling multiple faceZones
inletFaces
{
type surfaces;
libs (sampling);
writeControl onEnd;
region cabin;
surfaceFormat vtk;
fields (H2O U);
surfaces
{
inletFaces
{
type faceZone;
zones (inletFaces f1Zone);
interpolate false;
}
}
}
}

View File

@ -44,6 +44,28 @@ actions
type cellSet;
action invert;
}
{
name f0tmp;
type faceSet;
action new;
source patchToFace;
patch inlet;
}
{
name inletFaces;
type faceZoneSet;
action new;
source setToFaceZone;
faceSet f0tmp;
}
{
name f0tmp;
type faceSet;
action remove;
}
);

View File

@ -3,15 +3,12 @@ cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
# create mesh
runApplication blockMesh
restore0Dir
# initialise with potentialFoam solution
runApplication blockMesh
runApplication potentialFoam
# run the solver
runApplication $(getApplication)
#------------------------------------------------------------------------------

View File

@ -0,0 +1,16 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
restore0Dir
runApplication blockMesh
runApplication potentialFoam
runApplication decomposePar
runParallel $(getApplication)
#------------------------------------------------------------------------------

View File

@ -53,7 +53,27 @@ maxDeltaT 1e-03;
functions
{
surfaceFieldValue1
avgInlets
{
type surfaceFieldValue;
libs (fieldFunctionObjects);
enabled yes;
writeControl writeTime;
log yes;
writeFields off;
surfaceFormat vtk;
regionType patch;
// name inlet;
names (inlet "inlet.*");
operation weightedAverage;
weightField phi;
fields
(
T
);
}
avgOutlets
{
type surfaceFieldValue;
libs (fieldFunctionObjects);

View File

@ -0,0 +1,30 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object decomposeParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
numberOfSubdomains 4;
method scotch;
coeffs
{
n (2 2 1);
}
distributed no;
roots ( );
// ************************************************************************* //

View File

@ -48,7 +48,7 @@ runTimeModifiable yes;
functions
{
surfaceFieldValue1
avgOutlets
{
type surfaceFieldValue;
libs (fieldFunctionObjects);