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

This commit is contained in:
Sergio Ferraris
2013-06-07 10:11:37 +01:00
520 changed files with 102879 additions and 3007 deletions

View File

@ -56,6 +56,7 @@ primitives/Vector/lists/vectorListIOList.C
primitives/Tensor2D/tensor2D/tensor2D.C
primitives/SphericalTensor2D/sphericalTensor2D/sphericalTensor2D.C
primitives/SymmTensor2D/symmTensor2D/symmTensor2D.C
primitives/Vector2D/vector2D/vector2D.C
primitives/complex/complex.C
@ -554,6 +555,7 @@ $(Fields)/sphericalTensorField/sphericalTensorField.C
$(Fields)/diagTensorField/diagTensorField.C
$(Fields)/symmTensorField/symmTensorField.C
$(Fields)/tensorField/tensorField.C
$(Fields)/quaternionField/quaternionField.C
$(Fields)/triadField/triadField.C
$(Fields)/complexFields/complexFields.C
@ -573,6 +575,7 @@ $(Fields)/symmTensorField/symmTensorIOField.C
$(Fields)/symmTensorField/symmTensorFieldIOField.C
$(Fields)/tensorField/tensorIOField.C
$(Fields)/tensorField/tensorFieldIOField.C
$(Fields)/quaternionField/quaternionIOField.C
$(Fields)/triadField/triadIOField.C
$(Fields)/transformField/transformField.C

View File

@ -2677,8 +2677,8 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
}
template <class Type>
template <class FindNearestOp>
template<class Type>
template<class FindNearestOp>
Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
(
const point& sample,
@ -2728,8 +2728,8 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
}
template <class Type>
template <class FindNearestOp>
template<class Type>
template<class FindNearestOp>
Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
(
const linePointRef& ln,
@ -2799,8 +2799,8 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLineAny
// Find nearest intersection
template <class Type>
template <class FindIntersectOp>
template<class Type>
template<class FindIntersectOp>
Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
(
const point& start,
@ -2813,8 +2813,8 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
// Find nearest intersection
template <class Type>
template <class FindIntersectOp>
template<class Type>
template<class FindIntersectOp>
Foam::pointIndexHit Foam::indexedOctree<Type>::findLineAny
(
const point& start,

View File

@ -203,7 +203,7 @@ private:
// Query
//- Find nearest point to line.
template <class FindNearestOp>
template<class FindNearestOp>
void findNearest
(
const label nodeI,
@ -288,7 +288,7 @@ private:
// intersection point.
// findAny=true : return any intersection
// findAny=false: return nearest (to start) intersection
template <class FindIntersectOp>
template<class FindIntersectOp>
void traverseNode
(
const bool findAny,
@ -307,7 +307,7 @@ private:
) const;
//- Find any or nearest intersection
template <class FindIntersectOp>
template<class FindIntersectOp>
pointIndexHit findLine
(
const bool findAny,
@ -327,7 +327,7 @@ private:
// ) const;
//- Find any or nearest intersection of line between start and end.
template <class FindIntersectOp>
template<class FindIntersectOp>
pointIndexHit findLine
(
const bool findAny,
@ -540,7 +540,7 @@ public:
// - bool : any point found nearer than nearestDistSqr
// - label: index in shapes
// - point: actual nearest point found
template <class FindNearestOp>
template<class FindNearestOp>
pointIndexHit findNearest
(
const point& sample,
@ -563,7 +563,7 @@ public:
// ) const;
//- Low level: calculate nearest starting from subnode.
template <class FindNearestOp>
template<class FindNearestOp>
void findNearest
(
const label nodeI,
@ -590,7 +590,7 @@ public:
point& linePoint
) const;
template <class FindNearestOp>
template<class FindNearestOp>
pointIndexHit findNearest
(
const linePointRef& ln,
@ -615,7 +615,7 @@ public:
) const;
//- Find nearest intersection of line between start and end.
template <class FindIntersectOp>
template<class FindIntersectOp>
pointIndexHit findLine
(
const point& start,
@ -624,7 +624,7 @@ public:
) const;
//- Find any intersection of line between start and end.
template <class FindIntersectOp>
template<class FindIntersectOp>
pointIndexHit findLineAny
(
const point& start,

View File

@ -618,11 +618,11 @@ bool Foam::Time::writeObject
)
);
timeDict.add("value", value());
timeDict.add("value", timeToUserTime(value()));
timeDict.add("name", string(tmName));
timeDict.add("index", timeIndex_);
timeDict.add("deltaT", deltaT_);
timeDict.add("deltaT0", deltaT0_);
timeDict.add("deltaT", timeToUserTime(deltaT_));
timeDict.add("deltaT0", timeToUserTime(deltaT0_));
timeDict.regIOobject::writeObject(fmt, ver, cmp);
bool writeOK = objectRegistry::writeObject(fmt, ver, cmp);

View File

@ -0,0 +1,45 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "quaternionField.H"
#define TEMPLATE
#include "FieldFunctionsM.C"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "undefFieldFunctionsM.H"
// ************************************************************************* //

View File

@ -0,0 +1,64 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Typedef
Foam::quaternionField
Description
Specialisation of Field\<T\> for quaternion.
SourceFiles
quaternionField.C
\*---------------------------------------------------------------------------*/
#ifndef quaternionField_H
#define quaternionField_H
#include "Field.H"
#include "quaternion.H"
#define TEMPLATE
#include "FieldFunctionsM.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
typedef Field<quaternion> quaternionField;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "undefFieldFunctionsM.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,43 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
quaternionField with IO.
\*---------------------------------------------------------------------------*/
#include "quaternionIOField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTemplateTypeNameAndDebugWithName
(
quaternionIOField,
"quaternionField",
0
);
}
// ************************************************************************* //

View File

@ -0,0 +1,49 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Typedef
Foam::quaternionIOField
Description
quaternionField with IO.
\*---------------------------------------------------------------------------*/
#ifndef quaternionIOField_H
#define quaternionIOField_H
#include "quaternionField.H"
#include "IOField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
typedef IOField<quaternion> quaternionIOField;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -49,21 +49,16 @@ readField
<< endl;
}
// Patch or patch-groups. (using non-wild card entries of dictionaries)
// 1. Handle explicit patch names
forAllConstIter(dictionary, dict, iter)
{
if (iter().isDict() && !iter().keyword().isPattern())
{
const labelList patchIDs = bmesh_.findIndices
(
iter().keyword(),
true
);
label patchi = bmesh_.findPatchID(iter().keyword());
forAll(patchIDs, i)
if (patchi != -1)
{
label patchi = patchIDs[i];
this->set
(
patchi,
@ -78,7 +73,43 @@ readField
}
}
// Check for wildcard patch overrides
// 2. Patch-groups. (using non-wild card entries of dictionaries)
// (patchnames already matched above)
// Note: in order of entries in the dictionary (first patchGroups wins)
forAllConstIter(dictionary, dict, iter)
{
if (iter().isDict() && !iter().keyword().isPattern())
{
const labelList patchIDs = bmesh_.findIndices
(
iter().keyword(),
true // use patchGroups
);
forAll(patchIDs, i)
{
label patchi = patchIDs[i];
if (!this->set(patchi))
{
this->set
(
patchi,
PatchField<Type>::New
(
bmesh_[patchi],
field,
iter().dict()
)
);
}
}
}
}
// 3. Wildcard patch overrides
forAll(bmesh_, patchi)
{
if (!this->set(patchi))

View File

@ -165,7 +165,7 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
<< setw(8) << setprecision(4) << totFaceCellRatio/totNprocs
<< setw(8) << setprecision(4) << maxFaceCellRatio
<< " "
<< setw(8) << totNInt/totNprocs
<< setw(8) << scalar(totNInt)/totNprocs
<< setw(8) << maxNInt
<< " "
<< setw(8) << setprecision(4) << totRatio/totNprocs
@ -535,4 +535,123 @@ const Foam::labelListListList& Foam::GAMGAgglomeration::boundaryFaceMap
}
bool Foam::GAMGAgglomeration::checkRestriction
(
labelList& newRestrict,
label& nNewCoarse,
const lduAddressing& fineAddressing,
const labelUList& restrict,
const label nCoarse
)
{
if (fineAddressing.size() != restrict.size())
{
FatalErrorIn
(
"checkRestriction(..)"
) << "nCells:" << fineAddressing.size()
<< " agglom:" << restrict.size()
<< abort(FatalError);
}
// Seed (master) for every region
labelList master(identity(fineAddressing.size()));
// Now loop and transport master through region
const labelUList& lower = fineAddressing.lowerAddr();
const labelUList& upper = fineAddressing.upperAddr();
while (true)
{
label nChanged = 0;
forAll(lower, faceI)
{
label own = lower[faceI];
label nei = upper[faceI];
if (restrict[own] == restrict[nei])
{
// coarse-mesh-internal face
if (master[own] < master[nei])
{
master[nei] = master[own];
nChanged++;
}
else if (master[own] > master[nei])
{
master[own] = master[nei];
nChanged++;
}
}
}
reduce(nChanged, sumOp<label>());
if (nChanged == 0)
{
break;
}
}
// Count number of regions/masters per coarse cell
labelListList coarseToMasters(nCoarse);
nNewCoarse = 0;
forAll(restrict, cellI)
{
labelList& masters = coarseToMasters[restrict[cellI]];
if (findIndex(masters, master[cellI]) == -1)
{
masters.append(master[cellI]);
nNewCoarse++;
}
}
if (nNewCoarse > nCoarse)
{
//WarningIn("GAMGAgglomeration::checkRestriction(..)")
// << "Have " << nCoarse
// << " agglomerated cells but " << nNewCoarse
// << " disconnected regions" << endl;
// Keep coarseToMasters[0] the original coarse, allocate new ones
// for the others
labelListList coarseToNewCoarse(coarseToMasters.size());
nNewCoarse = nCoarse;
forAll(coarseToMasters, coarseI)
{
const labelList& masters = coarseToMasters[coarseI];
labelList& newCoarse = coarseToNewCoarse[coarseI];
newCoarse.setSize(masters.size());
newCoarse[0] = coarseI;
for (label i = 1; i < newCoarse.size(); i++)
{
newCoarse[i] = nNewCoarse++;
}
}
newRestrict.setSize(fineAddressing.size());
forAll(restrict, cellI)
{
label coarseI = restrict[cellI];
label index = findIndex(coarseToMasters[coarseI], master[cellI]);
newRestrict[cellI] = coarseToNewCoarse[coarseI][index];
}
return false;
}
else
{
return true;
}
}
// ************************************************************************* //

View File

@ -474,6 +474,16 @@ public:
const labelListListList& boundaryFaceMap(const label fineLeveli)
const;
//- Given restriction determines if coarse cells are connected.
// Return ok is so, otherwise creates new restriction that is
static bool checkRestriction
(
labelList& newRestrict,
label& nNewCoarse,
const lduAddressing& fineAddressing,
const labelUList& restrict,
const label nCoarse
);
};

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -58,6 +58,12 @@ Foam::pointBoundaryMesh::pointBoundaryMesh
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::pointBoundaryMesh::findPatchID(const word& patchName) const
{
return mesh()().boundaryMesh().findPatchID(patchName);
}
Foam::labelList Foam::pointBoundaryMesh::findIndices
(
const keyType& key,

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -96,6 +96,9 @@ public:
return mesh_;
}
//- Find patch index given a name
label findPatchID(const word& patchName) const;
//- Find patch indices given a name
labelList findIndices(const keyType&, const bool useGroups) const;

View File

@ -97,7 +97,6 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const faceList& fcs = mesh.faces();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
@ -154,31 +153,16 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
{
label faceI = pp.start() + i;
vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]];
skew[faceI] = primitiveMeshTools::boundaryFaceSkewness
(
mesh,
p,
fCtrs,
fAreas,
vector normal = fAreas[faceI];
normal /= mag(normal) + VSMALL;
vector d = normal*(normal & Cpf);
// Skewness vector
vector sv =
Cpf
- ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
vector svHat = sv/(mag(sv) + VSMALL);
// Normalisation distance calculated as the approximate distance
// from the face centre to the edge of the face in the direction
// of the skewness
scalar fd = 0.4*mag(d) + VSMALL;
const face& f = fcs[faceI];
forAll(f, pi)
{
fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
}
// Normalised skewness
skew[faceI] = mag(sv)/fd;
faceI,
cellCtrs[own[faceI]]
);
}
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -136,7 +136,7 @@ Foam::pointField Foam::coupledPolyPatch::getAnchorPoints
{
pointField anchors(faces.size());
if (transform == COINCIDENTFULLMATCH)
if (transform != COINCIDENTFULLMATCH)
{
// Return the first point
forAll(faces, faceI)
@ -344,7 +344,15 @@ void Foam::coupledPolyPatch::calcTransformTensors
Pout<< " error:" << error << endl;
}
if
if (transform == NOORDERING)
{
forwardT_.setSize(0);
reverseT_.setSize(0);
separation_.setSize(0);
collocated_ = boolList(1, true);
}
else if
(
transform == ROTATIONAL
|| (

View File

@ -193,7 +193,7 @@ void Foam::processorPolyPatch::calcGeometry(PstreamBuffers& pBufs)
// For small face area calculation the results of the area
// calculation have been found to only be accurate to ~1e-20
if (magSf < SMALL && nbrMagSf < SMALL)
if (magSf < SMALL || nbrMagSf < SMALL)
{
// Undetermined normal. Use dummy normal to force separation
// check.

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -276,6 +276,17 @@ public:
const polyBoundaryMesh& bm
);
//- Return a pointer to a new patch created on freestore from
// dictionary
static autoPtr<polyPatch> New
(
const word& patchType,
const word& name,
const dictionary& dict,
const label index,
const polyBoundaryMesh& bm
);
//- Destructor
virtual ~polyPatch();

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -95,6 +95,26 @@ Foam::autoPtr<Foam::polyPatch> Foam::polyPatch::New
word patchType(dict.lookup("type"));
dict.readIfPresent("geometricType", patchType);
return polyPatch::New(patchType, name, dict, index, bm);
}
Foam::autoPtr<Foam::polyPatch> Foam::polyPatch::New
(
const word& patchType,
const word& name,
const dictionary& dict,
const label index,
const polyBoundaryMesh& bm
)
{
if (debug)
{
Info<< "polyPatch::New(const word&, const word&, const dictionary&, "
"const label, const polyBoundaryMesh&) : constructing polyPatch"
<< endl;
}
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(patchType);

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -30,6 +30,7 @@ License
#include "PatchToolsGatherAndMerge.C"
#include "PatchToolsSearch.C"
#include "PatchToolsSortEdges.C"
#include "PatchToolsSortPoints.C"
#include "PatchToolsNormals.C"
#include "PatchToolsMatch.C"

View File

@ -170,6 +170,18 @@ public:
const PrimitivePatch<Face, FaceList, PointField, PointType>&
);
//- Return point-edge addressing sorted by order around the point.
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
static labelListList sortedPointEdges
(
const PrimitivePatch<Face, FaceList, PointField, PointType>&
);
//- If 2 face neighbours: label of face where ordering of edge
// is consistent with righthand walk.

View File

@ -0,0 +1,118 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "PatchTools.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
Foam::labelListList
Foam::PatchTools::sortedPointEdges
(
const PrimitivePatch<Face, FaceList, PointField, PointType>& p
)
{
// Now order the edges of each point according to whether they share a
// face
const labelListList& pointEdges = p.pointEdges();
const edgeList& edges = p.edges();
const labelListList& edgeFaces = p.edgeFaces();
const labelListList& faceEdges = p.faceEdges();
// create the lists for the various results. (resized on completion)
labelListList sortedPointEdges(pointEdges.size());
DynamicList<label> newEdgeList;
forAll(pointEdges, pointI)
{
const labelList& pEdges = pointEdges[pointI];
label edgeI = pEdges[0];
label prevFaceI = edgeFaces[edgeI][0];
newEdgeList.clear();
newEdgeList.setCapacity(pEdges.size());
do
{
newEdgeList.append(edgeI);
// Cross edge to next face
const labelList& eFaces = edgeFaces[edgeI];
if (eFaces.size() != 2)
{
break;
}
label faceI = eFaces[0];
if (faceI == prevFaceI)
{
faceI = eFaces[1];
}
// Cross face to next edge
const labelList& fEdges = faceEdges[faceI];
forAll(fEdges, feI)
{
const label nextEdgeI = fEdges[feI];
const edge& nextEdge = edges[nextEdgeI];
if
(
nextEdgeI != edgeI
&& (nextEdge.start() == pointI || nextEdge.end() == pointI)
)
{
edgeI = nextEdgeI;
break;
}
}
prevFaceI = faceI;
} while (edgeI != pEdges[0]);
if (newEdgeList.size() == pEdges.size())
{
sortedPointEdges[pointI] = newEdgeList;
}
}
return sortedPointEdges;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -75,76 +75,6 @@ calcPointEdges() const
<< "calcPointEdges() finished calculating pointEdges"
<< endl;
}
// Now order the edges of each point according to whether they share a
// face
// DynamicList<label> newEdgeList;
// forAll(pe, pointI)
// {
// const labelList& pEdges = pe[pointI];
// label edgeI = pEdges[0];
// label prevFaceI = edgeFaces()[edgeI][0];
// newEdgeList.clear();
// newEdgeList.setCapacity(pEdges.size());
// do
// {
// newEdgeList.append(edgeI);
// // Cross edge to next face
// const labelList& eFaces = edgeFaces()[edgeI];
// if (eFaces.size() != 2)
// {
// break;
// }
// label faceI = eFaces[0];
// if (faceI == prevFaceI)
// {
// faceI = eFaces[1];
// }
// // Cross face to next edge
// const labelList& fEdges = faceEdges()[faceI];
// forAll(fEdges, feI)
// {
// const label nextEdgeI = fEdges[feI];
// const edge& nextEdge = edges()[nextEdgeI];
// if
// (
// nextEdgeI != edgeI
// && (nextEdge.start() == pointI || nextEdge.end() == pointI)
// )
// {
// edgeI = nextEdgeI;
// break;
// }
// }
// prevFaceI = faceI;
// } while (edgeI != pEdges[0]);
// if (newEdgeList.size() == pEdges.size())
// {
// pe[pointI] = newEdgeList;
// }
// }
// if (debug)
// {
// Info<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
// << "calcPointEdges() finished ordering pointEdges"
// << endl;
// }
}

View File

@ -843,8 +843,9 @@ bool Foam::primitiveMesh::checkFaceFlatness
{
if (nSummed > 0)
{
Info<< " Face flatness (1 = flat, 0 = butterfly) : average = "
<< sumFlatness / nSummed << " min = " << minFlatness << endl;
Info<< " Face flatness (1 = flat, 0 = butterfly) : min = "
<< minFlatness << " average = " << sumFlatness / nSummed
<< endl;
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -64,6 +64,44 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness
return mag(sv)/fd;
}
Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label faceI,
const point& ownCc
)
{
vector Cpf = fCtrs[faceI] - ownCc;
vector normal = fAreas[faceI];
normal /= mag(normal) + VSMALL;
vector d = normal*(normal & Cpf);
// Skewness vector
vector sv =
Cpf
- ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
vector svHat = sv/(mag(sv) + VSMALL);
// Normalisation distance calculated as the approximate distance
// from the face centre to the edge of the face in the direction
// of the skewness
scalar fd = 0.4*mag(d) + VSMALL;
const face& f = mesh.faces()[faceI];
forAll(f, pi)
{
fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
}
// Normalised skewness
return mag(sv)/fd;
}
Foam::scalar Foam::primitiveMeshTools::faceOrthogonality
(
@ -119,7 +157,6 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const faceList& fcs = mesh.faces();
tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
scalarField& skew = tskew();
@ -145,31 +182,15 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
{
vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]];
vector normal = fAreas[faceI];
normal /= mag(normal) + VSMALL;
vector d = normal*(normal & Cpf);
// Skewness vector
vector sv =
Cpf
- ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
vector svHat = sv/(mag(sv) + VSMALL);
// Normalisation distance calculated as the approximate distance
// from the face centre to the edge of the face in the direction
// of the skewness
scalar fd = 0.4*mag(d) + VSMALL;
const face& f = fcs[faceI];
forAll(f, pi)
{
fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
}
// Normalised skewness
skew[faceI] = mag(sv)/fd;
skew[faceI] = boundaryFaceSkewness
(
mesh,
p,
fCtrs,
fAreas,
faceI,
cellCtrs[own[faceI]]
);
}
return tskew;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -48,32 +48,6 @@ namespace Foam
class primitiveMeshTools
{
protected:
// Protected Member Functions
//- Skewness of single face
static scalar faceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label faceI,
const point& ownCc,
const point& neiCc
);
//- Orthogonality of single face
static scalar faceOrthogonality
(
const point& ownCc,
const point& neiCc,
const vector& s
);
public:
//- Generate non-orthogonality field (internal faces only)
@ -154,6 +128,41 @@ public:
const PackedBoolList& internalOrCoupledFace
);
// Helpers: single face check
//- Skewness of single face
static scalar faceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label faceI,
const point& ownCc,
const point& neiCc
);
//- Skewness of single boundary face
static scalar boundaryFaceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label faceI,
const point& ownCc
);
//- Orthogonality of single face
static scalar faceOrthogonality
(
const point& ownCc,
const point& neiCc,
const vector& s
);
};

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -439,6 +439,14 @@ Foam::point Foam::plane::planePlaneIntersect
}
Foam::plane::side Foam::plane::sideOfPlane(const point& p) const
{
const scalar angle((p - basePoint_) & unitVector_);
return (angle < 0 ? FLIP : NORMAL);
}
void Foam::plane::writeDict(Ostream& os) const
{
os.writeKeyword("planeType") << "pointAndNormal"

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -69,6 +69,12 @@ public:
second() = s;
}
//- Construct from FixedList
inline Pair(const FixedList<Type, 2>& fl)
:
FixedList<Type, 2>(fl)
{}
//- Construct from Istream
inline Pair(Istream& is)
:

View File

@ -0,0 +1,147 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::SymmTensor2D
Description
Templated 2D symmetric tensor derived from VectorSpace adding construction
from 4 components, element access using xx(), xy() etc. member functions
and the inner-product (dot-product) and outer-product of two Vectors
(tensor-product) operators.
SourceFiles
SymmTensor2DI.H
\*---------------------------------------------------------------------------*/
#ifndef SymmTensor2D_H
#define SymmTensor2D_H
#include "VectorSpace.H"
#include "SphericalTensor2D.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class SymmTensor2D Declaration
\*---------------------------------------------------------------------------*/
template<class Cmpt>
class SymmTensor2D
:
public VectorSpace<SymmTensor2D<Cmpt>, Cmpt, 3>
{
public:
//- Equivalent type of labels used for valid component indexing
typedef SymmTensor2D<label> labelType;
// Member constants
enum
{
rank = 2 // Rank of SymmTensor2D is 2
};
// Static data members
static const char* const typeName;
static const char* componentNames[];
static const SymmTensor2D zero;
static const SymmTensor2D one;
static const SymmTensor2D max;
static const SymmTensor2D min;
static const SymmTensor2D I;
//- Component labeling enumeration
enum components { XX, XY, YY };
// Constructors
//- Construct null
inline SymmTensor2D();
//- Construct given VectorSpace
inline SymmTensor2D(const VectorSpace<SymmTensor2D<Cmpt>, Cmpt, 3>&);
//- Construct given SphericalTensor
inline SymmTensor2D(const SphericalTensor2D<Cmpt>&);
//- Construct given the three components
inline SymmTensor2D
(
const Cmpt txx, const Cmpt txy,
const Cmpt tyy
);
//- Construct from Istream
SymmTensor2D(Istream&);
// Member Functions
// Access
inline const Cmpt& xx() const;
inline const Cmpt& xy() const;
inline const Cmpt& yy() const;
inline Cmpt& xx();
inline Cmpt& xy();
inline Cmpt& yy();
//- Transpose
inline const SymmTensor2D<Cmpt>& T() const;
// Member Operators
//- Construct given SphericalTensor2D
inline void operator=(const SphericalTensor2D<Cmpt>&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Include inline implementations
#include "SymmTensor2DI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,503 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "Vector2D.H"
#include "Tensor2D.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Cmpt>
inline SymmTensor2D<Cmpt>::SymmTensor2D()
{}
template<class Cmpt>
inline SymmTensor2D<Cmpt>::SymmTensor2D
(
const VectorSpace<SymmTensor2D<Cmpt>, Cmpt, 3>& vs
)
:
VectorSpace<SymmTensor2D<Cmpt>, Cmpt, 3>(vs)
{}
template<class Cmpt>
inline SymmTensor2D<Cmpt>::SymmTensor2D(const SphericalTensor2D<Cmpt>& st)
{
this->v_[XX] = st.ii(); this->v_[XY] = 0;
this->v_[YY] = st.ii();
}
template<class Cmpt>
inline SymmTensor2D<Cmpt>::SymmTensor2D
(
const Cmpt txx, const Cmpt txy,
const Cmpt tyy
)
{
this->v_[XX] = txx; this->v_[XY] = txy;
this->v_[YY] = tyy;
}
template<class Cmpt>
inline SymmTensor2D<Cmpt>::SymmTensor2D(Istream& is)
:
VectorSpace<SymmTensor2D<Cmpt>, Cmpt, 3>(is)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Cmpt>
inline const Cmpt& SymmTensor2D<Cmpt>::xx() const
{
return this->v_[XX];
}
template<class Cmpt>
inline const Cmpt& SymmTensor2D<Cmpt>::xy() const
{
return this->v_[XY];
}
template<class Cmpt>
inline const Cmpt& SymmTensor2D<Cmpt>::yy() const
{
return this->v_[YY];
}
template<class Cmpt>
inline Cmpt& SymmTensor2D<Cmpt>::xx()
{
return this->v_[XX];
}
template<class Cmpt>
inline Cmpt& SymmTensor2D<Cmpt>::xy()
{
return this->v_[XY];
}
template<class Cmpt>
inline Cmpt& SymmTensor2D<Cmpt>::yy()
{
return this->v_[YY];
}
template<class Cmpt>
inline const SymmTensor2D<Cmpt>& SymmTensor2D<Cmpt>::T() const
{
return *this;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Cmpt>
inline void SymmTensor2D<Cmpt>::operator=(const SphericalTensor2D<Cmpt>& st)
{
this->v_[XX] = st.ii(); this->v_[XY] = 0;
this->v_[YY] = st.ii();
}
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
//- Inner-product between two symmetric tensors
template<class Cmpt>
inline Tensor2D<Cmpt>
operator&(const SymmTensor2D<Cmpt>& st1, const SymmTensor2D<Cmpt>& st2)
{
return Tensor2D<Cmpt>
(
st1.xx()*st2.xx() + st1.xy()*st2.xy(),
st1.xx()*st2.xy() + st1.xy()*st2.yy(),
st1.xy()*st2.xx() + st1.yy()*st2.xy(),
st1.xy()*st2.xy() + st1.yy()*st2.yy()
);
}
//- Double-dot-product between a symmetric tensor and a symmetric tensor
template<class Cmpt>
inline Cmpt
operator&&(const SymmTensor2D<Cmpt>& st1, const SymmTensor2D<Cmpt>& st2)
{
return
(
st1.xx()*st2.xx() + 2*st1.xy()*st2.xy()
+ st1.yy()*st2.yy()
);
}
//- Inner-product between a symmetric tensor and a vector
template<class Cmpt>
inline Vector2D<Cmpt>
operator&(const SymmTensor2D<Cmpt>& st, const Vector2D<Cmpt>& v)
{
return Vector2D<Cmpt>
(
st.xx()*v.x() + st.xy()*v.y(),
st.xy()*v.x() + st.yy()*v.y()
);
}
//- Inner-product between a vector and a symmetric tensor
template<class Cmpt>
inline Vector2D<Cmpt>
operator&(const Vector2D<Cmpt>& v, const SymmTensor2D<Cmpt>& st)
{
return Vector2D<Cmpt>
(
v.x()*st.xx() + v.y()*st.xy(),
v.x()*st.xy() + v.y()*st.yy()
);
}
template<class Cmpt>
inline Cmpt magSqr(const SymmTensor2D<Cmpt>& st)
{
return
(
magSqr(st.xx()) + 2*magSqr(st.xy())
+ magSqr(st.yy())
);
}
//- Return the trace of a symmetric tensor
template<class Cmpt>
inline Cmpt tr(const SymmTensor2D<Cmpt>& st)
{
return st.xx() + st.yy();
}
//- Return the spherical part of a symmetric tensor
template<class Cmpt>
inline SphericalTensor2D<Cmpt> sph(const SymmTensor2D<Cmpt>& st)
{
return (1.0/2.0)*tr(st);
}
//- Return the symmetric part of a symmetric tensor, i.e. itself
template<class Cmpt>
inline const SymmTensor2D<Cmpt>& symm(const SymmTensor2D<Cmpt>& st)
{
return st;
}
//- Return twice the symmetric part of a symmetric tensor
template<class Cmpt>
inline SymmTensor2D<Cmpt> twoSymm(const SymmTensor2D<Cmpt>& st)
{
return 2*st;
}
//- Return the deviatoric part of a symmetric tensor
template<class Cmpt>
inline SymmTensor2D<Cmpt> dev(const SymmTensor2D<Cmpt>& st)
{
return st - SphericalTensor2D<Cmpt>::oneThirdI*tr(st);
}
//- Return the deviatoric part of a symmetric tensor
template<class Cmpt>
inline SymmTensor2D<Cmpt> dev2(const SymmTensor2D<Cmpt>& st)
{
return st - SphericalTensor2D<Cmpt>::twoThirdsI*tr(st);
}
//- Return the determinant of a symmetric tensor
template<class Cmpt>
inline Cmpt det(const SymmTensor2D<Cmpt>& st)
{
return
(
st.xx()*st.yy() - st.xy()*st.xy()
);
}
//- Return the cofactor symmetric tensor of a symmetric tensor
template<class Cmpt>
inline SymmTensor2D<Cmpt> cof(const SymmTensor2D<Cmpt>& st)
{
return SymmTensor2D<Cmpt>
(
st.yy(), -st.xy(),
st.xx()
);
}
//- Return the inverse of a symmetric tensor give the determinant
template<class Cmpt>
inline SymmTensor2D<Cmpt> inv(const SymmTensor2D<Cmpt>& st, const Cmpt detst)
{
return cof(st)/detst;
}
//- Return the inverse of a symmetric tensor
template<class Cmpt>
inline SymmTensor2D<Cmpt> inv(const SymmTensor2D<Cmpt>& st)
{
return inv(st, det(st));
}
//- Return the 1st invariant of a symmetric tensor
template<class Cmpt>
inline Cmpt invariantI(const SymmTensor2D<Cmpt>& st)
{
return tr(st);
}
//- Return the 2nd invariant of a symmetric tensor
template<class Cmpt>
inline Cmpt invariantII(const SymmTensor2D<Cmpt>& st)
{
return
(
0.5*sqr(tr(st))
- 0.5*
(
st.xx()*st.xx() + st.xy()*st.xy()
+ st.xy()*st.xy() + st.yy()*st.yy()
)
);
}
//- Return the 3rd invariant of a symmetric tensor
template<class Cmpt>
inline Cmpt invariantIII(const SymmTensor2D<Cmpt>& st)
{
return det(st);
}
template<class Cmpt>
inline SymmTensor2D<Cmpt>
operator+(const SphericalTensor2D<Cmpt>& spt1, const SymmTensor2D<Cmpt>& st2)
{
return SymmTensor2D<Cmpt>
(
spt1.ii() + st2.xx(), st2.xy(),
spt1.ii() + st2.yy()
);
}
template<class Cmpt>
inline SymmTensor2D<Cmpt>
operator+(const SymmTensor2D<Cmpt>& st1, const SphericalTensor2D<Cmpt>& spt2)
{
return SymmTensor2D<Cmpt>
(
st1.xx() + spt2.ii(), st1.xy(),
st1.yy() + spt2.ii()
);
}
template<class Cmpt>
inline SymmTensor2D<Cmpt>
operator-(const SphericalTensor2D<Cmpt>& spt1, const SymmTensor2D<Cmpt>& st2)
{
return SymmTensor2D<Cmpt>
(
spt1.ii() - st2.xx(), -st2.xy(),
spt1.ii() - st2.yy()
);
}
template<class Cmpt>
inline SymmTensor2D<Cmpt>
operator-(const SymmTensor2D<Cmpt>& st1, const SphericalTensor2D<Cmpt>& spt2)
{
return SymmTensor2D<Cmpt>
(
st1.xx() - spt2.ii(), st1.xy(),
st1.yy() - spt2.ii()
);
}
//- Inner-product between a spherical symmetric tensor and a symmetric tensor
template<class Cmpt>
inline SymmTensor2D<Cmpt>
operator&(const SphericalTensor2D<Cmpt>& spt1, const SymmTensor2D<Cmpt>& st2)
{
return SymmTensor2D<Cmpt>
(
spt1.ii()*st2.xx(), spt1.ii()*st2.xy(),
spt1.ii()*st2.yy()
);
}
//- Inner-product between a tensor and a spherical tensor
template<class Cmpt>
inline SymmTensor2D<Cmpt>
operator&(const SymmTensor2D<Cmpt>& st1, const SphericalTensor2D<Cmpt>& spt2)
{
return SymmTensor2D<Cmpt>
(
st1.xx()*spt2.ii(), st1.xy()*spt2.ii(),
st1.yy()*spt2.ii()
);
}
//- Double-dot-product between a spherical tensor and a symmetric tensor
template<class Cmpt>
inline Cmpt
operator&&(const SphericalTensor2D<Cmpt>& spt1, const SymmTensor2D<Cmpt>& st2)
{
return(spt1.ii()*st2.xx() + spt1.ii()*st2.yy());
}
//- Double-dot-product between a tensor and a spherical tensor
template<class Cmpt>
inline Cmpt
operator&&(const SymmTensor2D<Cmpt>& st1, const SphericalTensor2D<Cmpt>& spt2)
{
return(st1.xx()*spt2.ii() + st1.yy()*spt2.ii());
}
template<class Cmpt>
inline SymmTensor2D<Cmpt> sqr(const Vector2D<Cmpt>& v)
{
return SymmTensor2D<Cmpt>
(
v.x()*v.x(), v.x()*v.y(),
v.y()*v.y()
);
}
template<class Cmpt>
class outerProduct<SymmTensor2D<Cmpt>, Cmpt>
{
public:
typedef SymmTensor2D<Cmpt> type;
};
template<class Cmpt>
class outerProduct<Cmpt, SymmTensor2D<Cmpt> >
{
public:
typedef SymmTensor2D<Cmpt> type;
};
template<class Cmpt>
class innerProduct<SymmTensor2D<Cmpt>, SymmTensor2D<Cmpt> >
{
public:
typedef SymmTensor2D<Cmpt> type;
};
template<class Cmpt>
class innerProduct<SymmTensor2D<Cmpt>, Vector2D<Cmpt> >
{
public:
typedef Vector2D<Cmpt> type;
};
template<class Cmpt>
class innerProduct<Vector2D<Cmpt>, SymmTensor2D<Cmpt> >
{
public:
typedef Vector2D<Cmpt> type;
};
template<class Cmpt>
class typeOfSum<SphericalTensor2D<Cmpt>, SymmTensor2D<Cmpt> >
{
public:
typedef SymmTensor2D<Cmpt> type;
};
template<class Cmpt>
class typeOfSum<SymmTensor2D<Cmpt>, SphericalTensor2D<Cmpt> >
{
public:
typedef SymmTensor2D<Cmpt> type;
};
template<class Cmpt>
class innerProduct<SphericalTensor2D<Cmpt>, SymmTensor2D<Cmpt> >
{
public:
typedef SymmTensor2D<Cmpt> type;
};
template<class Cmpt>
class innerProduct<SymmTensor2D<Cmpt>, SphericalTensor2D<Cmpt> >
{
public:
typedef SymmTensor2D<Cmpt> type;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,85 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "symmTensor2D.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template<>
const char* const symmTensor2D::typeName = "symmTensor2D";
template<>
const char* symmTensor2D::componentNames[] =
{
"xx", "xy",
"yy"
};
template<>
const symmTensor2D symmTensor2D::zero
(
0, 0,
0
);
template<>
const symmTensor2D symmTensor2D::one
(
1, 1,
1
);
template<>
const symmTensor2D symmTensor2D::max
(
VGREAT, VGREAT,
VGREAT
);
template<>
const symmTensor2D symmTensor2D::min
(
-VGREAT, -VGREAT,
-VGREAT
);
template<>
const symmTensor2D symmTensor2D::I
(
1, 0,
1
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,63 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Typedef
Foam::symmTensor2D
Description
SymmTensor2D of scalars.
SourceFiles
symmTensor2D.C
\*---------------------------------------------------------------------------*/
#ifndef symmTensor2D_H
#define symmTensor2D_H
#include "SymmTensor2D.H"
#include "contiguous.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
typedef SymmTensor2D<scalar> symmTensor2D;
//- Data associated with symmTensor2D type are contiguous
template<>
inline bool contiguous<symmTensor2D>() {return true;}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -46,6 +46,9 @@ SourceFiles
namespace Foam
{
template<class Cmpt>
class SymmTensor2D;
/*---------------------------------------------------------------------------*\
Class Tensor2D Declaration
\*---------------------------------------------------------------------------*/
@ -90,6 +93,9 @@ public:
//- Construct given VectorSpace
inline Tensor2D(const VectorSpace<Tensor2D<Cmpt>, Cmpt, 4>&);
//- Construct given SymmTensor2D
inline Tensor2D(const SymmTensor2D<Cmpt>&);
//- Construct given SphericalTensor2D
inline Tensor2D(const SphericalTensor2D<Cmpt>&);
@ -136,6 +142,9 @@ public:
// Member Operators
//- Copy SymmTensor2D
inline void operator=(const SymmTensor2D<Cmpt>&);
//- Copy SphericalTensor2D
inline void operator=(const SphericalTensor2D<Cmpt>&);
};

View File

@ -42,6 +42,14 @@ inline Tensor2D<Cmpt>::Tensor2D(const VectorSpace<Tensor2D<Cmpt>, Cmpt, 4>& vs)
{}
template<class Cmpt>
inline Tensor2D<Cmpt>::Tensor2D(const SymmTensor2D<Cmpt>& st)
{
this->v_[XX] = st.xx(); this->v_[XY] = st.xy();
this->v_[YX] = st.xy(); this->v_[YY] = st.yy();
}
template<class Cmpt>
inline Tensor2D<Cmpt>::Tensor2D(const SphericalTensor2D<Cmpt>& st)
{
@ -159,6 +167,14 @@ inline Tensor2D<Cmpt> Tensor2D<Cmpt>::T() const
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Cmpt>
inline void Tensor2D<Cmpt>::operator=(const SymmTensor2D<Cmpt>& st)
{
this->v_[XX] = st.xx(); this->v_[XY] = st.xy();
this->v_[YX] = st.xy(); this->v_[YY] = st.yy();
}
template<class Cmpt>
inline void Tensor2D<Cmpt>::operator=(const SphericalTensor2D<Cmpt>& st)
{
@ -240,23 +256,23 @@ inline SphericalTensor2D<Cmpt> sph(const Tensor2D<Cmpt>& t)
//- Return the symmetric part of a tensor
template<class Cmpt>
inline Tensor2D<Cmpt> symm(const Tensor2D<Cmpt>& t)
inline SymmTensor2D<Cmpt> symm(const Tensor2D<Cmpt>& t)
{
return Tensor2D<Cmpt>
return SymmTensor2D<Cmpt>
(
t.xx(), 0.5*(t.xy() + t.yx()),
0.5*(t.yx() + t.xy()), t.yy()
t.yy()
);
}
//- Return the twice the symmetric part of a tensor
template<class Cmpt>
inline Tensor2D<Cmpt> twoSymm(const Tensor2D<Cmpt>& t)
inline SymmTensor2D<Cmpt> twoSymm(const Tensor2D<Cmpt>& t)
{
return Tensor2D<Cmpt>
return SymmTensor2D<Cmpt>
(
t.xx() + t.xx(), t.xy() + t.yx(),
t.yx() + t.xy(), t.yy() + t.yy()
t.yy() + t.yy()
);
}
@ -356,7 +372,7 @@ inline Cmpt invariantIII(const Tensor2D<Cmpt>& t)
}
// * * * * * * * * * Mixed Tensor SphericalTensor Operators * * * * * * * * //
template<class Cmpt>
inline Tensor2D<Cmpt>
@ -455,6 +471,114 @@ operator&&(const Tensor2D<Cmpt>& t1, const SphericalTensor2D<Cmpt>& st2)
}
// * * * * * * * * * * Mixed Tensor SymmTensor Operators * * * * * * * * * * //
template<class Cmpt>
inline Tensor2D<Cmpt>
operator+(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
{
return Tensor2D<Cmpt>
(
st1.xx() + t2.xx(), st1.xy() + t2.xy(),
st1.xy() + t2.yx(), st1.yy() + t2.yy()
);
}
template<class Cmpt>
inline Tensor2D<Cmpt>
operator+(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
{
return Tensor2D<Cmpt>
(
t1.xx() + st2.xx(), t1.xy() + st2.xy(),
t1.yx() + st2.xy(), t1.yy() + st2.yy()
);
}
template<class Cmpt>
inline Tensor2D<Cmpt>
operator-(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
{
return Tensor2D<Cmpt>
(
st1.xx() - t2.xx(), st1.xy() - t2.xy(),
st1.xy() - t2.yx(), st1.yy() - t2.yy()
);
}
template<class Cmpt>
inline Tensor2D<Cmpt>
operator-(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
{
return Tensor2D<Cmpt>
(
t1.xx() - st2.xx(), t1.xy() - st2.xy(),
t1.yx() - st2.xy(), t1.yy() - st2.yy()
);
}
//- Inner-product between a spherical tensor and a tensor
template<class Cmpt>
inline Tensor2D<Cmpt>
operator&(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
{
return Tensor2D<Cmpt>
(
st1.xx()*t2.xx() + st1.xy()*t2.yx(),
st1.xx()*t2.xy() + st1.xy()*t2.yy(),
st1.xy()*t2.xx() + st1.yy()*t2.yx(),
st1.xy()*t2.xy() + st1.yy()*t2.yy()
);
}
//- Inner-product between a tensor and a spherical tensor
template<class Cmpt>
inline Tensor2D<Cmpt>
operator&(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
{
return Tensor2D<Cmpt>
(
t1.xx()*st2.xx() + t1.xy()*st2.xy(),
t1.xx()*st2.xy() + t1.xy()*st2.yy(),
t1.yx()*st2.xx() + t1.yy()*st2.xy(),
t1.yx()*st2.xy() + t1.yy()*st2.yy()
);
}
//- Double-dot-product between a spherical tensor and a tensor
template<class Cmpt>
inline Cmpt
operator&&(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
{
return
(
st1.xx()*t2.xx() + st1.xy()*t2.xy()
+ st1.xy()*t2.yx() + st1.yy()*t2.yy()
);
}
//- Double-dot-product between a tensor and a spherical tensor
template<class Cmpt>
inline Cmpt
operator&&(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
{
return
(
t1.xx()*st2.xx() + t1.xy()*st2.xy()
+ t1.yx()*st2.xy() + t1.yy()*st2.yy()
);
}
template<class Cmpt>
class typeOfSum<SphericalTensor2D<Cmpt>, Tensor2D<Cmpt> >
{

View File

@ -58,21 +58,67 @@ Foam::quaternion Foam::slerp
const scalar t
)
{
// Calculate angle between the quaternions
scalar cosHalfTheta = qa & qb;
label sign = 1;
if (mag(cosHalfTheta) >= 1)
if ((qa & qb) < 0)
{
return qa;
sign = -1;
}
scalar halfTheta = acos(cosHalfTheta);
scalar sinHalfTheta = sqrt(1.0 - sqr(cosHalfTheta));
return qa*pow((inv(qa)*sign*qb), t);
}
scalar wa = sin((1 - t)*halfTheta)/sinHalfTheta;
scalar wb = sin(t*halfTheta)/sinHalfTheta;
return wa*qa + wb*qb;
Foam::quaternion Foam::exp(const quaternion& q)
{
const scalar magV = mag(q.v());
if (magV == 0)
{
return quaternion(1, vector::zero);
}
const scalar expW = exp(q.w());
return quaternion
(
expW*cos(magV),
expW*sin(magV)*q.v()/magV
);
}
Foam::quaternion Foam::pow(const quaternion& q, const label power)
{
const scalar magQ = mag(q);
const scalar magV = mag(q.v());
quaternion powq(q.v());
if (magV != 0 && magQ != 0)
{
powq /= magV;
powq *= power*acos(q.w()/magQ);
}
return pow(magQ, power)*exp(powq);
}
Foam::quaternion Foam::pow(const quaternion& q, const scalar power)
{
const scalar magQ = mag(q);
const scalar magV = mag(q.v());
quaternion powq(q.v());
if (magV != 0 && magQ != 0)
{
powq /= magV;
powq *= power*acos(q.w()/magQ);
}
return pow(magQ, power)*exp(powq);
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -149,6 +149,10 @@ public:
//- The rotation tensor corresponding the quaternion
inline tensor R() const;
//- Return a vector of euler angles (rotations in radians about
// the x, y and z axes.
inline vector eulerAngles(const quaternion& q) const;
inline quaternion normalized() const;
@ -209,7 +213,7 @@ inline scalar mag(const quaternion& q);
//- Return the conjugate of the given quaternion
inline quaternion conjugate(const quaternion& q);
//- Return the normailzed (unit) quaternion of the given quaternion
//- Return the normalized (unit) quaternion of the given quaternion
inline quaternion normalize(const quaternion& q);
//- Return the inverse of the given quaternion
@ -226,6 +230,15 @@ quaternion slerp
const scalar t
);
//- Exponent of a quaternion
quaternion exp(const quaternion& q);
//- Power of a quaternion
quaternion pow(const quaternion& q, const label power);
//- Power of a quaternion
quaternion pow(const quaternion& q, const scalar power);
//- Data associated with quaternion type are contiguous
template<>
inline bool contiguous<quaternion>() {return true;}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -232,6 +232,58 @@ inline Foam::quaternion Foam::quaternion::invTransform
}
inline Foam::tensor Foam::quaternion::R() const
{
scalar w2 = sqr(w());
scalar x2 = sqr(v().x());
scalar y2 = sqr(v().y());
scalar z2 = sqr(v().z());
scalar txy = 2*v().x()*v().y();
scalar twz = 2*w()*v().z();
scalar txz = 2*v().x()*v().z();
scalar twy = 2*w()*v().y();
scalar tyz = 2*v().y()*v().z();
scalar twx = 2*w()*v().x();
return tensor
(
w2 + x2 - y2 - z2, txy - twz, txz + twy,
txy + twz, w2 - x2 + y2 - z2, tyz - twx,
txz - twy, tyz + twx, w2 - x2 - y2 + z2
);
}
inline Foam::vector Foam::quaternion::eulerAngles(const quaternion& q) const
{
vector angles(vector::zero);
const scalar& w = q.w();
const vector& v = q.v();
angles[0] = Foam::atan2
(
2*(w*v.x() + v.y()*v.z()),
1 - 2*(sqr(v.x()) + sqr(v.y()))
);
angles[1] = Foam::asin(2*(w*v.y() - v.z()*v.x()));
angles[2] = Foam::atan2
(
2*(w*v.z() + v.x()*v.y()),
1 - 2*(sqr(v.y()) + sqr(v.z()))
);
// Convert to degrees
// forAll(angles, aI)
// {
// angles[aI] = Foam::radToDeg(angles[aI]);
// }
return angles;
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void Foam::quaternion::operator=(const quaternion& q)
@ -323,29 +375,6 @@ inline Foam::quaternion Foam::normalize(const quaternion& q)
}
inline Foam::tensor Foam::quaternion::R() const
{
scalar w2 = sqr(w());
scalar x2 = sqr(v().x());
scalar y2 = sqr(v().y());
scalar z2 = sqr(v().z());
scalar txy = 2*v().x()*v().y();
scalar twz = 2*w()*v().z();
scalar txz = 2*v().x()*v().z();
scalar twy = 2*w()*v().y();
scalar tyz = 2*v().y()*v().z();
scalar twx = 2*w()*v().x();
return tensor
(
w2 + x2 - y2 - z2, txy - twz, txz + twy,
txy + twz, w2 - x2 + y2 - z2, tyz - twx,
txz - twy, tyz + twx, w2 - x2 - y2 + z2
);
}
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
inline bool Foam::operator==(const quaternion& q1, const quaternion& q2)

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,6 +25,7 @@ License
#include "cachedRandom.H"
#include "OSspecific.H"
#include "PstreamReduceOps.H"
// * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * //
@ -142,6 +143,78 @@ Foam::scalar Foam::cachedRandom::position
}
template<>
Foam::label Foam::cachedRandom::globalSample01()
{
scalar value = -GREAT;
if (Pstream::master())
{
value = scalar01();
}
reduce(value, maxOp<scalar>());
return round(value);
}
template<>
Foam::scalar Foam::cachedRandom::globalSample01()
{
scalar value = -GREAT;
if (Pstream::master())
{
value = scalar01();
}
reduce(value, maxOp<scalar>());
return value;
}
template<>
Foam::label Foam::cachedRandom::globalPosition
(
const label& start,
const label& end
)
{
label value = labelMin;
if (Pstream::master())
{
value = round(scalar01()*(end - start));
}
reduce(value, maxOp<label>());
return start + value;
}
template<>
Foam::scalar Foam::cachedRandom::globalPosition
(
const scalar& start,
const scalar& end
)
{
scalar value = -GREAT;
if (Pstream::master())
{
value = scalar01()*(end - start);
}
reduce(value, maxOp<scalar>());
return start + value;
}
void Foam::cachedRandom::operator=(const cachedRandom& cr)
{
seed_ = cr.seed_;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -121,17 +121,33 @@ public:
// Evaluation
//- Return a sample whose components lie in the range 0-1
template<class Type>
Type sample01();
// Random numbers
//- Return a sample whose components lie in the range 0-1
template<class Type>
Type sample01();
//- Return a sample between start and end
template<class Type>
Type position(const Type& start, const Type& end);
//- Return a sample between start and end
template<class Type>
Type position(const Type& start, const Type& end);
//- Randomise value in the range 0-1
template<class Type>
void randomise01(Type& value);
//- Randomise value in the range 0-1
template<class Type>
void randomise01(Type& value);
// Global random numbers - consistent across all processors
//- Return a sample whose components lie in the range 0-1
template<class Type>
Type globalSample01();
//- Return a sample between start and end
template<class Type>
Type globalPosition(const Type& start, const Type& end);
//- Randomise value in the range 0-1
template<class Type>
void globalRandomise01(Type& value);
// Operators
@ -160,6 +176,22 @@ scalar cachedRandom::position<scalar>
const scalar& end
);
template<>
label cachedRandom::globalSample01<label>();
template<>
scalar cachedRandom::globalSample01<scalar>();
template<>
label cachedRandom::globalPosition<label>(const label& start, const label& end);
template<>
scalar cachedRandom::globalPosition<scalar>
(
const scalar& start,
const scalar& end
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "cachedRandom.H"
#include "Pstream.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -61,4 +62,50 @@ void Foam::cachedRandom::randomise01(Type& value)
}
template<class Type>
Type Foam::cachedRandom::globalSample01()
{
Type value = -GREAT*pTraits<Type>::one;
if (Pstream::master())
{
value = sample01<Type>();
}
reduce(value, maxOp<Type>());
return value;
}
template<class Type>
Type Foam::cachedRandom::globalPosition(const Type& start, const Type& end)
{
Type value = -GREAT*pTraits<Type>::one;
if (Pstream::master())
{
value = position<Type>(start, end);
}
reduce(value, maxOp<Type>());
return value;
}
template<class Type>
void Foam::cachedRandom::globalRandomise01(Type& value)
{
value = -GREAT*pTraits<Type>::one;
if (Pstream::master())
{
value = sample01<Type>();
}
reduce(value, maxOp<Type>());
}
// ************************************************************************* //

View File

@ -124,49 +124,49 @@ void Foam::triad::orthogonalize()
{
for (int i=0; i<2; i++)
{
scalar o01 = mag(operator[](0) & operator[](1));
scalar o02 = mag(operator[](0) & operator[](2));
scalar o12 = mag(operator[](1) & operator[](2));
scalar o01 = mag(operator[](0) & operator[](1));
scalar o02 = mag(operator[](0) & operator[](2));
scalar o12 = mag(operator[](1) & operator[](2));
if (o01 < o02 && o01 < o12)
{
operator[](2) = orthogonal(operator[](0), operator[](1));
if (o01 < o02 && o01 < o12)
{
operator[](2) = orthogonal(operator[](0), operator[](1));
// if (o02 < o12)
// {
// operator[](1) = orthogonal(operator[](0), operator[](2));
// }
// else
// {
// operator[](0) = orthogonal(operator[](1), operator[](2));
// }
}
else if (o02 < o12)
{
operator[](1) = orthogonal(operator[](0), operator[](2));
// if (o02 < o12)
// {
// operator[](1) = orthogonal(operator[](0), operator[](2));
// }
// else
// {
// operator[](0) = orthogonal(operator[](1), operator[](2));
// }
}
else if (o02 < o12)
{
operator[](1) = orthogonal(operator[](0), operator[](2));
// if (o01 < o12)
// {
// operator[](2) = orthogonal(operator[](0), operator[](1));
// }
// else
// {
// operator[](0) = orthogonal(operator[](1), operator[](2));
// }
}
else
{
operator[](0) = orthogonal(operator[](1), operator[](2));
// if (o01 < o12)
// {
// operator[](2) = orthogonal(operator[](0), operator[](1));
// }
// else
// {
// operator[](0) = orthogonal(operator[](1), operator[](2));
// }
}
else
{
operator[](0) = orthogonal(operator[](1), operator[](2));
// if (o02 < o01)
// {
// operator[](1) = orthogonal(operator[](0), operator[](2));
// }
// else
// {
// operator[](2) = orthogonal(operator[](0), operator[](1));
// }
}
// if (o02 < o01)
// {
// operator[](1) = orthogonal(operator[](0), operator[](2));
// }
// else
// {
// operator[](2) = orthogonal(operator[](0), operator[](1));
// }
}
}
}
}
@ -174,19 +174,19 @@ void Foam::triad::orthogonalize()
void Foam::triad::operator+=(const triad& t2)
{
if (t2.set(0) && !set(0))
{
operator[](0) = t2.operator[](0);
}
bool preset[3];
if (t2.set(1) && !set(1))
for (direction i=0; i<3; i++)
{
operator[](1) = t2.operator[](1);
}
if (t2.set(2) && !set(2))
{
operator[](2) = t2.operator[](2);
if (t2.set(i) && !set(i))
{
operator[](i) = t2.operator[](i);
preset[i] = true;
}
else
{
preset[i] = false;
}
}
if (set() && t2.set())
@ -196,6 +196,12 @@ void Foam::triad::operator+=(const triad& t2)
for (direction i=0; i<3; i++)
{
if (preset[i])
{
signd[i] = 0;
continue;
}
scalar mostAligned = -1;
for (direction j=0; j<3; j++)
{
@ -222,10 +228,7 @@ void Foam::triad::operator+=(const triad& t2)
}
}
}
}
for (direction i=0; i<3; i++)
{
operator[](i) += signd[i]*t2.operator[](correspondance[i]);
}
}
@ -373,4 +376,33 @@ void Foam::triad::operator=(const tensor& t)
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::scalar Foam::diff(const triad& A, const triad& B)
{
triad tmpA = A.sortxyz();
triad tmpB = B.sortxyz();
scalar sumDifference = 0;
for (direction dir = 0; dir < 3; dir++)
{
if (!tmpA.set(dir) || !tmpB.set(dir))
{
continue;
}
scalar cosPhi =
(tmpA[dir] & tmpB[dir])
/(mag(tmpA[dir])*mag(tmpA[dir]) + SMALL);
cosPhi = min(max(cosPhi, -1), 1);
sumDifference += mag(cosPhi - 1);
}
return (sumDifference/3);
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -154,6 +154,9 @@ public:
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
//- Return a quantity of the difference between two triads
scalar diff(const triad& A, const triad& B);
//- Data associated with quaternion type are contiguous
template<>
inline bool contiguous<triad>() {return true;}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -239,6 +239,7 @@ bool Foam::motionSmoother::checkMesh
maxIntSkew,
maxBounSkew,
mesh,
mesh.points(),
mesh.cellCentres(),
mesh.faceCentres(),
mesh.faceAreas(),
@ -481,14 +482,14 @@ bool Foam::motionSmoother::checkMesh
(
readScalar(dict.lookup("minArea", true))
);
const scalar maxIntSkew
(
readScalar(dict.lookup("maxInternalSkewness", true))
);
const scalar maxBounSkew
(
readScalar(dict.lookup("maxBoundarySkewness", true))
);
//const scalar maxIntSkew
//(
// readScalar(dict.lookup("maxInternalSkewness", true))
//);
//const scalar maxBounSkew
//(
// readScalar(dict.lookup("maxBoundarySkewness", true))
//);
const scalar minWeight
(
readScalar(dict.lookup("minFaceWeight", true))
@ -615,27 +616,30 @@ bool Foam::motionSmoother::checkMesh
nWrongFaces = nNewWrongFaces;
}
if (maxIntSkew > 0 || maxBounSkew > 0)
{
meshGeom.checkFaceSkewness
(
report,
maxIntSkew,
maxBounSkew,
checkFaces,
baffles,
&wrongFaces
);
label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
Info<< " faces with skewness > "
<< setw(3) << maxIntSkew
<< " (internal) or " << setw(3) << maxBounSkew
<< " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
//- Note: cannot check the skewness without the points and don't want
// to store them on polyMeshGeometry.
//if (maxIntSkew > 0 || maxBounSkew > 0)
//{
// meshGeom.checkFaceSkewness
// (
// report,
// maxIntSkew,
// maxBounSkew,
// checkFaces,
// baffles,
// &wrongFaces
// );
//
// label nNewWrongFaces = returnReduce(wrongFaces.size(),sumOp<label>());
//
// Info<< " faces with skewness > "
// << setw(3) << maxIntSkew
// << " (internal) or " << setw(3) << maxBounSkew
// << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
//
// nWrongFaces = nNewWrongFaces;
//}
if (minWeight >= 0 && minWeight < 1)
{

View File

@ -29,6 +29,7 @@ License
#include "tetrahedron.H"
#include "syncTools.H"
#include "unitConversion.H"
#include "primitiveMeshTools.H"
namespace Foam
{
@ -286,29 +287,6 @@ Foam::scalar Foam::polyMeshGeometry::checkNonOrtho
}
Foam::scalar Foam::polyMeshGeometry::calcSkewness
(
const point& ownCc,
const point& neiCc,
const point& fc
)
{
scalar dOwn = mag(fc - ownCc);
scalar dNei = mag(fc - neiCc);
point faceIntersection =
ownCc*dNei/(dOwn+dNei+VSMALL)
+ neiCc*dOwn/(dOwn+dNei+VSMALL);
return
mag(fc - faceIntersection)
/ (
mag(neiCc-ownCc)
+ VSMALL
);
}
// Create the neighbour pyramid - it will have positive volume
bool Foam::polyMeshGeometry::checkFaceTet
(
@ -1008,6 +986,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
const scalar internalSkew,
const scalar boundarySkew,
const polyMesh& mesh,
const pointField& points,
const vectorField& cellCentres,
const vectorField& faceCentres,
const vectorField& faceAreas,
@ -1024,14 +1003,8 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Calculate coupled cell centre
pointField neiCc(mesh.nFaces()-mesh.nInternalFaces());
for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
{
neiCc[faceI-mesh.nInternalFaces()] = cellCentres[own[faceI]];
}
syncTools::swapBoundaryFacePositions(mesh, neiCc);
pointField neiCc;
syncTools::swapBoundaryCellPositions(mesh, cellCentres, neiCc);
scalar maxSkew = 0;
@ -1043,11 +1016,16 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
if (mesh.isInternalFace(faceI))
{
scalar skewness = calcSkewness
scalar skewness = primitiveMeshTools::faceSkewness
(
mesh,
points,
faceCentres,
faceAreas,
faceI,
cellCentres[own[faceI]],
cellCentres[nei[faceI]],
faceCentres[faceI]
cellCentres[nei[faceI]]
);
// Check if the skewness vector is greater than the PN vector.
@ -1073,11 +1051,16 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
}
else if (patches[patches.whichPatch(faceI)].coupled())
{
scalar skewness = calcSkewness
scalar skewness = primitiveMeshTools::faceSkewness
(
mesh,
points,
faceCentres,
faceAreas,
faceI,
cellCentres[own[faceI]],
neiCc[faceI-mesh.nInternalFaces()],
faceCentres[faceI]
neiCc[faceI-mesh.nInternalFaces()]
);
// Check if the skewness vector is greater than the PN vector.
@ -1103,21 +1086,17 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
}
else
{
// Boundary faces: consider them to have only skewness error.
// (i.e. treat as if mirror cell on other side)
scalar skewness = primitiveMeshTools::boundaryFaceSkewness
(
mesh,
points,
faceCentres,
faceAreas,
vector faceNormal = faceAreas[faceI];
faceNormal /= mag(faceNormal) + ROOTVSMALL;
faceI,
cellCentres[own[faceI]]
);
vector dOwn = faceCentres[faceI] - cellCentres[own[faceI]];
vector dWall = faceNormal*(faceNormal & dOwn);
point faceIntersection = cellCentres[own[faceI]] + dWall;
scalar skewness =
mag(faceCentres[faceI] - faceIntersection)
/(2*mag(dWall) + ROOTVSMALL);
// Check if the skewness vector is greater than the PN vector.
// This does not cause trouble but is a good indication of a poor
@ -1148,12 +1127,18 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
label face1 = baffles[i].second();
const point& ownCc = cellCentres[own[face0]];
const point& neiCc = cellCentres[own[face1]];
scalar skewness = calcSkewness
scalar skewness = primitiveMeshTools::faceSkewness
(
mesh,
points,
faceCentres,
faceAreas,
face0,
ownCc,
cellCentres[own[face1]],
faceCentres[face0]
neiCc
);
// Check if the skewness vector is greater than the PN vector.
@ -2361,30 +2346,30 @@ bool Foam::polyMeshGeometry::checkFaceTets
}
bool Foam::polyMeshGeometry::checkFaceSkewness
(
const bool report,
const scalar internalSkew,
const scalar boundarySkew,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
) const
{
return checkFaceSkewness
(
report,
internalSkew,
boundarySkew,
mesh_,
cellCentres_,
faceCentres_,
faceAreas_,
checkFaces,
baffles,
setPtr
);
}
//bool Foam::polyMeshGeometry::checkFaceSkewness
//(
// const bool report,
// const scalar internalSkew,
// const scalar boundarySkew,
// const labelList& checkFaces,
// const List<labelPair>& baffles,
// labelHashSet* setPtr
//) const
//{
// return checkFaceSkewness
// (
// report,
// internalSkew,
// boundarySkew,
// mesh_,
// cellCentres_,
// faceCentres_,
// faceAreas_,
// checkFaces,
// baffles,
// setPtr
// );
//}
bool Foam::polyMeshGeometry::checkFaceWeights

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -229,6 +229,7 @@ public:
const scalar internalSkew,
const scalar boundarySkew,
const polyMesh& mesh,
const pointField& points,
const vectorField& cellCentres,
const vectorField& faceCentres,
const vectorField& faceAreas,
@ -372,16 +373,6 @@ public:
labelHashSet* setPtr
) const;
bool checkFaceSkewness
(
const bool report,
const scalar internalSkew,
const scalar boundarySkew,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
) const;
bool checkFaceWeights
(
const bool report,

View File

@ -32,6 +32,7 @@ License
#include "polyTopoChange.H"
#include "globalIndex.H"
#include "PackedBoolList.H"
#include "faceZone.H"
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
@ -492,8 +493,37 @@ Foam::polyMeshFilter::polyMeshFilter(const fvMesh& mesh)
<< edgeReductionFactor_ << nl
<< endl;
Info<< "Collapse faces with reduction factor = " << faceReductionFactor_
<< endl;
if (collapseFacesCoeffDict_.empty())
{
Info<< "Face collapsing is off" << endl;
}
else
{
Info<< "Face collapsing is on" << endl;
Info<< " Initial face length factor = "<< initialFaceLengthFactor_
<< endl;
}
Info<< "Control mesh quality = " << controlMeshQuality_.asText() << endl;
if (controlMeshQuality_)
{
Info<< " Minimum edge length reduction factor = "
<< edgeReductionFactor_ << nl
<< " Minimum face area reduction factor = "
<< faceReductionFactor_ << endl;
Info<< " Maximum number of collapse iterations = " << maxIterations_
<< endl;
Info<< " Maximum number of edge/face reduction factor smoothing "
<< "iterations = " << maxSmoothIters_ << endl;
Info<< " Maximum number of times a point can contribute to bad "
<< "faces across " << nl
<< " collapse iterations = " << maxPointErrorCount_
<< endl;
}
Info<< "Selectively disabling wanted collapses until resulting quality"
<< " satisfies constraints in system/meshQualityDict" << nl
@ -520,6 +550,16 @@ Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
// Maintain the number of times a point has been part of a bad face
labelList pointErrorCount(mesh_.nPoints(), 0);
PackedBoolList newErrorPoint(mesh_.nPoints());
edgeCollapser::checkMeshQuality
(
mesh_,
meshQualityCoeffDict_,
newErrorPoint
);
bool newBadFaces = true;
// Main loop
// ~~~~~~~~~
// It tries and do some collapses, checks the resulting mesh and
@ -529,7 +569,8 @@ Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
while
(
nOuterIterations < maxIterations_
&& nBadFaces > nOriginalBadFaces
//&& nBadFaces > nOriginalBadFaces
&& newBadFaces
)
{
Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
@ -872,6 +913,21 @@ Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
isErrorPoint,
pointErrorCount
);
newBadFaces = false;
forAll(mesh_.points(), pI)
{
if (isErrorPoint[origToCurrentPointMap[pI]])
{
if (!newErrorPoint[pI])
{
newBadFaces = true;
break;
}
}
}
reduce(newBadFaces, orOp<bool>());
}
else
{
@ -1105,83 +1161,330 @@ Foam::label Foam::polyMeshFilter::filterEdges
}
Foam::label Foam::polyMeshFilter::filterIndirectPatchFaces()
Foam::label Foam::polyMeshFilter::filterFaceZone(const faceZone& fZone)
{
newMeshPtr_ = copyMesh(mesh_);
fvMesh& newMesh = newMeshPtr_();
const label nOriginalBadFaces = 0;
label nIterations = 0;
label nBadFaces = 0;
label nBadFaces = labelMax;
label nOuterIterations = 0;
while (true)
minEdgeLen_.resize(mesh_.nEdges(), minLen_);
faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor_);
// Maintain the number of times a point has been part of a bad face
labelList pointErrorCount(mesh_.nPoints(), 0);
// Main loop
// ~~~~~~~~~
// It tries and do some collapses, checks the resulting mesh and
// 'freezes' some edges (by marking in minEdgeLen) and tries again.
// This will iterate ultimately to the situation where every edge is
// frozen and nothing gets collapsed.
while
(
nOuterIterations < maxIterations_
&& nBadFaces > nOriginalBadFaces
)
{
Info<< nl << indent << "Iteration = "
<< nIterations++ << nl << incrIndent << endl;
Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
<< endl;
// Per edge collapse status
PackedBoolList collapseEdge(newMesh.nEdges());
printScalarFieldStats("Edge Filter Factor", minEdgeLen_);
printScalarFieldStats("Face Filter Factor", faceFilterFactor_);
Map<point> collapsePointToLocation(newMesh.nPoints());
// Reset the new mesh to the old mesh
newMeshPtr_ = copyMesh(mesh_);
fvMesh& newMesh = newMeshPtr_();
labelList boundaryPoint(newMesh.nPoints());
scalarField newMeshFaceFilterFactor = faceFilterFactor_;
edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
collapser.markIndirectPatchFaces
(
collapseEdge,
collapsePointToLocation
);
// Merge edge collapses into consistent collapse-network.
// Make sure no cells get collapsed.
List<pointEdgeCollapse> allPointInfo;
const globalIndex globalPoints(newMesh.nPoints());
collapser.consistentCollapse
(
globalPoints,
boundaryPoint,
collapsePointToLocation,
collapseEdge,
allPointInfo
);
label nLocalCollapse = collapseEdge.count();
reduce(nLocalCollapse, sumOp<label>());
Info<< nl << indent << "Collapsing " << nLocalCollapse
<< " edges after synchronisation and PointEdgeWave" << endl;
if (nLocalCollapse == 0)
labelList origToCurrentPointMap(identity(newMesh.nPoints()));
{
break;
label nInnerIterations = 0;
label nPrevLocalCollapse = labelMax;
Info<< incrIndent;
while (true)
{
Info<< nl << indent << "Inner iteration = "
<< nInnerIterations++ << nl << incrIndent << endl;
// Per edge collapse status
PackedBoolList collapseEdge(newMesh.nEdges());
Map<point> collapsePointToLocation(newMesh.nPoints());
// Mark points on boundary
const labelList boundaryPoint = findBoundaryPoints(newMesh);
edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
{
// Collapse faces
labelPair nCollapsedPtEdge = collapser.markFaceZoneEdges
(
fZone,
newMeshFaceFilterFactor,
boundaryPoint,
collapseEdge,
collapsePointToLocation
);
label nCollapsed = 0;
forAll(nCollapsedPtEdge, collapseTypeI)
{
nCollapsed += nCollapsedPtEdge[collapseTypeI];
}
reduce(nCollapsed, sumOp<label>());
Info<< indent
<< "Collapsing " << nCollapsed << " faces"
<< " (to point = "
<< returnReduce
(
nCollapsedPtEdge.first(),
sumOp<label>()
)
<< ", to edge = "
<< returnReduce
(
nCollapsedPtEdge.second(),
sumOp<label>()
)
<< ")" << endl;
if (nCollapsed == 0)
{
Info<< decrIndent;
Info<< decrIndent;
break;
}
}
// Merge edge collapses into consistent collapse-network.
// Make sure no cells get collapsed.
List<pointEdgeCollapse> allPointInfo;
const globalIndex globalPoints(newMesh.nPoints());
collapser.consistentCollapse
(
globalPoints,
boundaryPoint,
collapsePointToLocation,
collapseEdge,
allPointInfo
);
label nLocalCollapse = collapseEdge.count();
reduce(nLocalCollapse, sumOp<label>());
Info<< nl << indent << "Collapsing " << nLocalCollapse
<< " edges after synchronisation and PointEdgeWave" << endl;
if (nLocalCollapse >= nPrevLocalCollapse)
{
Info<< decrIndent;
Info<< decrIndent;
break;
}
else
{
nPrevLocalCollapse = nLocalCollapse;
}
{
// Apply collapses to current mesh
polyTopoChange newMeshMod(newMesh);
// Insert mesh refinement into polyTopoChange.
collapser.setRefinement(allPointInfo, newMeshMod);
Info<< indent << "Apply changes to the current mesh"
<< decrIndent << endl;
// Apply changes to current mesh
autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
(
newMesh,
false
);
const mapPolyMesh& newMap = newMapPtr();
// Update fields
newMesh.updateMesh(newMap);
if (newMap.hasMotionPoints())
{
newMesh.movePoints(newMap.preMotionPoints());
}
// Relabel the boundary points
// labelList newBoundaryPoints(newMesh.nPoints(), -1);
// forAll(newBoundaryPoints, pI)
// {
// const label newToOldptI= map.pointMap()[pI];
// newBoundaryPoints[pI] = boundaryIOPts[newToOldptI];
// }
// boundaryIOPts = newBoundaryPoints;
mapOldMeshFaceFieldToNewMesh
(
newMesh,
newMap.faceMap(),
newMeshFaceFilterFactor
);
updateOldToNewPointMap
(
newMap.reversePointMap(),
origToCurrentPointMap
);
}
}
}
// Apply collapses to current mesh
polyTopoChange newMeshMod(newMesh);
// Insert mesh refinement into polyTopoChange.
collapser.setRefinement(allPointInfo, newMeshMod);
scalarField newMeshMinEdgeLen = minEdgeLen_;
Info<< indent << "Apply changes to the current mesh"
<< decrIndent << endl;
label nInnerIterations = 0;
label nPrevLocalCollapse = labelMax;
// Apply changes to current mesh
autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
(
newMesh,
false
);
const mapPolyMesh& newMap = newMapPtr();
// Update fields
newMesh.updateMesh(newMap);
if (newMap.hasMotionPoints())
while (true)
{
newMesh.movePoints(newMap.preMotionPoints());
Info<< nl << indent << "Inner iteration = "
<< nInnerIterations++ << nl << incrIndent << endl;
// Per edge collapse status
PackedBoolList collapseEdge(newMesh.nEdges());
Map<point> collapsePointToLocation(newMesh.nPoints());
// Mark points on boundary
const labelList boundaryPoint = findBoundaryPoints
(
newMesh//,
// boundaryIOPts
);
edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
// Work out which edges to collapse
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// This is by looking at minEdgeLen (to avoid frozen edges)
// and marking in collapseEdge.
label nSmallCollapsed = collapser.markSmallEdges
(
newMeshMinEdgeLen,
boundaryPoint,
collapseEdge,
collapsePointToLocation
);
reduce(nSmallCollapsed, sumOp<label>());
Info<< indent << "Collapsing " << nSmallCollapsed
<< " small edges" << endl;
// Merge inline edges
label nMerged = collapser.markMergeEdges
(
maxCos_,
boundaryPoint,
collapseEdge,
collapsePointToLocation
);
reduce(nMerged, sumOp<label>());
Info<< indent << "Collapsing " << nMerged << " in line edges"
<< endl;
if (nMerged + nSmallCollapsed == 0)
{
Info<< decrIndent;
break;
}
// Merge edge collapses into consistent collapse-network.
// Make sure no cells get collapsed.
List<pointEdgeCollapse> allPointInfo;
const globalIndex globalPoints(newMesh.nPoints());
collapser.consistentCollapse
(
globalPoints,
boundaryPoint,
collapsePointToLocation,
collapseEdge,
allPointInfo
);
label nLocalCollapse = collapseEdge.count();
reduce(nLocalCollapse, sumOp<label>());
Info<< nl << indent << "Collapsing " << nLocalCollapse
<< " edges after synchronisation and PointEdgeWave" << endl;
if (nLocalCollapse >= nPrevLocalCollapse)
{
Info<< decrIndent;
break;
}
else
{
nPrevLocalCollapse = nLocalCollapse;
}
// Apply collapses to current mesh
polyTopoChange newMeshMod(newMesh);
// Insert mesh refinement into polyTopoChange.
collapser.setRefinement(allPointInfo, newMeshMod);
Info<< indent << "Apply changes to the current mesh"
<< decrIndent << endl;
// Apply changes to current mesh
autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
(
newMesh,
false
);
const mapPolyMesh& newMap = newMapPtr();
// Update fields
newMesh.updateMesh(newMap);
if (newMap.hasMotionPoints())
{
newMesh.movePoints(newMap.preMotionPoints());
}
// Relabel the boundary points
// labelList newBoundaryPoints(newMesh.nPoints(), -1);
// forAll(newBoundaryPoints, pI)
// {
// const label newToOldptI= map.pointMap()[pI];
// newBoundaryPoints[pI] = boundaryIOPts[newToOldptI];
// }
// boundaryIOPts = newBoundaryPoints;
// Synchronise the factors
mapOldMeshEdgeFieldToNewMesh
(
newMesh,
newMap.pointMap(),
newMeshMinEdgeLen
);
updateOldToNewPointMap
(
newMap.reversePointMap(),
origToCurrentPointMap
);
}
// Mesh check
// ~~~~~~~~~~~~~~~~~~
// Do not allow collapses in regions of error.
@ -1201,6 +1504,33 @@ Foam::label Foam::polyMeshFilter::filterIndirectPatchFaces()
<< " Number of marked points : "
<< returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
<< endl;
updatePointErrorCount
(
isErrorPoint,
origToCurrentPointMap,
pointErrorCount
);
checkMeshEdgesAndRelaxEdges
(
newMesh,
origToCurrentPointMap,
isErrorPoint,
pointErrorCount
);
checkMeshFacesAndRelaxEdges
(
newMesh,
origToCurrentPointMap,
isErrorPoint,
pointErrorCount
);
}
else
{
return -1;
}
}

View File

@ -52,6 +52,7 @@ namespace Foam
class polyMesh;
class fvMesh;
class PackedBoolList;
class faceZone;
/*---------------------------------------------------------------------------*\
Class polyMeshFilter Declaration
@ -238,7 +239,7 @@ public:
label filterEdges(const label nOriginalBadFaces);
//- Filter all faces that are in the face zone indirectPatchFaces
label filterIndirectPatchFaces();
label filterFaceZone(const faceZone& fZone);
};

View File

@ -32,6 +32,7 @@ License
#include "globalIndex.H"
#include "removePoints.H"
#include "motionSmoother.H"
#include "OFstream.H"
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
@ -1242,7 +1243,19 @@ Foam::edgeCollapser::edgeCollapser
(
dict.lookupOrDefault<scalar>("allowEarlyCollapseCoeff", 0)
)
{}
{
if (debug)
{
Info<< "Edge Collapser Settings:" << nl
<< " Guard Fraction = " << guardFraction_ << nl
<< " Max collapse face to point side length = "
<< maxCollapseFaceToPointSideLengthCoeff_ << nl
<< " " << (allowEarlyCollapseToPoint_ ? "Allow" : "Do not allow")
<< " early collapse to point" << nl
<< " Early collapse coeff = " << allowEarlyCollapseCoeff_
<< endl;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -2017,94 +2030,151 @@ Foam::labelPair Foam::edgeCollapser::markSmallSliverFaces
}
void Foam::edgeCollapser::markIndirectPatchFaces
Foam::labelPair Foam::edgeCollapser::markFaceZoneEdges
(
const faceZone& fZone,
const scalarField& faceFilterFactor,
const labelList& pointPriority,
PackedBoolList& collapseEdge,
Map<point>& collapsePointToLocation
) const
{
const faceZone& indirectFaceZone = mesh_.faceZones()["indirectPatchFaces"];
const faceList& faces = mesh_.faces();
const edgeList& edges = mesh_.edges();
const pointField& points = mesh_.points();
const labelListList& edgeFaces = mesh_.edgeFaces();
const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
const scalarField targetFaceSizes = calcTargetFaceSizes();
forAll(edges, eI)
// Calculate number of faces that will be collapsed to a point or an edge
label nCollapseToPoint = 0;
label nCollapseToEdge = 0;
forAll(faces, fI)
{
const edge& e = edges[eI];
const labelList& eFaces = edgeFaces[eI];
bool keepEdge = false;
label nInternalFaces = 0;
label nPatchFaces = 0;
label nIndirectFaces = 0;
bool coupled = false;
forAll(eFaces, eFaceI)
if (fZone.whichFace(fI) == -1)
{
const label eFaceIndex = eFaces[eFaceI];
if (mesh_.isInternalFace(eFaceIndex))
{
nInternalFaces++;
}
else
{
const label patchIndex = bMesh.whichPatch(eFaceIndex);
const polyPatch& pPatch = bMesh[patchIndex];
if (pPatch.coupled())
{
coupled = true;
nInternalFaces++;
}
else
{
// Keep the edge if an attached face is not in the face zone
if (indirectFaceZone.whichFace(eFaceIndex) == -1)
{
nPatchFaces++;
}
else
{
nIndirectFaces++;
}
}
}
continue;
}
if (eFaces.size() != nInternalFaces + nPatchFaces + nIndirectFaces)
const face& f = faces[fI];
if (faceFilterFactor[fI] == 0)
{
Pout<< eFaces.size() << " ("
<< nInternalFaces << "/" << nPatchFaces << "/" << nIndirectFaces
<< ")" << endl;
continue;
}
if
collapseType flagCollapseFace = collapseFace
(
eFaces.size() == nInternalFaces
|| nIndirectFaces < (coupled ? 1 : 2)
)
pointPriority,
f,
fI,
targetFaceSizes[fI],
collapseEdge,
collapsePointToLocation,
faceFilterFactor
);
if (flagCollapseFace == noCollapse)
{
keepEdge = true;
continue;
}
if (!keepEdge)
else if (flagCollapseFace == toPoint)
{
collapseEdge[eI] = true;
const Foam::point collapsePoint =
0.5*(points[e.end()] + points[e.start()]);
collapsePointToLocation.insert(e.start(), collapsePoint);
collapsePointToLocation.insert(e.end(), collapsePoint);
nCollapseToPoint++;
}
else if (flagCollapseFace == toEdge)
{
nCollapseToEdge++;
}
else
{
FatalErrorIn("collapseFaces(const polyMesh&, List<labelPair>&)")
<< "Face is marked to be collapsed to " << flagCollapseFace
<< ". Currently can only collapse to point/edge."
<< abort(FatalError);
}
}
return labelPair(nCollapseToPoint, nCollapseToEdge);
// const edgeList& edges = mesh_.edges();
// const pointField& points = mesh_.points();
// const labelListList& edgeFaces = mesh_.edgeFaces();
// const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
//
// forAll(edges, eI)
// {
// const edge& e = edges[eI];
//
// const labelList& eFaces = edgeFaces[eI];
//
// bool keepEdge = false;
//
// label nInternalFaces = 0;
// label nPatchFaces = 0;
// label nIndirectFaces = 0;
//
// bool coupled = false;
//
// forAll(eFaces, eFaceI)
// {
// const label eFaceIndex = eFaces[eFaceI];
//
// if (mesh_.isInternalFace(eFaceIndex))
// {
// nInternalFaces++;
// }
// else
// {
// const label patchIndex = bMesh.whichPatch(eFaceIndex);
// const polyPatch& pPatch = bMesh[patchIndex];
//
// if (pPatch.coupled())
// {
// coupled = true;
// nInternalFaces++;
// }
// else
// {
// // Keep the edge if an attached face is not in the zone
// if (fZone.whichFace(eFaceIndex) == -1)
// {
// nPatchFaces++;
// }
// else
// {
// nIndirectFaces++;
// }
// }
// }
// }
//
// if (eFaces.size() != nInternalFaces + nPatchFaces + nIndirectFaces)
// {
// Pout<< eFaces.size() << " ("
// << nInternalFaces << "/" << nPatchFaces << "/"
// << nIndirectFaces << ")" << endl;
// }
//
// if
// (
// eFaces.size() == nInternalFaces
// || nIndirectFaces < (coupled ? 1 : 2)
// )
// {
// keepEdge = true;
// }
//
// if (!keepEdge)
// {
// collapseEdge[eI] = true;
//
// const Foam::point collapsePoint =
// 0.5*(points[e.end()] + points[e.start()]);
//
// collapsePointToLocation.insert(e.start(), collapsePoint);
// collapsePointToLocation.insert(e.end(), collapsePoint);
// }
// }
// OFstream str
// (
// mesh_.time().path()

View File

@ -339,8 +339,11 @@ public:
) const;
//- Marks edges in the faceZone indirectPatchFaces for collapse
void markIndirectPatchFaces
labelPair markFaceZoneEdges
(
const faceZone& fZone,
const scalarField& faceFilterFactor,
const labelList& pointPriority,
PackedBoolList& collapseEdge,
Map<point>& collapsePointToLocation
) const;

View File

@ -68,6 +68,19 @@ namespace Foam
"multiple",
"none"
};
template<>
const char* Foam::NamedEnum
<
Foam::extendedFeatureEdgeMesh::sideVolumeType,
4
>::names[] =
{
"inside",
"outside",
"both",
"neither"
};
}
const Foam::NamedEnum<Foam::extendedFeatureEdgeMesh::pointStatus, 4>
@ -76,6 +89,9 @@ const Foam::NamedEnum<Foam::extendedFeatureEdgeMesh::pointStatus, 4>
const Foam::NamedEnum<Foam::extendedFeatureEdgeMesh::edgeStatus, 6>
Foam::extendedFeatureEdgeMesh::edgeStatusNames_;
const Foam::NamedEnum<Foam::extendedFeatureEdgeMesh::sideVolumeType, 4>
Foam::extendedFeatureEdgeMesh::sideVolumeTypeNames_;
Foam::scalar Foam::extendedFeatureEdgeMesh::cosNormalAngleTol_ =
Foam::cos(degToRad(0.1));
@ -106,7 +122,9 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh(const IOobject& io)
openStart_(0),
multipleStart_(0),
normals_(0),
normalVolumeTypes_(0),
edgeDirections_(0),
normalDirections_(0),
edgeNormals_(0),
featurePointNormals_(0),
featurePointEdges_(0),
@ -144,6 +162,8 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh(const IOobject& io)
>> openStart_
>> multipleStart_
>> normals_
>> normalVolumeTypes_
>> normalDirections_
>> edgeNormals_
>> featurePointNormals_
>> featurePointEdges_
@ -165,7 +185,7 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh(const IOobject& io)
edgeDirections_[eI] = eds[eI].vec(pts);
}
edgeDirections_ /= mag(edgeDirections_);
edgeDirections_ /= (mag(edgeDirections_) + SMALL);
}
}
@ -196,7 +216,9 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
openStart_(fem.openStart()),
multipleStart_(fem.multipleStart()),
normals_(fem.normals()),
normalVolumeTypes_(fem.normalVolumeTypes()),
edgeDirections_(fem.edgeDirections()),
normalDirections_(fem.normalDirections()),
edgeNormals_(fem.edgeNormals()),
featurePointNormals_(fem.featurePointNormals()),
featurePointEdges_(fem.featurePointEdges()),
@ -224,7 +246,9 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
openStart_(0),
multipleStart_(0),
normals_(0),
normalVolumeTypes_(0),
edgeDirections_(0),
normalDirections_(0),
edgeNormals_(0),
featurePointNormals_(0),
featurePointEdges_(0),
@ -263,7 +287,9 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
openStart_(-1),
multipleStart_(-1),
normals_(0),
normalVolumeTypes_(0),
edgeDirections_(0),
normalDirections_(0),
edgeNormals_(0),
featurePointNormals_(0),
featurePointEdges_(0),
@ -287,6 +313,36 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
regionFeatureEdges,
featurePoints
);
const labelListList& edgeFaces = surf.edgeFaces();
normalVolumeTypes_.setSize(normals_.size());
// Noting when the normal of a face has been used so not to duplicate
labelList faceMap(surf.size(), -1);
label nAdded = 0;
forAll(featureEdges, i)
{
label sFEI = featureEdges[i];
// Pick up the faces adjacent to the feature edge
const labelList& eFaces = edgeFaces[sFEI];
forAll(eFaces, j)
{
label eFI = eFaces[j];
// Check to see if the points have been already used
if (faceMap[eFI] == -1)
{
normalVolumeTypes_[nAdded++] = INSIDE;
faceMap[eFI] = nAdded - 1;
}
}
}
}
@ -309,7 +365,9 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
openStart_(-1),
multipleStart_(-1),
normals_(0),
normalVolumeTypes_(0),
edgeDirections_(0),
normalDirections_(0),
edgeNormals_(0),
featurePointNormals_(0),
featurePointEdges_(0),
@ -341,7 +399,9 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
label openStart,
label multipleStart,
const vectorField& normals,
const PackedList<2>& normalVolumeTypes,
const vectorField& edgeDirections,
const labelListList& normalDirections,
const labelListList& edgeNormals,
const labelListList& featurePointNormals,
const labelListList& featurePointEdges,
@ -358,7 +418,9 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
openStart_(openStart),
multipleStart_(multipleStart),
normals_(normals),
normalVolumeTypes_(normalVolumeTypes),
edgeDirections_(edgeDirections),
normalDirections_(normalDirections),
edgeNormals_(edgeNormals),
featurePointNormals_(featurePointNormals),
featurePointEdges_(featurePointEdges),
@ -1300,6 +1362,10 @@ bool Foam::extendedFeatureEdgeMesh::writeData(Ostream& os) const
<< multipleStart_ << nl
<< "// normals" << nl
<< normals_ << nl
<< "// normal volume types" << nl
<< normalVolumeTypes_ << nl
<< "// normalDirections" << nl
<< normalDirections_ << nl
<< "// edgeNormals" << nl
<< edgeNormals_ << nl
<< "// featurePointNormals" << nl

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -61,6 +61,7 @@ SourceFiles
#include "treeDataEdge.H"
#include "treeDataPoint.H"
#include "PrimitivePatch.H"
#include "PackedList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -108,6 +109,17 @@ public:
static const Foam::NamedEnum<edgeStatus, 6> edgeStatusNames_;
//- Normals point to the outside
enum sideVolumeType
{
INSIDE = 0, // mesh inside
OUTSIDE = 1, // mesh outside
BOTH = 2, // e.g. a baffle
NEITHER = 3 // not sure when this may be used
};
static const Foam::NamedEnum<sideVolumeType, 4> sideVolumeTypeNames_;
private:
@ -150,9 +162,17 @@ private:
// points and edges, unsorted
vectorField normals_;
//-
PackedList<2> normalVolumeTypes_;
//- Flat and open edges require the direction of the edge
vectorField edgeDirections_;
//- Starting directions for the edges.
// This vector points to the half of the plane defined by the first
// edge normal.
labelListList normalDirections_;
//- Indices of the normals that are adjacent to the feature edges
labelListList edgeNormals_;
@ -267,7 +287,9 @@ public:
label openStart,
label multipleStart,
const vectorField& normals,
const PackedList<2>& normalVolumeTypes,
const vectorField& edgeDirections,
const labelListList& normalDirections,
const labelListList& edgeNormals,
const labelListList& featurePointNormals,
const labelListList& featurePointEdges,
@ -369,9 +391,15 @@ public:
// and points
inline const vectorField& normals() const;
//- Return
inline const PackedList<2>& normalVolumeTypes() const;
//- Return the edgeDirection vectors
inline const vectorField& edgeDirections() const;
//-
inline const labelListList& normalDirections() const;
//- Return the direction of edgeI, pointing away from ptI
inline vector edgeDirection(label edgeI, label ptI) const;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -90,12 +90,24 @@ inline const Foam::vectorField& Foam::extendedFeatureEdgeMesh::normals() const
return normals_;
}
inline const Foam::PackedList<2>&
Foam::extendedFeatureEdgeMesh::normalVolumeTypes() const
{
return normalVolumeTypes_;
}
inline const Foam::vectorField& Foam::extendedFeatureEdgeMesh::edgeDirections()
const
{
return edgeDirections_;
}
inline const Foam::labelListList&
Foam::extendedFeatureEdgeMesh::normalDirections() const
{
return normalDirections_;
}
inline Foam::vector Foam::extendedFeatureEdgeMesh::edgeDirection
(
@ -205,10 +217,7 @@ inline const Foam::labelList& Foam::extendedFeatureEdgeMesh::regionEdges() const
inline Foam::extendedFeatureEdgeMesh::pointStatus
Foam::extendedFeatureEdgeMesh::getPointStatus
(
label ptI
) const
Foam::extendedFeatureEdgeMesh::getPointStatus(label ptI) const
{
if (ptI < concaveStart_)
{
@ -230,10 +239,7 @@ Foam::extendedFeatureEdgeMesh::getPointStatus
inline Foam::extendedFeatureEdgeMesh::edgeStatus
Foam::extendedFeatureEdgeMesh::getEdgeStatus
(
label edgeI
) const
Foam::extendedFeatureEdgeMesh::getEdgeStatus(label edgeI) const
{
if (edgeI < internalStart_)
{

View File

@ -27,6 +27,8 @@ License
#include "ListListOps.H"
#include "unitConversion.H"
#include "PackedBoolList.H"
#include "PatchTools.H"
#include "searchableBox.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -41,9 +43,9 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
{
const pointField& sFeatLocalPts(surf.localPoints());
const edgeList& sFeatEds(surf.edges());
const labelListList& edgeFaces = surf.edgeFaces();
const labelListList edgeFaces = PatchTools::sortedEdgeFaces(surf);
const vectorField& faceNormals = surf.faceNormals();
const labelListList& pointEdges = surf.pointEdges();
const labelListList pointEdges = PatchTools::sortedPointEdges(surf);
// Extract and reorder the data from surfaceFeatures
@ -59,6 +61,7 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
DynamicList<vector> norms;
vectorField edgeDirections(nFeatEds);
labelListList edgeNormals(nFeatEds);
labelListList normalDirections(nFeatEds);
DynamicList<label> regionEdges;
// Keep track of the ordered feature point feature edges
@ -103,7 +106,9 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
edgeMap[sFEI] = i;
const edge& fE(sFeatEds[sFEI]);
const edge& fE = sFeatEds[sFEI];
edgeDirections[i] = fE.vec(sFeatLocalPts);
// Check to see if the points have been already used
if (pointMap[fE.start()] == -1)
@ -128,6 +133,7 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
const labelList& eFaces = edgeFaces[sFEI];
edgeNormals[i].setSize(eFaces.size());
normalDirections[i].setSize(eFaces.size());
forAll(eFaces, j)
{
@ -142,6 +148,22 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
}
edgeNormals[i][j] = faceMap[eFI];
const vector cross = (faceNormals[eFI] ^ edgeDirections[i]);
const vector fC0tofE0 =
surf[eFI].centre(surf.points())
- sFeatLocalPts[fE.start()];
normalDirections[i][j] =
(
(
(cross/(mag(cross) + VSMALL))
& (fC0tofE0/(mag(fC0tofE0)+ VSMALL))
)
> 0.0
? 1
: -1
);
}
vector fC0tofC1(vector::zero);
@ -155,8 +177,6 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
edStatus[i] = classifyEdge(norms, edgeNormals[i], fC0tofC1);
edgeDirections[i] = fE.vec(sFeatLocalPts);
if (isRegionFeatureEdge[i])
{
regionEdges.append(i);
@ -175,6 +195,7 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
forAll(fpfe, eI)
{
const label oldEdgeIndex = fpfe[eI];
const label newFeatureEdgeIndex = edgeMap[oldEdgeIndex];
if (newFeatureEdgeIndex != -1)
@ -187,7 +208,6 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
}
// Reorder the edges by classification
List<DynamicList<label> > allEds(nEdgeTypes);
DynamicList<label>& externalEds(allEds[0]);
@ -257,6 +277,7 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
inplaceReorder(edMap, edStatus);
inplaceReorder(edMap, edgeDirections);
inplaceReorder(edMap, edgeNormals);
inplaceReorder(edMap, normalDirections);
inplaceRenumber(edMap, regionEdges);
forAll(featurePointFeatureEdges, pI)
@ -272,11 +293,13 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
// Initialise sorted edge related data
edgeDirections_ = edgeDirections/(mag(edgeDirections) + VSMALL);
edgeNormals_ = edgeNormals;
normalDirections_ = normalDirections;
regionEdges_ = regionEdges;
// Normals are all now found and indirectly addressed, can also be stored
normals_ = vectorField(norms);
// Reorder the feature points by classification
List<DynamicList<label> > allPts(3);
@ -354,7 +377,6 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
// renumbered edges
reset(xferMove(pts), xferMove(eds));
// Generate the featurePointNormals
labelListList featurePointNormals(nonFeatureStart_);

View File

@ -18,6 +18,7 @@ $(basicFvPatches)/generic/genericFvPatch.C
constraintFvPatches = $(fvPatches)/constraint
$(constraintFvPatches)/cyclic/cyclicFvPatch.C
$(constraintFvPatches)/cyclicAMI/cyclicAMIFvPatch.C
$(constraintFvPatches)/cyclicACMI/cyclicACMIFvPatch.C
$(constraintFvPatches)/cyclicSlip/cyclicSlipFvPatch.C
$(constraintFvPatches)/empty/emptyFvPatch.C
$(constraintFvPatches)/nonuniformTransformCyclic/nonuniformTransformCyclicFvPatch.C
@ -109,6 +110,7 @@ $(basicFvPatchFields)/zeroGradient/zeroGradientFvPatchFields.C
constraintFvPatchFields = $(fvPatchFields)/constraint
$(constraintFvPatchFields)/cyclic/cyclicFvPatchFields.C
$(constraintFvPatchFields)/cyclicAMI/cyclicAMIFvPatchFields.C
$(constraintFvPatchFields)/cyclicACMI/cyclicACMIFvPatchFields.C
$(constraintFvPatchFields)/cyclicSlip/cyclicSlipFvPatchFields.C
$(constraintFvPatchFields)/empty/emptyFvPatchFields.C
$(constraintFvPatchFields)/jumpCyclic/jumpCyclicFvPatchFields.C
@ -203,6 +205,7 @@ $(basicFvsPatchFields)/sliced/slicedFvsPatchFields.C
constraintFvsPatchFields = $(fvsPatchFields)/constraint
$(constraintFvsPatchFields)/cyclic/cyclicFvsPatchFields.C
$(constraintFvsPatchFields)/cyclicAMI/cyclicAMIFvsPatchFields.C
$(constraintFvsPatchFields)/cyclicACMI/cyclicACMIFvsPatchFields.C
$(constraintFvsPatchFields)/cyclicSlip/cyclicSlipFvsPatchFields.C
$(constraintFvsPatchFields)/empty/emptyFvsPatchFields.C
$(constraintFvsPatchFields)/nonuniformTransformCyclic/nonuniformTransformCyclicFvsPatchFields.C

View File

@ -130,8 +130,8 @@ void Foam::porosityModels::DarcyForchheimer::calcForce
vectorField& force
) const
{
scalarField Udiag(U.size());
vectorField Usource(U.size());
scalarField Udiag(U.size(), 0.0);
vectorField Usource(U.size(), vector::zero);
const scalarField& V = mesh_.V();
apply(Udiag, Usource, V, rho, mu, U);

View File

@ -183,8 +183,8 @@ void Foam::porosityModels::fixedCoeff::calcForce
vectorField& force
) const
{
scalarField Udiag(U.size());
vectorField Usource(U.size());
scalarField Udiag(U.size(), 0.0);
vectorField Usource(U.size(), vector::zero);
const scalarField& V = mesh_.V();
scalar rhoRef = readScalar(coeffs_.lookup("rhoRef"));

View File

@ -53,7 +53,7 @@ namespace porosityModels
{
/*---------------------------------------------------------------------------*\
Class fixedCoeff Declaration
Class fixedCoeff Declaration
\*---------------------------------------------------------------------------*/
class fixedCoeff

View File

@ -155,7 +155,7 @@ Foam::tmp<Foam::vectorField> Foam::porosityModel::porosityModel::force
}
void Foam::porosityModel::porosityModel::addResistance
void Foam::porosityModel::addResistance
(
fvVectorMatrix& UEqn
) const
@ -169,7 +169,7 @@ void Foam::porosityModel::porosityModel::addResistance
}
void Foam::porosityModel::porosityModel::addResistance
void Foam::porosityModel::addResistance
(
fvVectorMatrix& UEqn,
const volScalarField& rho,
@ -185,7 +185,7 @@ void Foam::porosityModel::porosityModel::addResistance
}
void Foam::porosityModel::porosityModel::addResistance
void Foam::porosityModel::addResistance
(
const fvVectorMatrix& UEqn,
volTensorField& AU,
@ -209,14 +209,14 @@ void Foam::porosityModel::porosityModel::addResistance
}
bool Foam::porosityModel::porosityModel::movePoints()
bool Foam::porosityModel::movePoints()
{
// no updates necessary; all member data independent of mesh
return true;
}
void Foam::porosityModel::porosityModel::updateMesh(const mapPolyMesh& mpm)
void Foam::porosityModel::updateMesh(const mapPolyMesh& mpm)
{
// no updates necessary; all member data independent of mesh
}

View File

@ -74,7 +74,7 @@ void Foam::porosityModels::powerLaw::calcForce
vectorField& force
) const
{
scalarField Udiag(U.size());
scalarField Udiag(U.size(), 0.0);
const scalarField& V = mesh_.V();
apply(Udiag, V, rho, U);

View File

@ -0,0 +1,411 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cyclicACMIFvPatchField.H"
#include "transformField.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
Foam::cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField
(
const fvPatch& p,
const DimensionedField<Type, volMesh>& iF
)
:
cyclicACMILduInterfaceField(),
coupledFvPatchField<Type>(p, iF),
cyclicACMIPatch_(refCast<const cyclicACMIFvPatch>(p))
{}
template<class Type>
Foam::cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField
(
const cyclicACMIFvPatchField<Type>& ptf,
const fvPatch& p,
const DimensionedField<Type, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
cyclicACMILduInterfaceField(),
coupledFvPatchField<Type>(ptf, p, iF, mapper),
cyclicACMIPatch_(refCast<const cyclicACMIFvPatch>(p))
{
if (!isA<cyclicACMIFvPatch>(this->patch()))
{
FatalErrorIn
(
"cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField"
"("
"const cyclicACMIFvPatchField<Type>& ,"
"const fvPatch&, "
"const DimensionedField<Type, volMesh>&, "
"const fvPatchFieldMapper&"
")"
) << " patch type '" << p.type()
<< "' not constraint type '" << typeName << "'"
<< "\n for patch " << p.name()
<< " of field " << this->dimensionedInternalField().name()
<< " in file " << this->dimensionedInternalField().objectPath()
<< exit(FatalIOError);
}
}
template<class Type>
Foam::cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField
(
const fvPatch& p,
const DimensionedField<Type, volMesh>& iF,
const dictionary& dict
)
:
cyclicACMILduInterfaceField(),
coupledFvPatchField<Type>(p, iF, dict),
cyclicACMIPatch_(refCast<const cyclicACMIFvPatch>(p))
{
if (!isA<cyclicACMIFvPatch>(p))
{
FatalIOErrorIn
(
"cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField"
"("
"const fvPatch&, "
"const DimensionedField<Type, volMesh>&, "
"const dictionary&"
")",
dict
) << " patch type '" << p.type()
<< "' not constraint type '" << typeName << "'"
<< "\n for patch " << p.name()
<< " of field " << this->dimensionedInternalField().name()
<< " in file " << this->dimensionedInternalField().objectPath()
<< exit(FatalIOError);
}
if (!dict.found("value") && this->coupled())
{
this->evaluate(Pstream::blocking);
}
}
template<class Type>
Foam::cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField
(
const cyclicACMIFvPatchField<Type>& ptf
)
:
cyclicACMILduInterfaceField(),
coupledFvPatchField<Type>(ptf),
cyclicACMIPatch_(ptf.cyclicACMIPatch_)
{}
template<class Type>
Foam::cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField
(
const cyclicACMIFvPatchField<Type>& ptf,
const DimensionedField<Type, volMesh>& iF
)
:
cyclicACMILduInterfaceField(),
coupledFvPatchField<Type>(ptf, iF),
cyclicACMIPatch_(ptf.cyclicACMIPatch_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
bool Foam::cyclicACMIFvPatchField<Type>::coupled() const
{
return cyclicACMIPatch_.coupled();
}
template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::cyclicACMIFvPatchField<Type>::patchNeighbourField() const
{
const Field<Type>& iField = this->internalField();
const labelUList& nbrFaceCellsCoupled =
cyclicACMIPatch_.cyclicACMIPatch().neighbPatch().faceCells();
const labelUList& nbrFaceCellsNonOverlap =
cyclicACMIPatch_.cyclicACMIPatch().nonOverlapPatch().faceCells();
Field<Type> pnfCoupled(iField, nbrFaceCellsCoupled);
Field<Type> pnfNonOverlap(iField, nbrFaceCellsNonOverlap);
tmp<Field<Type> > tpnf
(
new Field<Type>
(
cyclicACMIPatch_.interpolate
(
pnfCoupled,
pnfNonOverlap
)
)
);
if (doTransform())
{
tpnf() = transform(forwardT(), tpnf());
}
return tpnf;
}
template<class Type>
const Foam::cyclicACMIFvPatchField<Type>&
Foam::cyclicACMIFvPatchField<Type>::neighbourPatchField() const
{
const GeometricField<Type, fvPatchField, volMesh>& fld =
static_cast<const GeometricField<Type, fvPatchField, volMesh>&>
(
this->internalField()
);
return refCast<const cyclicACMIFvPatchField<Type> >
(
fld.boundaryField()[cyclicACMIPatch_.neighbPatchID()]
);
}
template<class Type>
const Foam::fvPatchField<Type>&
Foam::cyclicACMIFvPatchField<Type>::nonOverlapPatchField() const
{
const GeometricField<Type, fvPatchField, volMesh>& fld =
static_cast<const GeometricField<Type, fvPatchField, volMesh>&>
(
this->internalField()
);
return fld.boundaryField()[cyclicACMIPatch_.nonOverlapPatchID()];
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
(
scalarField& result,
const scalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
) const
{
// note: only applying coupled contribution
const labelUList& nbrFaceCellsCoupled =
cyclicACMIPatch_.cyclicACMIPatch().neighbPatch().faceCells();
scalarField pnf(psiInternal, nbrFaceCellsCoupled);
// Transform according to the transformation tensors
transformCoupleField(pnf, cmpt);
const labelUList& faceCells = cyclicACMIPatch_.faceCells();
pnf = cyclicACMIPatch_.interpolate(pnf);
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
(
Field<Type>& result,
const Field<Type>& psiInternal,
const scalarField& coeffs,
const Pstream::commsTypes
) const
{
// note: only applying coupled contribution
const labelUList& nbrFaceCellsCoupled =
cyclicACMIPatch_.cyclicACMIPatch().neighbPatch().faceCells();
Field<Type> pnf(psiInternal, nbrFaceCellsCoupled);
// Transform according to the transformation tensors
transformCoupleField(pnf);
const labelUList& faceCells = cyclicACMIPatch_.faceCells();
pnf = cyclicACMIPatch_.interpolate(pnf);
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::cyclicACMIFvPatchField<Type>::snGrad
(
const scalarField& deltaCoeffs
) const
{
// note: only applying coupled contribution
return coupledFvPatchField<Type>::snGrad(deltaCoeffs);
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::updateCoeffs()
{
const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
const fvPatchField<Type>& npf = nonOverlapPatchField();
const_cast<fvPatchField<Type>&>(npf).updateCoeffs(mask);
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::evaluate
(
const Pstream::commsTypes comms
)
{
// blend contrubutions from the coupled and non-overlap patches
const fvPatchField<Type>& npf = nonOverlapPatchField();
const_cast<fvPatchField<Type>&>(npf).evaluate();
coupledFvPatchField<Type>::evaluate(comms);
const Field<Type>& cpf = *this;
const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
Field<Type>::operator=(mask*cpf + (1.0 - mask)*npf);
fvPatchField<Type>::evaluate();
}
template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::cyclicACMIFvPatchField<Type>::valueInternalCoeffs
(
const tmp<scalarField>& w
) const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which area already scaled
return coupledFvPatchField<Type>::valueInternalCoeffs(w);
}
template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::cyclicACMIFvPatchField<Type>::valueBoundaryCoeffs
(
const tmp<scalarField>& w
) const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which area already scaled
return coupledFvPatchField<Type>::valueBoundaryCoeffs(w);
}
template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::cyclicACMIFvPatchField<Type>::gradientInternalCoeffs
(
const scalarField& deltaCoeffs
) const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which area already scaled
return coupledFvPatchField<Type>::gradientInternalCoeffs(deltaCoeffs);
}
template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::cyclicACMIFvPatchField<Type>::gradientInternalCoeffs() const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which area already scaled
return coupledFvPatchField<Type>::gradientInternalCoeffs();
}
template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::cyclicACMIFvPatchField<Type>::gradientBoundaryCoeffs
(
const scalarField& deltaCoeffs
) const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which area already scaled
return coupledFvPatchField<Type>::gradientBoundaryCoeffs(deltaCoeffs);
}
template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::cyclicACMIFvPatchField<Type>::gradientBoundaryCoeffs() const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which area already scaled
return coupledFvPatchField<Type>::gradientBoundaryCoeffs();
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::manipulateMatrix
(
fvMatrix<Type>& matrix
)
{
// blend contrubutions from the coupled and non-overlap patches
const fvPatchField<Type>& npf = nonOverlapPatchField();
const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
const_cast<fvPatchField<Type>&>(npf).manipulateMatrix(matrix, mask);
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::write(Ostream& os) const
{
fvPatchField<Type>::write(os);
this->writeEntry("value", os);
}
// ************************************************************************* //

View File

@ -0,0 +1,304 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::cyclicACMIFvPatchField
Group
grpCoupledBoundaryConditions
Description
This boundary condition enforces a cyclic condition between a pair of
boundaries, whereby communication between the patches is performed using
an arbitrarily coupled mesh interface (ACMI) interpolation.
\heading Patch usage
Example of the boundary condition specification:
\verbatim
myPatch
{
type cyclicACMI;
}
\endverbatim
SeeAlso
Foam::AMIInterpolation
SourceFiles
cyclicACMIFvPatchField.C
\*---------------------------------------------------------------------------*/
#ifndef cyclicACMIFvPatchField_H
#define cyclicACMIFvPatchField_H
#include "coupledFvPatchField.H"
#include "cyclicACMILduInterfaceField.H"
#include "cyclicACMIFvPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class cyclicACMIFvPatchField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class cyclicACMIFvPatchField
:
virtual public cyclicACMILduInterfaceField,
public coupledFvPatchField<Type>
{
// Private data
//- Local reference cast into the cyclic patch
const cyclicACMIFvPatch& cyclicACMIPatch_;
// Private Member Functions
//- Return neighbour side field given internal fields
template<class Type2>
tmp<Field<Type2> > neighbourSideField
(
const Field<Type2>&
) const;
public:
//- Runtime type information
TypeName(cyclicACMIFvPatch::typeName_());
// Constructors
//- Construct from patch and internal field
cyclicACMIFvPatchField
(
const fvPatch&,
const DimensionedField<Type, volMesh>&
);
//- Construct from patch, internal field and dictionary
cyclicACMIFvPatchField
(
const fvPatch&,
const DimensionedField<Type, volMesh>&,
const dictionary&
);
//- Construct by mapping given cyclicACMIFvPatchField onto a new patch
cyclicACMIFvPatchField
(
const cyclicACMIFvPatchField<Type>&,
const fvPatch&,
const DimensionedField<Type, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
cyclicACMIFvPatchField(const cyclicACMIFvPatchField<Type>&);
//- Construct and return a clone
virtual tmp<fvPatchField<Type> > clone() const
{
return tmp<fvPatchField<Type> >
(
new cyclicACMIFvPatchField<Type>(*this)
);
}
//- Construct as copy setting internal field reference
cyclicACMIFvPatchField
(
const cyclicACMIFvPatchField<Type>&,
const DimensionedField<Type, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchField<Type> > clone
(
const DimensionedField<Type, volMesh>& iF
) const
{
return tmp<fvPatchField<Type> >
(
new cyclicACMIFvPatchField<Type>(*this, iF)
);
}
// Member functions
// Access
//- Return local reference cast into the cyclic AMI patch
const cyclicACMIFvPatch& cyclicACMIPatch() const
{
return cyclicACMIPatch_;
}
// Evaluation functions
//- Return true if coupled. Note that the underlying patch
// is not coupled() - the points don't align.
virtual bool coupled() const;
//- Return neighbour coupled internal cell data
virtual tmp<Field<Type> > patchNeighbourField() const;
//- Return reference to neighbour patchField
const cyclicACMIFvPatchField<Type>& neighbourPatchField() const;
//- Return reference to non-overlapping patchField
const fvPatchField<Type>& nonOverlapPatchField() const;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
const scalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
Field<Type>&,
const Field<Type>&,
const scalarField&,
const Pstream::commsTypes commsType
) const;
//- Return patch-normal gradient
virtual tmp<Field<Type> > snGrad
(
const scalarField& deltaCoeffs
) const;
//- Update the coefficients associated with the patch field
void updateCoeffs();
//- Evaluate the patch field
virtual void evaluate
(
const Pstream::commsTypes commsType
);
//- Return the matrix diagonal coefficients corresponding to the
// evaluation of the value of this patchField with given weights
virtual tmp<Field<Type> > valueInternalCoeffs
(
const tmp<scalarField>&
) const;
//- Return the matrix source coefficients corresponding to the
// evaluation of the value of this patchField with given weights
virtual tmp<Field<Type> > valueBoundaryCoeffs
(
const tmp<scalarField>&
) const;
//- Return the matrix diagonal coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type> > gradientInternalCoeffs
(
const scalarField& deltaCoeffs
) const;
//- Return the matrix diagonal coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type> > gradientInternalCoeffs() const;
//- Return the matrix source coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type> > gradientBoundaryCoeffs
(
const scalarField& deltaCoeffs
) const;
//- Return the matrix source coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type> > gradientBoundaryCoeffs() const;
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<Type>& matrix);
// Cyclic AMI coupled interface functions
//- Does the patch field perform the transformation
virtual bool doTransform() const
{
return
!(cyclicACMIPatch_.parallel() || pTraits<Type>::rank == 0);
}
//- Return face transformation tensor
virtual const tensorField& forwardT() const
{
return cyclicACMIPatch_.forwardT();
}
//- Return neighbour-cell transformation tensor
virtual const tensorField& reverseT() const
{
return cyclicACMIPatch_.reverseT();
}
//- Return rank of component for transform
virtual int rank() const
{
return pTraits<Type>::rank;
}
// I-O
//- Write
virtual void write(Ostream& os) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "cyclicACMIFvPatchField.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,43 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cyclicACMIFvPatchFields.H"
#include "addToRunTimeSelectionTable.H"
#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
makePatchFields(cyclicACMI);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,49 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#ifndef cyclicACMIFvPatchFields_H
#define cyclicACMIFvPatchFields_H
#include "cyclicACMIFvPatchField.H"
#include "fieldTypes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeFieldTypedefs(cyclicACMI);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,50 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#ifndef cyclicACMIFvPatchFieldsFwd_H
#define cyclicACMIFvPatchFieldsFwd_H
#include "fieldTypes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type> class cyclicACMIFvPatchField;
makePatchTypeFieldTypedefs(cyclicACMI);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -42,6 +42,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(p),
internalField_(iF),
updated_(false),
manipulatedMatrix_(false),
patchType_(word::null)
{}
@ -58,6 +59,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(p),
internalField_(iF),
updated_(false),
manipulatedMatrix_(false),
patchType_(word::null)
{}
@ -75,6 +77,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(p),
internalField_(iF),
updated_(false),
manipulatedMatrix_(false),
patchType_(ptf.patchType_)
{}
@ -92,6 +95,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(p),
internalField_(iF),
updated_(false),
manipulatedMatrix_(false),
patchType_(dict.lookupOrDefault<word>("patchType", word::null))
{
if (dict.found("value"))
@ -133,6 +137,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(ptf.patch_),
internalField_(ptf.internalField_),
updated_(false),
manipulatedMatrix_(false),
patchType_(ptf.patchType_)
{}
@ -148,6 +153,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(ptf.patch_),
internalField_(iF),
updated_(false),
manipulatedMatrix_(false),
patchType_(ptf.patchType_)
{}
@ -267,6 +273,28 @@ void Foam::fvPatchField<Type>::rmap
}
template<class Type>
void Foam::fvPatchField<Type>::updateCoeffs()
{
updated_ = true;
}
template<class Type>
void Foam::fvPatchField<Type>::updateCoeffs(const scalarField& weights)
{
if (!updated_)
{
updateCoeffs();
Field<Type>& fld = *this;
fld *= weights;
updated_ = true;
}
}
template<class Type>
void Foam::fvPatchField<Type>::evaluate(const Pstream::commsTypes)
{
@ -276,13 +304,25 @@ void Foam::fvPatchField<Type>::evaluate(const Pstream::commsTypes)
}
updated_ = false;
manipulatedMatrix_ = false;
}
template<class Type>
void Foam::fvPatchField<Type>::manipulateMatrix(fvMatrix<Type>& matrix)
{
// do nothing
manipulatedMatrix_ = true;
}
template<class Type>
void Foam::fvPatchField<Type>::manipulateMatrix
(
fvMatrix<Type>& matrix,
const scalarField& weights
)
{
manipulatedMatrix_ = true;
}

View File

@ -93,6 +93,10 @@ class fvPatchField
// the construction of the matrix
bool updated_;
//- Update index used so that manipulateMatrix is called only once
// during the construction of the matrix
bool manipulatedMatrix_;
//- Optional patch type, used to allow specified boundary conditions
// to be applied to constraint patches by providing the constraint
// patch type as 'patchType'
@ -327,6 +331,12 @@ public:
return updated_;
}
//- Return true if the matrix has already been manipulated
bool manipulatedMatrix() const
{
return manipulatedMatrix_;
}
// Mapping functions
@ -365,10 +375,12 @@ public:
//- Update the coefficients associated with the patch field
// Sets Updated to true
virtual void updateCoeffs()
{
updated_ = true;
}
virtual void updateCoeffs();
//- Update the coefficients associated with the patch field
// and apply weight field
// Sets Updated to true
virtual void updateCoeffs(const scalarField& weights);
//- Return internal field next to patch as patch field
virtual tmp<Field<Type> > patchInternalField() const;
@ -479,6 +491,13 @@ public:
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<Type>& matrix);
//- Manipulate matrix with given weights
virtual void manipulateMatrix
(
fvMatrix<Type>& matrix,
const scalarField& weights
);
// I-O

View File

@ -0,0 +1,149 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cyclicACMIFvsPatchField.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
Foam::cyclicACMIFvsPatchField<Type>::cyclicACMIFvsPatchField
(
const fvPatch& p,
const DimensionedField<Type, surfaceMesh>& iF
)
:
coupledFvsPatchField<Type>(p, iF),
cyclicACMIPatch_(refCast<const cyclicACMIFvPatch>(p))
{}
template<class Type>
Foam::cyclicACMIFvsPatchField<Type>::cyclicACMIFvsPatchField
(
const cyclicACMIFvsPatchField<Type>& ptf,
const fvPatch& p,
const DimensionedField<Type, surfaceMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
coupledFvsPatchField<Type>(ptf, p, iF, mapper),
cyclicACMIPatch_(refCast<const cyclicACMIFvPatch>(p))
{
if (!isA<cyclicACMIFvPatch>(this->patch()))
{
FatalErrorIn
(
"cyclicACMIFvsPatchField<Type>::cyclicACMIFvsPatchField\n"
"("
"const cyclicACMIFvsPatchField<Type>&, "
"const fvPatch&, "
"const DimensionedField<Type, surfaceMesh>&, "
"const fvPatchFieldMapper&"
")"
) << "Field type does not correspond to patch type for patch "
<< this->patch().index() << "." << endl
<< "Field type: " << typeName << endl
<< "Patch type: " << this->patch().type()
<< exit(FatalError);
}
}
template<class Type>
Foam::cyclicACMIFvsPatchField<Type>::cyclicACMIFvsPatchField
(
const fvPatch& p,
const DimensionedField<Type, surfaceMesh>& iF,
const dictionary& dict
)
:
coupledFvsPatchField<Type>(p, iF, dict),
cyclicACMIPatch_(refCast<const cyclicACMIFvPatch>(p))
{
if (!isA<cyclicACMIFvPatch>(p))
{
FatalIOErrorIn
(
"cyclicACMIFvsPatchField<Type>::cyclicACMIFvsPatchField"
"("
"const fvPatch&, "
"const Field<Type>&, "
"const dictionary&"
")",
dict
) << "patch " << this->patch().index() << " not cyclicACMI type. "
<< "Patch type = " << p.type()
<< exit(FatalIOError);
}
}
template<class Type>
Foam::cyclicACMIFvsPatchField<Type>::cyclicACMIFvsPatchField
(
const cyclicACMIFvsPatchField<Type>& ptf
)
:
coupledFvsPatchField<Type>(ptf),
cyclicACMIPatch_(ptf.cyclicACMIPatch_)
{}
template<class Type>
Foam::cyclicACMIFvsPatchField<Type>::cyclicACMIFvsPatchField
(
const cyclicACMIFvsPatchField<Type>& ptf,
const DimensionedField<Type, surfaceMesh>& iF
)
:
coupledFvsPatchField<Type>(ptf, iF),
cyclicACMIPatch_(ptf.cyclicACMIPatch_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
bool Foam::cyclicACMIFvsPatchField<Type>::coupled() const
{
if
(
Pstream::parRun()
|| (
this->cyclicACMIPatch_.size()
&& this->cyclicACMIPatch_.cyclicACMIPatch().neighbPatch().size()
)
)
{
return true;
}
else
{
return false;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,150 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::cyclicACMIFvsPatchField
Description
Foam::cyclicACMIFvsPatchField
SourceFiles
cyclicACMIFvsPatchField.C
\*---------------------------------------------------------------------------*/
#ifndef cyclicACMIFvsPatchField_H
#define cyclicACMIFvsPatchField_H
#include "coupledFvsPatchField.H"
#include "cyclicACMIFvPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class cyclicACMIFvsPatchField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class cyclicACMIFvsPatchField
:
public coupledFvsPatchField<Type>
{
// Private data
//- Local reference cast into the cyclic patch
const cyclicACMIFvPatch& cyclicACMIPatch_;
public:
//- Runtime type information
TypeName(cyclicACMIFvPatch::typeName_());
// Constructors
//- Construct from patch and internal field
cyclicACMIFvsPatchField
(
const fvPatch&,
const DimensionedField<Type, surfaceMesh>&
);
//- Construct from patch, internal field and dictionary
cyclicACMIFvsPatchField
(
const fvPatch&,
const DimensionedField<Type, surfaceMesh>&,
const dictionary&
);
//- Construct by mapping given cyclicACMIFvsPatchField onto a new patch
cyclicACMIFvsPatchField
(
const cyclicACMIFvsPatchField<Type>&,
const fvPatch&,
const DimensionedField<Type, surfaceMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
cyclicACMIFvsPatchField
(
const cyclicACMIFvsPatchField<Type>&
);
//- Construct and return a clone
virtual tmp<fvsPatchField<Type> > clone() const
{
return tmp<fvsPatchField<Type> >
(
new cyclicACMIFvsPatchField<Type>(*this)
);
}
//- Construct as copy setting internal field reference
cyclicACMIFvsPatchField
(
const cyclicACMIFvsPatchField<Type>&,
const DimensionedField<Type, surfaceMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvsPatchField<Type> > clone
(
const DimensionedField<Type, surfaceMesh>& iF
) const
{
return tmp<fvsPatchField<Type> >
(
new cyclicACMIFvsPatchField<Type>(*this, iF)
);
}
// Member functions
// Access
//- Return true if running parallel
virtual bool coupled() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "cyclicACMIFvsPatchField.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,43 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cyclicACMIFvsPatchFields.H"
#include "fvsPatchFields.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
makeFvsPatchFields(cyclicACMI);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,49 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#ifndef cyclicACMIFvsPatchFields_H
#define cyclicACMIFvsPatchFields_H
#include "cyclicACMIFvsPatchField.H"
#include "fieldTypes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makeFvsPatchTypeFieldTypedefs(cyclicACMI);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,50 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#ifndef cyclicACMIFvsPatchFieldsFwd_H
#define cyclicACMIFvsPatchFieldsFwd_H
#include "fieldTypes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type> class cyclicACMIFvsPatchField;
makeFvsPatchTypeFieldTypedefs(cyclicACMI);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -235,7 +235,7 @@ Foam::fv::faceMDLimitedGrad<Foam::vector>::calcGrad
g[own],
maxFace - vvfOwn,
minFace - vvfOwn,
Cf[facei] - C[nei]
Cf[facei] - C[own]
);

View File

@ -0,0 +1,202 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cyclicACMIFvPatch.H"
#include "addToRunTimeSelectionTable.H"
#include "fvMesh.H"
#include "transform.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(cyclicACMIFvPatch, 0);
addToRunTimeSelectionTable(fvPatch, cyclicACMIFvPatch, polyPatch);
}
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void Foam::cyclicACMIFvPatch::updateAreas() const
{
if (cyclicACMIPolyPatch_.updated())
{
// Set Sf and magSf for both sides' coupled and non-overlapping patches
// owner couple
const_cast<vectorField&>(Sf()) = patch().faceAreas();
const_cast<scalarField&>(magSf()) = mag(patch().faceAreas());
// owner non-overlapping
const fvPatch& nonOverlapPatch = this->nonOverlapPatch();
const_cast<vectorField&>(nonOverlapPatch.Sf()) =
nonOverlapPatch.patch().faceAreas();
const_cast<scalarField&>(nonOverlapPatch.magSf()) =
mag(nonOverlapPatch.patch().faceAreas());
// neighbour couple
const cyclicACMIFvPatch& nbrACMI = neighbPatch();
const_cast<vectorField&>(nbrACMI.Sf()) =
nbrACMI.patch().faceAreas();
const_cast<scalarField&>(nbrACMI.magSf()) =
mag(nbrACMI.patch().faceAreas());
// neighbour non-overlapping
const fvPatch& nbrNonOverlapPatch = nbrACMI.nonOverlapPatch();
const_cast<vectorField&>(nbrNonOverlapPatch.Sf()) =
nbrNonOverlapPatch.patch().faceAreas();
const_cast<scalarField&>(nbrNonOverlapPatch.magSf()) =
mag(nbrNonOverlapPatch.patch().faceAreas());
// set the updated flag
cyclicACMIPolyPatch_.setUpdated(false);
}
}
void Foam::cyclicACMIFvPatch::makeWeights(scalarField& w) const
{
if (coupled())
{
const cyclicACMIFvPatch& nbrPatch = neighbFvPatch();
const fvPatch& nbrPatchNonOverlap = nonOverlapPatch();
const scalarField deltas(nf() & fvPatch::delta());
const scalarField nbrDeltas
(
interpolate
(
nbrPatch.nf() & nbrPatch.fvPatch::delta(),
nbrPatchNonOverlap.nf() & nbrPatchNonOverlap.delta()
)
);
forAll(deltas, faceI)
{
scalar di = deltas[faceI];
scalar dni = nbrDeltas[faceI];
w[faceI] = dni/(di + dni);
}
}
else
{
// Behave as uncoupled patch
fvPatch::makeWeights(w);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::cyclicACMIFvPatch::coupled() const
{
return Pstream::parRun() || (this->size() && neighbFvPatch().size());
}
Foam::tmp<Foam::vectorField> Foam::cyclicACMIFvPatch::delta() const
{
if (coupled())
{
const cyclicACMIFvPatch& nbrPatchCoupled = neighbFvPatch();
const fvPatch& nbrPatchNonOverlap = nonOverlapPatch();
const vectorField patchD(fvPatch::delta());
vectorField nbrPatchD
(
interpolate
(
nbrPatchCoupled.fvPatch::delta(),
nbrPatchNonOverlap.delta()
)
);
const vectorField nbrPatchD0
(
interpolate
(
vectorField(nbrPatchCoupled.size(), vector::zero),
nbrPatchNonOverlap.delta()()
)
);
nbrPatchD -= nbrPatchD0;
tmp<vectorField> tpdv(new vectorField(patchD.size()));
vectorField& pdv = tpdv();
// do the transformation if necessary
if (parallel())
{
forAll(patchD, faceI)
{
const vector& ddi = patchD[faceI];
const vector& dni = nbrPatchD[faceI];
pdv[faceI] = ddi - dni;
}
}
else
{
forAll(patchD, faceI)
{
const vector& ddi = patchD[faceI];
const vector& dni = nbrPatchD[faceI];
pdv[faceI] = ddi - transform(forwardT()[0], dni);
}
}
return tpdv;
}
else
{
return fvPatch::delta();
}
}
Foam::tmp<Foam::labelField> Foam::cyclicACMIFvPatch::interfaceInternalField
(
const labelUList& internalData
) const
{
return patchInternalField(internalData);
}
Foam::tmp<Foam::labelField> Foam::cyclicACMIFvPatch::internalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF
) const
{
return neighbFvPatch().patchInternalField(iF);
}
// ************************************************************************* //

View File

@ -0,0 +1,268 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::cyclicACMIFvPatch
Description
Cyclic patch for Arbitrarily Coupled Mesh Interface (ACMI)
SourceFiles
cyclicACMIFvPatch.C
\*---------------------------------------------------------------------------*/
#ifndef cyclicACMIFvPatch_H
#define cyclicACMIFvPatch_H
#include "coupledFvPatch.H"
#include "cyclicACMILduInterface.H"
#include "cyclicACMIPolyPatch.H"
#include "fvBoundaryMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class cyclicACMIFvPatch Declaration
\*---------------------------------------------------------------------------*/
class cyclicACMIFvPatch
:
public coupledFvPatch,
public cyclicACMILduInterface
{
// Private data
const cyclicACMIPolyPatch& cyclicACMIPolyPatch_;
protected:
// Protected Member functions
//- Update the patch areas after AMI update
void updateAreas() const;
//- Make patch weighting factors
void makeWeights(scalarField&) const;
public:
//- Runtime type information
TypeName(cyclicACMIPolyPatch::typeName_());
// Constructors
//- Construct from polyPatch
cyclicACMIFvPatch(const polyPatch& patch, const fvBoundaryMesh& bm)
:
coupledFvPatch(patch, bm),
cyclicACMILduInterface(),
cyclicACMIPolyPatch_(refCast<const cyclicACMIPolyPatch>(patch))
{}
// Member functions
// Access
//- Return local reference cast into the cyclic patch
const cyclicACMIPolyPatch& cyclicACMIPatch() const
{
return cyclicACMIPolyPatch_;
}
//- Return neighbour
virtual label neighbPatchID() const
{
return cyclicACMIPolyPatch_.neighbPatchID();
}
virtual bool owner() const
{
return cyclicACMIPolyPatch_.owner();
}
//- Return neighbour fvPatch
virtual const cyclicACMIFvPatch& neighbPatch() const
{
return refCast<const cyclicACMIFvPatch>
(
this->boundaryMesh()[cyclicACMIPolyPatch_.neighbPatchID()]
);
}
//- Return neighbour
virtual label nonOverlapPatchID() const
{
return cyclicACMIPolyPatch_.nonOverlapPatchID();
}
//- Return non-overlapping fvPatch
virtual const fvPatch& nonOverlapPatch() const
{
return this->boundaryMesh()[nonOverlapPatchID()];
}
//- Return a reference to the AMI interpolator
virtual const AMIPatchToPatchInterpolation& AMI() const
{
const AMIPatchToPatchInterpolation& AMI =
cyclicACMIPolyPatch_.AMI();
updateAreas();
return AMI;
}
//- Are the cyclic planes parallel
virtual bool parallel() const
{
return cyclicACMIPolyPatch_.parallel();
}
//- Return face transformation tensor
virtual const tensorField& forwardT() const
{
return cyclicACMIPolyPatch_.forwardT();
}
//- Return neighbour-cell transformation tensor
virtual const tensorField& reverseT() const
{
return cyclicACMIPolyPatch_.reverseT();
}
const cyclicACMIFvPatch& neighbFvPatch() const
{
return refCast<const cyclicACMIFvPatch>
(
this->boundaryMesh()[cyclicACMIPolyPatch_.neighbPatchID()]
);
}
//- Return true if this patch is coupled. This is equivalent
// to the coupledPolyPatch::coupled() if parallel running or
// both sides present, false otherwise
virtual bool coupled() const;
//- Return delta (P to N) vectors across coupled patch
virtual tmp<vectorField> delta() const;
template<class Type>
tmp<Field<Type> > interpolate
(
const Field<Type>& fldCoupled
) const
{
updateAreas();
return
cyclicACMIPolyPatch_.cyclicAMIPolyPatch::interpolate
(
fldCoupled
);
}
template<class Type>
tmp<Field<Type> > interpolate
(
const tmp<Field<Type> >& tfldCoupled
) const
{
updateAreas();
return
cyclicACMIPolyPatch_.cyclicAMIPolyPatch::interpolate
(
tfldCoupled
);
}
template<class Type>
tmp<Field<Type> > interpolate
(
const Field<Type>& fldCoupled,
const Field<Type>& fldNonOverlap
) const
{
updateAreas();
return
cyclicACMIPolyPatch_.interpolate
(
fldCoupled,
fldNonOverlap
);
}
template<class Type>
tmp<Field<Type> > interpolate
(
const tmp<Field<Type> >& tFldCoupled,
const tmp<Field<Type> >& tFldNonOverlap
) const
{
updateAreas();
return
cyclicACMIPolyPatch_.interpolate
(
tFldCoupled,
tFldNonOverlap
);
}
// Interface transfer functions
//- Return the values of the given internal data adjacent to
// the interface as a field
virtual tmp<labelField> interfaceInternalField
(
const labelUList& internalData
) const;
//- Return neighbour field
virtual tmp<labelField> internalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& internalData
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -111,16 +111,7 @@ public:
mesh.gradScheme(gradSchemeName_)
)
)
{
if (!schemeData.eof())
{
IOWarningIn("linearUpwind(const fvMesh&, Istream&)", schemeData)
<< "unexpected additional entries in stream." << nl
<< " Only the name of the gradient scheme in the"
" 'gradSchemes' dictionary should be specified."
<< endl;
}
}
{}
//- Construct from faceFlux and Istream
linearUpwind
@ -140,21 +131,7 @@ public:
mesh.gradScheme(gradSchemeName_)
)
)
{
if (!schemeData.eof())
{
IOWarningIn
(
"linearUpwind(const fvMesh&, "
"const surfaceScalarField& faceFlux, Istream&)",
schemeData
) << "unexpected additional entries in stream." << nl
<< " Only the name of the gradient scheme in the"
" 'gradSchemes' dictionary should be specified."
<< endl;
}
}
{}
// Member Functions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -110,19 +110,7 @@ public:
mesh.gradScheme(gradSchemeName_)
)
)
{
if (!schemeData.eof())
{
IOWarningIn
(
"linearUpwindV(const fvMesh&, Istream&)",
schemeData
) << "unexpected additional entries in stream." << nl
<< " Only the name of the gradient scheme in the"
" 'gradSchemes' dictionary should be specified."
<< endl;
}
}
{}
//- Construct from faceFlux and Istream
linearUpwindV
@ -142,20 +130,7 @@ public:
mesh.gradScheme(gradSchemeName_)
)
)
{
if (!schemeData.eof())
{
IOWarningIn
(
"linearUpwindV(const fvMesh&, "
"const surfaceScalarField& faceFlux, Istream&)",
schemeData
) << "unexpected additional entries in stream." << nl
<< " Only the name of the gradient scheme in the"
" 'gradSchemes' dictionary should be specified."
<< endl;
}
}
{}
// Member Functions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -156,6 +156,26 @@ Foam::tmp<Foam::labelField> Foam::MGridGenGAMGAgglomeration::agglomerate
finalAgglom.begin()
);
{
label nNewCoarseCells = 0;
labelList newRestrictAddr;
bool ok = checkRestriction
(
newRestrictAddr,
nNewCoarseCells
,
fineAddressing,
finalAgglom,
nCoarseCells
);
if (!ok)
{
nCoarseCells = nNewCoarseCells;
finalAgglom.transfer(newRestrictAddr);
}
}
return tmp<labelField>(new labelField(finalAgglom));
}

View File

@ -566,7 +566,11 @@ void Foam::KinematicCloud<CloudType>::checkParcelProperties
{
const scalar carrierDt = mesh_.time().deltaTValue();
parcel.stepFraction() = (carrierDt - lagrangianDt)/carrierDt;
parcel.typeId() = constProps_.parcelTypeId();
if (parcel.typeId() == -1)
{
parcel.typeId() = constProps_.parcelTypeId();
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -95,7 +95,9 @@ Foam::forceSuSp Foam::SRFForce<CloudType>::calcNonCoupled
const vector& r = p.position();
// Coriolis and centrifugal acceleration terms
value.Su() = mass*(2.0*(p.U() ^ omega) + (omega ^ (r ^ omega)));
value.Su() =
mass*(1.0 - p.rhoc()/p.rho())
*(2.0*(p.U() ^ omega) + (omega ^ (r ^ omega)));
return value;
}

View File

@ -336,7 +336,7 @@ void Foam::autoLayerDriver::smoothNormals
) const
{
// Get smoothly varying internal normals field.
Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
Info<< "shrinkMeshDistance : Smoothing normals in interior ..." << endl;
const fvMesh& mesh = meshRefiner_.mesh();
const edgeList& edges = mesh.edges();
@ -1161,6 +1161,11 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
scalar mDist = medialDist[pointI];
if (wDist2 < sqr(SMALL) && mDist < SMALL)
//- Note: maybe less strict:
//(
// wDist2 < sqr(meshRefiner_.mergeDistance())
// && mDist < meshRefiner_.mergeDistance()
//)
{
medialRatio[pointI] = 0.0;
}

View File

@ -35,6 +35,7 @@ License
#include "shellSurfaces.H"
#include "mapDistributePolyMesh.H"
#include "unitConversion.H"
#include "snapParameters.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -610,9 +611,16 @@ Foam::label Foam::autoRefineDriver::danglingCellRefine
<< " cells (out of " << mesh.globalData().nTotalCells()
<< ')' << endl;
// Stop when no cells to refine. No checking of minRefineCells since
// too few cells
if (nCellsToRefine == 0)
// Stop when no cells to refine. After a few iterations check if too
// few cells
if
(
nCellsToRefine == 0
|| (
iter >= 1
&& nCellsToRefine <= refineParams.minRefineCells()
)
)
{
Info<< "Stopping refining since too few cells selected."
<< nl << endl;
@ -870,6 +878,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
void Foam::autoRefineDriver::baffleAndSplitMesh
(
const refinementParameters& refineParams,
const snapParameters& snapParams,
const bool handleSnapProblems,
const dictionary& motionDict
)
@ -887,10 +896,17 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
meshRefiner_.baffleAndSplitMesh
(
handleSnapProblems, // detect&remove potential snap problem
// Snap problem cell detection
snapParams,
refineParams.useTopologicalSnapDetection(),
false, // perpendicular edge connected cells
scalarField(0), // per region perpendicular angle
// Free standing baffles
!handleSnapProblems, // merge free standing baffles?
refineParams.planarAngle(),
motionDict,
const_cast<Time&>(mesh.time()),
globalToMasterPatch_,
@ -951,6 +967,7 @@ void Foam::autoRefineDriver::zonify
void Foam::autoRefineDriver::splitAndMergeBaffles
(
const refinementParameters& refineParams,
const snapParameters& snapParams,
const bool handleSnapProblems,
const dictionary& motionDict
)
@ -973,10 +990,17 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
meshRefiner_.baffleAndSplitMesh
(
handleSnapProblems,
// Snap problem cell detection
snapParams,
refineParams.useTopologicalSnapDetection(),
handleSnapProblems, // remove perp edge connected cells
perpAngle, // perp angle
// Free standing baffles
true, // merge free standing baffles?
refineParams.planarAngle(), // planar angle
motionDict,
const_cast<Time&>(mesh.time()),
globalToMasterPatch_,
@ -1081,6 +1105,7 @@ void Foam::autoRefineDriver::doRefine
(
const dictionary& refineDict,
const refinementParameters& refineParams,
const snapParameters& snapParams,
const bool prepareForSnapping,
const dictionary& motionDict
)
@ -1146,13 +1171,25 @@ void Foam::autoRefineDriver::doRefine
// Introduce baffles at surface intersections. Remove sections unreachable
// from keepPoint.
baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);
baffleAndSplitMesh
(
refineParams,
snapParams,
prepareForSnapping,
motionDict
);
// Mesh is at its finest. Do optional zoning.
zonify(refineParams);
// Pull baffles apart
splitAndMergeBaffles(refineParams, prepareForSnapping, motionDict);
splitAndMergeBaffles
(
refineParams,
snapParams,
prepareForSnapping,
motionDict
);
// Do something about cells with refined faces on the boundary
if (prepareForSnapping)

View File

@ -43,6 +43,8 @@ namespace Foam
// Forward declaration of classes
class refinementParameters;
class snapParameters;
class meshRefinement;
class decompositionMethod;
class fvMeshDistribute;
@ -121,6 +123,7 @@ class autoRefineDriver
void baffleAndSplitMesh
(
const refinementParameters& refineParams,
const snapParameters& snapParams,
const bool handleSnapProblems,
const dictionary& motionDict
);
@ -131,6 +134,7 @@ class autoRefineDriver
void splitAndMergeBaffles
(
const refinementParameters& refineParams,
const snapParameters& snapParams,
const bool handleSnapProblems,
const dictionary& motionDict
);
@ -176,6 +180,7 @@ public:
(
const dictionary& refineDict,
const refinementParameters& refineParams,
const snapParameters& snapParams,
const bool prepareForSnapping,
const dictionary& motionDict
);

View File

@ -121,7 +121,7 @@ Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
(
const motionSmoother& meshMover,
const List<labelPair>& baffles
) const
)
{
const indirectPrimitivePatch& pp = meshMover.patch();
@ -615,14 +615,14 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles
Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
(
const fvMesh& mesh,
const snapParameters& snapParams,
const indirectPrimitivePatch& pp
) const
)
{
const edgeList& edges = pp.edges();
const labelListList& pointEdges = pp.pointEdges();
const pointField& localPoints = pp.localPoints();
const fvMesh& mesh = meshRefiner_.mesh();
scalarField maxEdgeLen(localPoints.size(), -GREAT);
@ -655,13 +655,14 @@ Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
void Foam::autoSnapDriver::preSmoothPatch
(
const meshRefinement& meshRefiner,
const snapParameters& snapParams,
const label nInitErrors,
const List<labelPair>& baffles,
motionSmoother& meshMover
) const
)
{
const fvMesh& mesh = meshRefiner_.mesh();
const fvMesh& mesh = meshRefiner.mesh();
labelList checkFaces;
@ -724,11 +725,11 @@ void Foam::autoSnapDriver::preSmoothPatch
{
const_cast<Time&>(mesh.time())++;
Info<< "Writing patch smoothed mesh to time "
<< meshRefiner_.timeName() << '.' << endl;
meshRefiner_.write
<< meshRefiner.timeName() << '.' << endl;
meshRefiner.write
(
debug,
mesh.time().path()/meshRefiner_.timeName()
mesh.time().path()/meshRefiner.timeName()
);
Info<< "Dumped mesh in = "
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
@ -742,12 +743,11 @@ void Foam::autoSnapDriver::preSmoothPatch
// Get (pp-local) indices of points that are both on zone and on patched surface
Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
(
const fvMesh& mesh,
const indirectPrimitivePatch& pp,
const word& zoneName
) const
)
{
const fvMesh& mesh = meshRefiner_.mesh();
label zoneI = mesh.faceZones().findZoneID(zoneName);
if (zoneI == -1)
@ -755,7 +755,7 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
FatalErrorIn
(
"autoSnapDriver::getZoneSurfacePoints"
"(const indirectPrimitivePatch&, const word&)"
"(const fvMesh&, const indirectPrimitivePatch&, const word&)"
) << "Cannot find zone " << zoneName
<< exit(FatalError);
}
@ -793,17 +793,17 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
(
const meshRefinement& meshRefiner,
const scalarField& snapDist,
motionSmoother& meshMover
) const
const indirectPrimitivePatch& pp
)
{
Info<< "Calculating patchDisplacement as distance to nearest surface"
<< " point ..." << endl;
const indirectPrimitivePatch& pp = meshMover.patch();
const pointField& localPoints = pp.localPoints();
const refinementSurfaces& surfaces = meshRefiner_.surfaces();
const fvMesh& mesh = meshRefiner_.mesh();
const refinementSurfaces& surfaces = meshRefiner.surfaces();
const fvMesh& mesh = meshRefiner.mesh();
// Displacement per patch point
vectorField patchDisp(localPoints.size(), vector::zero);
@ -815,9 +815,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
// Divide surfaces into zoned and unzoned
labelList zonedSurfaces =
meshRefiner_.surfaces().getNamedSurfaces();
meshRefiner.surfaces().getNamedSurfaces();
labelList unzonedSurfaces =
meshRefiner_.surfaces().getUnnamedSurfaces();
meshRefiner.surfaces().getUnnamedSurfaces();
// 1. All points to non-interface surfaces
@ -870,6 +870,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
(
getZoneSurfacePoints
(
mesh,
pp,
faceZoneNames[zoneSurfI]
)
@ -966,6 +967,254 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
}
////XXXXXXXXX
//// Get (pp-local) indices of points that are on both patches
//Foam::labelList Foam::autoSnapDriver::getPatchSurfacePoints
//(
// const fvMesh& mesh,
// const indirectPrimitivePatch& allPp,
// const polyPatch& pp
//)
//{
// // Could use PrimitivePatch & localFaces to extract points but might just
// // as well do it ourselves.
//
// boolList pointOnZone(allPp.nPoints(), false);
//
// forAll(pp, i)
// {
// const face& f = pp[i];
//
// forAll(f, fp)
// {
// label meshPointI = f[fp];
//
// Map<label>::const_iterator iter =
// allPp.meshPointMap().find(meshPointI);
//
// if (iter != allPp.meshPointMap().end())
// {
// label pointI = iter();
// pointOnZone[pointI] = true;
// }
// }
// }
//
// return findIndices(pointOnZone, true);
//}
//Foam::vectorField Foam::autoSnapDriver::calcNearestLocalSurface
//(
// const meshRefinement& meshRefiner,
// const scalarField& snapDist,
// const indirectPrimitivePatch& pp
//)
//{
// Info<< "Calculating patchDisplacement as distance to nearest"
// << " local surface point ..." << endl;
//
// const pointField& localPoints = pp.localPoints();
// const refinementSurfaces& surfaces = meshRefiner.surfaces();
// const fvMesh& mesh = meshRefiner.mesh();
// const polyBoundaryMesh& pbm = mesh.boundaryMesh();
//
//
//// // Assume that all patch-internal points get attracted to their surface
//// // only. So we want to know if point is on multiple regions
////
//// labelList minPatch(mesh.nPoints(), labelMax);
//// labelList maxPatch(mesh.nPoints(), labelMin);
////
//// forAll(meshMover.adaptPatchIDs(), i)
//// {
//// label patchI = meshMover.adaptPatchIDs()[i];
//// const labelList& meshPoints = pbm[patchI].meshPoints();
////
//// forAll(meshPoints, meshPointI)
//// {
//// label meshPointI = meshPoints[meshPointI];
//// minPatch[meshPointI] = min(minPatch[meshPointI], patchI);
//// maxPatch[meshPointI] = max(maxPatch[meshPointI], patchI);
//// }
//// }
////
//// syncTools::syncPointList
//// (
//// mesh,
//// minPatch,
//// minEqOp<label>(), // combine op
//// labelMax // null value
//// );
//// syncTools::syncPointList
//// (
//// mesh,
//// maxPatch,
//// maxEqOp<label>(), // combine op
//// labelMin // null value
//// );
//
// // Now all points with minPatch != maxPatch will be on the outside of
// // the patch.
//
// // Displacement per patch point
// vectorField patchDisp(localPoints.size(), vector::zero);
// // Current best snap distance
// scalarField minSnapDist(snapDist);
// // Current surface snapped to
// labelList snapSurf(localPoints.size(), -1);
//
// const labelList& surfaceGeometry = surfaces.surfaces();
// forAll(surfaceGeometry, surfI)
// {
// label geomI = surfaceGeometry[surfI];
// const wordList& regNames = allGeometry.regionNames()[geomI];
// forAll(regNames, regionI)
// {
// label globalRegionI = surfaces.globalRegion(surfI, regionI);
// // Collect master patch points
// label masterPatchI = globalToMasterPatch_[globalRegionI];
// label slavePatchI = globalToSlavePatch_[globalRegionI];
//
// labelList patchPointIndices
// (
// getPatchSurfacePoints
// (
// mesh,
// pp,
// pbm[masterPatchI]
// )
// );
//
// // Find nearest for points both on faceZone and pp.
// List<pointIndexHit> hitInfo;
// surfaces.findNearest
// (
// surfI,
// regionI,
// pointField(localPoints, patchPointIndices),
// sqr(scalarField(minSnapDist, patchPointIndices)),
// hitInfo
// );
//
// forAll(hitInfo, i)
// {
// label pointI = patchPointIndices[i];
//
// if (hitInfo[i].hit())
// {
// const point& pt = hitInfo[i].hitPoint();
// patchDisp[pointI] = pt-localPoints[pointI];
// minSnapDist[pointI] = min
// (
// minSnapDist[pointI],
// mag(patchDisp[pointI])
// );
// snapSurf[pointI] = surfI;
// }
// }
//
// // Slave patch
// if (slavePatchI != masterPatchI)
// {
// labelList patchPointIndices
// (
// getPatchSurfacePoints
// (
// mesh,
// pp,
// pbm[slavePatchI]
// )
// );
//
// // Find nearest for points both on faceZone and pp.
// List<pointIndexHit> hitInfo;
// surfaces.findNearest
// (
// surfI,
// regionI,
// pointField(localPoints, patchPointIndices),
// sqr(scalarField(minSnapDist, patchPointIndices)),
// hitInfo
// );
//
// forAll(hitInfo, i)
// {
// label pointI = patchPointIndices[i];
//
// if (hitInfo[i].hit())
// {
// const point& pt = hitInfo[i].hitPoint();
// patchDisp[pointI] = pt-localPoints[pointI];
// minSnapDist[pointI] = min
// (
// minSnapDist[pointI],
// mag(patchDisp[pointI])
// );
// snapSurf[pointI] = surfI;
// }
// }
// }
// }
// }
//
//
// // Check if all points are being snapped
// forAll(snapSurf, pointI)
// {
// if (snapSurf[pointI] == -1)
// {
// WarningIn("autoSnapDriver::calcNearestLocalSurface(..)")
// << "For point:" << pointI
// << " coordinate:" << localPoints[pointI]
// << " did not find any surface within:"
// << minSnapDist[pointI]
// << " metre." << endl;
// }
// }
//
// {
// scalarField magDisp(mag(patchDisp));
//
// Info<< "Wanted displacement : average:"
// << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
// << " min:" << gMin(magDisp)
// << " max:" << gMax(magDisp) << endl;
// }
//
//
// Info<< "Calculated surface displacement in = "
// << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
//
//
// // Limit amount of movement.
// forAll(patchDisp, patchPointI)
// {
// scalar magDisp = mag(patchDisp[patchPointI]);
//
// if (magDisp > snapDist[patchPointI])
// {
// patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
//
// Pout<< "Limiting displacement for " << patchPointI
// << " from " << magDisp << " to " << snapDist[patchPointI]
// << endl;
// }
// }
//
// // Points on zones in one domain but only present as point on other
// // will not do condition 2 on all. Sync explicitly.
// syncTools::syncPointList
// (
// mesh,
// pp.meshPoints(),
// patchDisp,
// minMagSqrEqOp<point>(), // combine op
// vector(GREAT, GREAT, GREAT)// null value (note: cannot use VGREAT)
// );
//
// return patchDisp;
//}
////XXXXXXXXX
void Foam::autoSnapDriver::smoothDisplacement
(
const snapParameters& snapParams,
@ -1153,7 +1402,15 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
scalarField faceSnapDist(pp.size(), -GREAT);
{
// Distance to attract to nearest feature on surface
const scalarField snapDist(calcSnapDistance(snapParams, pp));
const scalarField snapDist
(
calcSnapDistance
(
mesh,
snapParams,
pp
)
);
const faceList& localFaces = pp.localFaces();
@ -1429,7 +1686,7 @@ void Foam::autoSnapDriver::doSnap
indirectPrimitivePatch& pp = ppPtr();
// Distance to attract to nearest feature on surface
const scalarField snapDist(calcSnapDistance(snapParams, pp));
const scalarField snapDist(calcSnapDistance(mesh, snapParams, pp));
// Construct iterative mesh mover.
@ -1467,7 +1724,14 @@ void Foam::autoSnapDriver::doSnap
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
// Pre-smooth patch vertices (so before determining nearest)
preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
preSmoothPatch
(
meshRefiner_,
snapParams,
nInitErrors,
baffles,
meshMover
);
for (label iter = 0; iter < nFeatIter; iter++)
@ -1478,7 +1742,12 @@ void Foam::autoSnapDriver::doSnap
// Calculate displacement at every patch point. Insert into
// meshMover.
vectorField disp = calcNearestSurface(snapDist, meshMover);
vectorField disp = calcNearestSurface
(
meshRefiner_,
snapDist,
meshMover.patch()
);
// Override displacement with feature edge attempt
if (doFeatures)

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -81,11 +81,11 @@ class autoSnapDriver
//- Calculate displacement per patch point to smooth out patch.
// Quite complicated in determining which points to move where.
pointField smoothPatchDisplacement
static pointField smoothPatchDisplacement
(
const motionSmoother&,
const List<labelPair>&
) const;
);
//- Check that face zones are synced
void checkCoupledFaceZones() const;
@ -406,37 +406,51 @@ public:
autoPtr<mapPolyMesh> mergeZoneBaffles(const List<labelPair>&);
//- Calculate edge length per patch point.
scalarField calcSnapDistance
static scalarField calcSnapDistance
(
const fvMesh& mesh,
const snapParameters& snapParams,
const indirectPrimitivePatch&
) const;
);
//- Smooth the mesh (patch and internal) to increase visibility
// of surface points (on castellated mesh) w.r.t. surface.
void preSmoothPatch
static void preSmoothPatch
(
const meshRefinement& meshRefiner,
const snapParameters& snapParams,
const label nInitErrors,
const List<labelPair>& baffles,
motionSmoother&
) const;
);
//- Get points both on patch and facezone.
labelList getZoneSurfacePoints
static labelList getZoneSurfacePoints
(
const fvMesh& mesh,
const indirectPrimitivePatch&,
const word& zoneName
) const;
);
//- Per patch point calculate point on nearest surface. Set as
// boundary conditions of motionSmoother displacement field. Return
// displacement of patch points.
vectorField calcNearestSurface
static vectorField calcNearestSurface
(
const meshRefinement& meshRefiner,
const scalarField& snapDist,
motionSmoother& meshMover
) const;
const indirectPrimitivePatch&
);
////- Per patch point calculate point on nearest surface. Set as
//// boundary conditions of motionSmoother displacement field.
//// Return displacement of patch points.
//static vectorField calcNearestLocalSurface
//(
// const meshRefinement& meshRefiner,
// const scalarField& snapDist,
// const indirectPrimitivePatch&
//);
//- Smooth the displacement field to the internal.
void smoothDisplacement

View File

@ -45,7 +45,11 @@ Foam::refinementParameters::refinementParameters
nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))),
keepPoints_(dict.lookup("keepPoints")),
allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
useTopologicalSnapDetection_
(
dict.lookupOrDefault<bool>("useTopologicalSnapDetection", true)
),
maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance", 0))
{}
@ -65,7 +69,11 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
keepPoints_(pointField(1, dict.lookup("locationInMesh"))),
allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
useTopologicalSnapDetection_
(
dict.lookupOrDefault<bool>("useTopologicalSnapDetection", true)
),
maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance", 0))
{
scalar featAngle(readScalar(dict.lookup("resolveFeatureAngle")));
@ -83,7 +91,7 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh)
const
const
{
// Force calculation of tet-diag decomposition (for use in findCell)
(void)mesh.tetBasePtIs();

View File

@ -80,6 +80,10 @@ class refinementParameters
// cellZone?
Switch allowFreeStandingZoneFaces_;
//- Use old topology based problem-cell removal (cells with 8 points
// on surface)
Switch useTopologicalSnapDetection_;
//- Allowed load unbalance
scalar maxLoadUnbalance_;
@ -157,6 +161,13 @@ public:
return allowFreeStandingZoneFaces_;
}
//- Use old topology based problem-cell removal
// (cells with 8 points on surface)
bool useTopologicalSnapDetection() const
{
return useTopologicalSnapDetection_;
}
//- Allowed load unbalance
scalar maxLoadUnbalance() const
{

View File

@ -72,6 +72,8 @@ class globalIndex;
class removePoints;
class localPointRegion;
class snapParameters;
/*---------------------------------------------------------------------------*\
Class meshRefinement Declaration
\*---------------------------------------------------------------------------*/
@ -470,6 +472,14 @@ private:
const labelList& globalToMasterPatch
) const;
//- Returns list with for every internal face -1 or the patch
// they should be baffled into.
labelList markFacesOnProblemCellsGeometric
(
const snapParameters& snapParams,
const dictionary& motionDict
) const;
// Baffle merging
@ -536,6 +546,8 @@ private:
//- Remove any loose standing cells
void handleSnapProblems
(
const snapParameters& snapParams,
const bool useTopologicalSnapDetection,
const bool removeEdgeConnectedCells,
const scalarField& perpendicularAngle,
const dictionary& motionDict,
@ -760,10 +772,17 @@ public:
void baffleAndSplitMesh
(
const bool handleSnapProblems,
// How to remove problem snaps
const snapParameters& snapParams,
const bool useTopologicalSnapDetection,
const bool removeEdgeConnectedCells,
const scalarField& perpendicularAngle,
// How to handle free-standing baffles
const bool mergeFreeStanding,
const scalar freeStandingAngle,
const dictionary& motionDict,
Time& runTime,
const labelList& globalToMasterPatch,

View File

@ -1872,6 +1872,8 @@ void Foam::meshRefinement::makeConsistentFaceIndex
void Foam::meshRefinement::handleSnapProblems
(
const snapParameters& snapParams,
const bool useTopologicalSnapDetection,
const bool removeEdgeConnectedCells,
const scalarField& perpendicularAngle,
const dictionary& motionDict,
@ -1885,44 +1887,39 @@ void Foam::meshRefinement::handleSnapProblems
<< "----------------------------------------------" << nl
<< endl;
labelList facePatch
(
markFacesOnProblemCells
labelList facePatch;
if (useTopologicalSnapDetection)
{
facePatch = markFacesOnProblemCells
(
motionDict,
removeEdgeConnectedCells,
perpendicularAngle,
globalToMasterPatch
)
);
);
}
else
{
facePatch = markFacesOnProblemCellsGeometric(snapParams, motionDict);
}
Info<< "Analyzed problem cells in = "
<< runTime.cpuTimeIncrement() << " s\n" << nl << endl;
if (debug&meshRefinement::MESH)
{
faceSet problemTopo(mesh_, "problemFacesTopo", 100);
faceSet problemFaces(mesh_, "problemFaces", 100);
const labelList facePatchTopo
(
markFacesOnProblemCells
(
motionDict,
removeEdgeConnectedCells,
perpendicularAngle,
globalToMasterPatch
)
);
forAll(facePatchTopo, faceI)
forAll(facePatch, faceI)
{
if (facePatchTopo[faceI] != -1)
if (facePatch[faceI] != -1)
{
problemTopo.insert(faceI);
problemFaces.insert(faceI);
}
}
problemTopo.instance() = timeName();
Pout<< "Dumping " << problemTopo.size()
<< " problem faces to " << problemTopo.objectPath() << endl;
problemTopo.write();
problemFaces.instance() = timeName();
Pout<< "Dumping " << problemFaces.size()
<< " problem faces to " << problemFaces.objectPath() << endl;
problemFaces.write();
}
Info<< "Introducing baffles to delete problem cells." << nl << endl;
@ -1962,6 +1959,8 @@ void Foam::meshRefinement::handleSnapProblems
void Foam::meshRefinement::baffleAndSplitMesh
(
const bool doHandleSnapProblems,
const snapParameters& snapParams,
const bool useTopologicalSnapDetection,
const bool removeEdgeConnectedCells,
const scalarField& perpendicularAngle,
const bool mergeFreeStanding,
@ -2029,83 +2028,10 @@ void Foam::meshRefinement::baffleAndSplitMesh
if (doHandleSnapProblems)
{
//Info<< nl
// << "Introducing baffles to block off problem cells" << nl
// << "----------------------------------------------" << nl
// << endl;
//
//labelList facePatch
//(
// markFacesOnProblemCells
// (
// motionDict,
// removeEdgeConnectedCells,
// perpendicularAngle,
// globalToMasterPatch
// )
// //markFacesOnProblemCellsGeometric(motionDict)
//);
//Info<< "Analyzed problem cells in = "
// << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
//
//if (debug&meshRefinement::MESH)
//{
// faceSet problemTopo(mesh_, "problemFacesTopo", 100);
//
// const labelList facePatchTopo
// (
// markFacesOnProblemCells
// (
// motionDict,
// removeEdgeConnectedCells,
// perpendicularAngle,
// globalToMasterPatch
// )
// );
// forAll(facePatchTopo, faceI)
// {
// if (facePatchTopo[faceI] != -1)
// {
// problemTopo.insert(faceI);
// }
// }
// problemTopo.instance() = timeName();
// Pout<< "Dumping " << problemTopo.size()
// << " problem faces to " << problemTopo.objectPath() << endl;
// problemTopo.write();
//}
//
//Info<< "Introducing baffles to delete problem cells." << nl << endl;
//
//if (debug)
//{
// runTime++;
//}
//
//// Create baffles with same owner and neighbour for now.
//createBaffles(facePatch, facePatch);
//
//if (debug)
//{
// // Debug:test all is still synced across proc patches
// checkData();
//}
//Info<< "Created baffles in = "
// << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
//
//printMeshInfo(debug, "After introducing baffles");
//
//if (debug&meshRefinement::MESH)
//{
// Pout<< "Writing extra baffled mesh to time "
// << timeName() << endl;
// write(debug, runTime.path()/"extraBaffles");
// Pout<< "Dumped debug data in = "
// << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
//}
handleSnapProblems
(
snapParams,
useTopologicalSnapDetection,
removeEdgeConnectedCells,
perpendicularAngle,
motionDict,
@ -2192,6 +2118,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
// and delete them
handleSnapProblems
(
snapParams,
useTopologicalSnapDetection,
removeEdgeConnectedCells,
perpendicularAngle,
motionDict,

View File

@ -36,6 +36,10 @@ License
#include "polyMeshGeometry.H"
#include "IOmanip.H"
#include "unitConversion.H"
#include "autoSnapDriver.H"
#include "snapParameters.H"
#include "motionSmoother.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -968,164 +972,241 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
}
//// Mark faces to be baffled to prevent snapping problems. Does
//// test to find nearest surface and checks which faces would get squashed.
//Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
//(
// const dictionary& motionDict
//) const
//{
// // Construct addressing engine.
// autoPtr<indirectPrimitivePatch> ppPtr
// (
// meshRefinement::makePatch
// (
// mesh_,
// meshedPatches()
// )
// );
// const indirectPrimitivePatch& pp = ppPtr();
// const pointField& localPoints = pp.localPoints();
// const labelList& meshPoints = pp.meshPoints();
//
// // Find nearest (non-baffle) surface
// pointField newPoints(mesh_.points());
// {
// List<pointIndexHit> hitInfo;
// labelList hitSurface;
// surfaces_.findNearest
// (
// surfaces_.getUnnamedSurfaces(),
// localPoints,
// scalarField(localPoints.size(), sqr(GREAT)),// sqr of attraction
// hitSurface,
// hitInfo
// );
//
// forAll(hitInfo, i)
// {
// if (hitInfo[i].hit())
// {
// //label pointI = meshPoints[i];
// //Pout<< " " << pointI << " moved from "
// // << mesh_.points()[pointI] << " by "
// // << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
// // << endl;
// newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
// }
// }
// }
//
// // Per face (internal or coupled!) the patch that the
// // baffle should get (or -1).
// labelList facePatch(mesh_.nFaces(), -1);
// // Count of baffled faces
// label nBaffleFaces = 0;
//
// {
// pointField oldPoints(mesh_.points());
// mesh_.movePoints(newPoints);
// faceSet wrongFaces(mesh_, "wrongFaces", 100);
// {
// //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
//
// // Just check the errors from squashing
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// const labelList allFaces(identity(mesh_.nFaces()));
// label nWrongFaces = 0;
//
// scalar minArea(readScalar(motionDict.lookup("minArea")));
// if (minArea > -SMALL)
// {
// polyMeshGeometry::checkFaceArea
// (
// false,
// minArea,
// mesh_,
// mesh_.faceAreas(),
// allFaces,
// &wrongFaces
// );
//
// label nNewWrongFaces = returnReduce
// (
// wrongFaces.size(),
// sumOp<label>()
// );
//
// Info<< " faces with area < "
// << setw(5) << minArea
// << " m^2 : "
// << nNewWrongFaces-nWrongFaces << endl;
//
// nWrongFaces = nNewWrongFaces;
// }
//
//// scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
// scalar minDet = 0.01;
// if (minDet > -1)
// {
// polyMeshGeometry::checkCellDeterminant
// (
// false,
// minDet,
// mesh_,
// mesh_.faceAreas(),
// allFaces,
// polyMeshGeometry::affectedCells(mesh_, allFaces),
// &wrongFaces
// );
//
// label nNewWrongFaces = returnReduce
// (
// wrongFaces.size(),
// sumOp<label>()
// );
//
// Info<< " faces on cells with determinant < "
// << setw(5) << minDet << " : "
// << nNewWrongFaces-nWrongFaces << endl;
//
// nWrongFaces = nNewWrongFaces;
// }
// }
//
//
// forAllConstIter(faceSet, wrongFaces, iter)
// {
// label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
//
// if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
// {
// facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
// nBaffleFaces++;
//
// //Pout<< " " << iter.key()
// // //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
// // << " is destined for patch " << facePatch[iter.key()]
// // << endl;
// }
// }
// // Restore points.
// mesh_.movePoints(oldPoints);
// }
//
//
// Info<< "markFacesOnProblemCellsGeometric : marked "
// << returnReduce(nBaffleFaces, sumOp<label>())
// << " additional internal and coupled faces"
// << " to be converted into baffles." << endl;
//
// syncTools::syncFaceList
// (
// mesh_,
// facePatch,
// maxEqOp<label>()
// );
//
// return facePatch;
//}
// Mark faces to be baffled to prevent snapping problems. Does
// test to find nearest surface and checks which faces would get squashed.
Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
(
const snapParameters& snapParams,
const dictionary& motionDict
) const
{
pointField oldPoints(mesh_.points());
// Repeat (most of) autoSnapDriver::doSnap
{
labelList adaptPatchIDs(meshedPatches());
// Construct addressing engine.
autoPtr<indirectPrimitivePatch> ppPtr
(
meshRefinement::makePatch
(
mesh_,
adaptPatchIDs
)
);
indirectPrimitivePatch& pp = ppPtr();
// Distance to attract to nearest feature on surface
const scalarField snapDist
(
autoSnapDriver::calcSnapDistance(mesh_, snapParams, pp)
);
// Construct iterative mesh mover.
Info<< "Constructing mesh displacer ..." << endl;
Info<< "Using mesh parameters " << motionDict << nl << endl;
const pointMesh& pMesh = pointMesh::New(mesh_);
motionSmoother meshMover
(
mesh_,
pp,
adaptPatchIDs,
meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs)(),
motionDict
);
// Check initial mesh
Info<< "Checking initial mesh ..." << endl;
labelHashSet wrongFaces(mesh_.nFaces()/100);
motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
const label nInitErrors = returnReduce
(
wrongFaces.size(),
sumOp<label>()
);
Info<< "Detected " << nInitErrors << " illegal faces"
<< " (concave, zero area or negative cell pyramid volume)"
<< endl;
Info<< "Checked initial mesh in = "
<< mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
// Pre-smooth patch vertices (so before determining nearest)
autoSnapDriver::preSmoothPatch
(
*this,
snapParams,
nInitErrors,
List<labelPair>(0), //baffles
meshMover
);
const vectorField disp
(
autoSnapDriver::calcNearestSurface
(
*this,
snapDist, // attraction
pp
)
);
const labelList& meshPoints = pp.meshPoints();
pointField newPoints(mesh_.points());
forAll(meshPoints, i)
{
newPoints[meshPoints[i]] += disp[i];
}
mesh_.movePoints(newPoints);
}
// Per face (internal or coupled!) the patch that the
// baffle should get (or -1).
labelList facePatch(mesh_.nFaces(), -1);
// Count of baffled faces
label nBaffleFaces = 0;
{
faceSet wrongFaces(mesh_, "wrongFaces", 100);
{
//motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
// Just check the errors from squashing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const labelList allFaces(identity(mesh_.nFaces()));
label nWrongFaces = 0;
//const scalar minV(readScalar(motionDict.lookup("minVol", true)));
//if (minV > -GREAT)
//{
// polyMeshGeometry::checkFacePyramids
// (
// false,
// minV,
// mesh_,
// mesh_.cellCentres(),
// mesh_.points(),
// allFaces,
// List<labelPair>(0),
// &wrongFaces
// );
//
// label nNewWrongFaces = returnReduce
// (
// wrongFaces.size(),
// sumOp<label>()
// );
//
// Info<< " faces with pyramid volume < "
// << setw(5) << minV
// << " m^3 : "
// << nNewWrongFaces-nWrongFaces << endl;
//
// nWrongFaces = nNewWrongFaces;
//}
scalar minArea(readScalar(motionDict.lookup("minArea")));
if (minArea > -SMALL)
{
polyMeshGeometry::checkFaceArea
(
false,
minArea,
mesh_,
mesh_.faceAreas(),
allFaces,
&wrongFaces
);
label nNewWrongFaces = returnReduce
(
wrongFaces.size(),
sumOp<label>()
);
Info<< " faces with area < "
<< setw(5) << minArea
<< " m^2 : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
if (minDet > -1)
{
polyMeshGeometry::checkCellDeterminant
(
false,
minDet,
mesh_,
mesh_.faceAreas(),
allFaces,
polyMeshGeometry::affectedCells(mesh_, allFaces),
&wrongFaces
);
label nNewWrongFaces = returnReduce
(
wrongFaces.size(),
sumOp<label>()
);
Info<< " faces on cells with determinant < "
<< setw(5) << minDet << " : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
}
forAllConstIter(faceSet, wrongFaces, iter)
{
label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
{
facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
nBaffleFaces++;
//Pout<< " " << iter.key()
// //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
// << " is destined for patch " << facePatch[iter.key()]
// << endl;
}
}
}
// Restore points.
mesh_.movePoints(oldPoints);
Info<< "markFacesOnProblemCellsGeometric : marked "
<< returnReduce(nBaffleFaces, sumOp<label>())
<< " additional internal and coupled faces"
<< " to be converted into baffles." << endl;
syncTools::syncFaceList
(
mesh_,
facePatch,
maxEqOp<label>()
);
return facePatch;
}
// ************************************************************************* //

View File

@ -77,7 +77,8 @@ const Foam::NamedEnum<Foam::refinementSurfaces::faceZoneType, 3>
Foam::refinementSurfaces::refinementSurfaces
(
const searchableSurfaces& allGeometry,
const dictionary& surfacesDict
const dictionary& surfacesDict,
const label gapLevelIncrement
)
:
allGeometry_(allGeometry),
@ -143,7 +144,7 @@ Foam::refinementSurfaces::refinementSurfaces
globalLevelIncr[surfI] = dict.lookupOrDefault
(
"gapLevelIncrement",
0
gapLevelIncrement
);
if
@ -274,7 +275,7 @@ Foam::refinementSurfaces::refinementSurfaces
label levelIncr = regionDict.lookupOrDefault
(
"gapLevelIncrement",
0
gapLevelIncrement
);
regionLevelIncr[surfI].insert(regionI, levelIncr);
@ -402,6 +403,49 @@ Foam::refinementSurfaces::refinementSurfaces
}
Foam::refinementSurfaces::refinementSurfaces
(
const searchableSurfaces& allGeometry,
const labelList& surfaces,
const wordList& names,
const wordList& faceZoneNames,
const wordList& cellZoneNames,
const List<areaSelectionAlgo>& zoneInside,
const pointField& zoneInsidePoints,
const List<faceZoneType>& faceType,
const labelList& regionOffset,
const labelList& minLevel,
const labelList& maxLevel,
const labelList& gapLevel,
const scalarField& perpendicularAngle,
const PtrList<dictionary>& patchInfo
)
:
allGeometry_(allGeometry),
surfaces_(surfaces),
names_(names),
faceZoneNames_(faceZoneNames),
cellZoneNames_(cellZoneNames),
zoneInside_(zoneInside),
zoneInsidePoints_(zoneInsidePoints),
faceType_(faceType),
regionOffset_(regionOffset),
minLevel_(minLevel),
maxLevel_(maxLevel),
gapLevel_(gapLevel),
perpendicularAngle_(perpendicularAngle),
patchInfo_(patchInfo.size())
{
forAll(patchInfo_, pI)
{
if (patchInfo.set(pI))
{
patchInfo_[pI] = patchInfo[pI];
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Get indices of unnamed surfaces (surfaces without faceZoneName)

View File

@ -110,7 +110,7 @@ private:
pointField zoneInsidePoints_;
//- Per 'interface' surface :
// Waht to do with outside
// What to do with outside
List<faceZoneType> faceType_;
//- From local region number to global region number
@ -149,7 +149,27 @@ public:
refinementSurfaces
(
const searchableSurfaces& allGeometry,
const dictionary&
const dictionary&,
const label gapLevelIncrement
);
//- Construct from components
refinementSurfaces
(
const searchableSurfaces& allGeometry,
const labelList& surfaces,
const wordList& names,
const wordList& faceZoneNames,
const wordList& cellZoneNames,
const List<areaSelectionAlgo>& zoneInside,
const pointField& zoneInsidePoints,
const List<faceZoneType>& faceType,
const labelList& regionOffset,
const labelList& minLevel,
const labelList& maxLevel,
const labelList& gapLevel,
const scalarField& perpendicularAngle,
const PtrList<dictionary>& patchInfo
);

View File

@ -56,6 +56,11 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolationMethodToWord
method = "faceAreaWeightAMI";
break;
}
case imPartialFaceAreaWeight:
{
method = "partialFaceAreaWeightAMI";
break;
}
default:
{
FatalErrorIn
@ -87,7 +92,15 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::wordTointerpolationMethod
wordList methods
(
IStringStream("(directAMI mapNearestAMI faceAreaWeightAMI)")()
IStringStream
(
"("
"directAMI "
"mapNearestAMI "
"faceAreaWeightAMI "
"partialFaceAreaWeightAMI"
")"
)()
);
if (im == "directAMI")
@ -102,6 +115,10 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::wordTointerpolationMethod
{
method = imFaceAreaWeight;
}
else if (im == "partialFaceAreaWeightAMI")
{
method = imPartialFaceAreaWeight;
}
else
{
FatalErrorIn
@ -183,6 +200,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
const labelListList& addr,
scalarListList& wght,
scalarField& wghtSum,
const bool conformal,
const bool output
)
{
@ -191,12 +209,19 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
forAll(wght, faceI)
{
scalarList& w = wght[faceI];
scalar denom = patchAreas[faceI];
scalar s = sum(w);
scalar t = s/patchAreas[faceI];
scalar t = s/denom;
if (conformal)
{
denom = s;
}
forAll(w, i)
{
w[i] /= s;
w[i] /= denom;
}
wghtSum[faceI] = t;
@ -476,6 +501,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
srcAddress,
srcWeights,
srcWeightsSum,
true,
false
);
}
@ -749,7 +775,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
tgtMagSf_[faceI] = tgtPatch[faceI].mag(tgtPatch.points());
}
// Calculate if patches present on multiple processors
singlePatchProc_ = calcDistribution(srcPatch, tgtPatch);
@ -794,32 +819,28 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
newTgtPoints
);
// calculate AMI interpolation
{
autoPtr<AMIMethod<SourcePatch, TargetPatch> > AMIPtr
autoPtr<AMIMethod<SourcePatch, TargetPatch> > AMIPtr
(
AMIMethod<SourcePatch, TargetPatch>::New
(
AMIMethod<SourcePatch, TargetPatch>::New
(
interpolationMethodToWord(method_),
srcPatch,
newTgtPatch,
srcMagSf_,
tgtMagSf_,
triMode_,
reverseTarget_
)
);
AMIPtr->calculate
(
srcAddress_,
srcWeights_,
tgtAddress_,
tgtWeights_
);
}
interpolationMethodToWord(method_),
srcPatch,
newTgtPatch,
srcMagSf_,
tgtMagSf_,
triMode_,
reverseTarget_
)
);
AMIPtr->calculate
(
srcAddress_,
srcWeights_,
tgtAddress_,
tgtWeights_
);
// Now
// ~~~
@ -886,6 +907,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
srcAddress_,
srcWeights_,
srcWeightsSum_,
AMIPtr->conformal(),
true
);
normaliseWeights
@ -895,6 +917,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
tgtAddress_,
tgtWeights_,
tgtWeightsSum_,
AMIPtr->conformal(),
true
);
@ -903,7 +926,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
srcMapPtr_.reset(new mapDistribute(globalSrcFaces, tgtAddress_, cMap));
tgtMapPtr_.reset(new mapDistribute(globalTgtFaces, srcAddress_, cMap));
if (debug)
{
writeFaceConnectivity(srcPatch, newTgtPatch, srcAddress_);
@ -913,29 +935,27 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
{
// calculate AMI interpolation
{
autoPtr<AMIMethod<SourcePatch, TargetPatch> > AMIPtr
autoPtr<AMIMethod<SourcePatch, TargetPatch> > AMIPtr
(
AMIMethod<SourcePatch, TargetPatch>::New
(
AMIMethod<SourcePatch, TargetPatch>::New
(
interpolationMethodToWord(method_),
srcPatch,
tgtPatch,
srcMagSf_,
tgtMagSf_,
triMode_,
reverseTarget_
)
);
interpolationMethodToWord(method_),
srcPatch,
tgtPatch,
srcMagSf_,
tgtMagSf_,
triMode_,
reverseTarget_
)
);
AMIPtr->calculate
(
srcAddress_,
srcWeights_,
tgtAddress_,
tgtWeights_
);
}
AMIPtr->calculate
(
srcAddress_,
srcWeights_,
tgtAddress_,
tgtWeights_
);
normaliseWeights
(
@ -944,6 +964,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
srcAddress_,
srcWeights_,
srcWeightsSum_,
AMIPtr->conformal(),
true
);
normaliseWeights
@ -953,6 +974,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
tgtAddress_,
tgtWeights_,
tgtWeightsSum_,
AMIPtr->conformal(),
true
);
}

View File

@ -88,7 +88,8 @@ public:
{
imDirect,
imMapNearest,
imFaceAreaWeight
imFaceAreaWeight,
imPartialFaceAreaWeight
};
//- Convert interpolationMethod to word representation
@ -236,6 +237,7 @@ private:
const labelListList& addr,
scalarListList& wght,
scalarField& wghtSum,
const bool conformal,
const bool output
);

View File

@ -41,24 +41,28 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::checkPatches() const
}
const scalar maxBoundsError = 0.05;
// check bounds of source and target
boundBox bbSrc(srcPatch_.points(), srcPatch_.meshPoints(), true);
boundBox bbTgt(tgtPatch_.points(), tgtPatch_.meshPoints(), true);
boundBox bbTgtInf(bbTgt);
bbTgtInf.inflate(maxBoundsError);
if (!bbTgtInf.contains(bbSrc))
if (conformal())
{
WarningIn("AMIMethod<SourcePatch, TargetPatch>::checkPatches()")
<< "Source and target patch bounding boxes are not similar" << nl
<< " source box span : " << bbSrc.span() << nl
<< " target box span : " << bbTgt.span() << nl
<< " source box : " << bbSrc << nl
<< " target box : " << bbTgt << nl
<< " inflated target box : " << bbTgtInf << endl;
const scalar maxBoundsError = 0.05;
// check bounds of source and target
boundBox bbSrc(srcPatch_.points(), srcPatch_.meshPoints(), true);
boundBox bbTgt(tgtPatch_.points(), tgtPatch_.meshPoints(), true);
boundBox bbTgtInf(bbTgt);
bbTgtInf.inflate(maxBoundsError);
if (!bbTgtInf.contains(bbSrc))
{
WarningIn("AMIMethod<SourcePatch, TargetPatch>::checkPatches()")
<< "Source and target patch bounding boxes are not similar"
<< nl
<< " source box span : " << bbSrc.span() << nl
<< " target box span : " << bbTgt.span() << nl
<< " source box : " << bbSrc << nl
<< " target box : " << bbTgt << nl
<< " inflated target box : " << bbTgtInf << endl;
}
}
}
@ -74,6 +78,8 @@ bool Foam::AMIMethod<SourcePatch, TargetPatch>::initialise
label& tgtFaceI
)
{
checkPatches();
// set initial sizes for weights and addressing - must be done even if
// returns false below
srcAddress.setSize(srcPatch_.size());
@ -241,17 +247,16 @@ Foam::label Foam::AMIMethod<SourcePatch, TargetPatch>::findTargetFace
pointIndexHit sample = treePtr_->findNearest(srcPt, 10.0*srcFaceArea);
if (debug)
{
Pout<< "Source point = " << srcPt << ", Sample point = "
<< sample.hitPoint() << ", Sample index = " << sample.index()
<< endl;
}
if (sample.hit())
{
targetFaceI = sample.index();
if (debug)
{
Pout<< "Source point = " << srcPt << ", Sample point = "
<< sample.hitPoint() << ", Sample index = " << sample.index()
<< endl;
}
}
return targetFaceI;
@ -334,8 +339,6 @@ Foam::AMIMethod<SourcePatch, TargetPatch>::AMIMethod
srcNonOverlap_(),
triMode_(triMode)
{
checkPatches();
label srcSize = returnReduce(srcPatch_.size(), sumOp<label>());
label tgtSize = returnReduce(tgtPatch_.size(), sumOp<label>());
@ -352,4 +355,13 @@ Foam::AMIMethod<SourcePatch, TargetPatch>::~AMIMethod()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
bool Foam::AMIMethod<SourcePatch, TargetPatch>::conformal() const
{
return true;
}
// ************************************************************************* //

View File

@ -207,6 +207,9 @@ public:
// Note: this should be empty for correct functioning
inline const labelList& srcNonOverlap() const;
//- Flag to indicate that interpolation patches are conformal
virtual bool conformal() const;
// Manipulation

View File

@ -25,7 +25,85 @@ License
#include "faceAreaWeightAMI.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
(
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght,
label srcFaceI,
label tgtFaceI
)
{
// construct weights and addressing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nFacesRemaining = srcAddr.size();
// list of tgt face neighbour faces
DynamicList<label> nbrFaces(10);
// list of faces currently visited for srcFaceI to avoid multiple hits
DynamicList<label> visitedFaces(10);
// list to keep track of tgt faces used to seed src faces
labelList seedFaces(nFacesRemaining, -1);
seedFaces[srcFaceI] = tgtFaceI;
// list to keep track of whether src face can be mapped
boolList mapFlag(nFacesRemaining, true);
// reset starting seed
label startSeedI = 0;
DynamicList<label> nonOverlapFaces;
do
{
// Do advancing front starting from srcFaceI,tgtFaceI
bool faceProcessed = processSourceFace
(
srcFaceI,
tgtFaceI,
nbrFaces,
visitedFaces,
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
mapFlag[srcFaceI] = false;
nFacesRemaining--;
if (!faceProcessed)
{
nonOverlapFaces.append(srcFaceI);
}
// choose new src face from current src face neighbour
if (nFacesRemaining > 0)
{
setNextFaces
(
startSeedI,
srcFaceI,
tgtFaceI,
mapFlag,
seedFaces,
visitedFaces
);
}
} while (nFacesRemaining > 0);
this->srcNonOverlap_.transfer(nonOverlapFaces);
}
template<class SourcePatch, class TargetPatch>
bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
@ -45,6 +123,11 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
List<DynamicList<scalar> >& tgtWght
)
{
if (tgtStartFaceI == -1)
{
return false;
}
nbrFaces.clear();
visitedFaces.clear();
@ -101,11 +184,15 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
label& tgtFaceI,
const boolList& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces
const DynamicList<label>& visitedFaces,
bool errorOnNotFound
) const
{
const labelList& srcNbrFaces = this->srcPatch_.faceFaces()[srcFaceI];
// initialise tgtFaceI
tgtFaceI = -1;
// set possible seeds for later use
bool valuesSet = false;
forAll(srcNbrFaces, i)
@ -195,19 +282,23 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
}
}
FatalErrorIn
(
"void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::"
"setNextFaces"
"("
"label&, "
"label&, "
"label&, "
"const boolList&, "
"labelList&, "
"const DynamicList<label>&"
") const"
) << "Unable to set source and target faces" << abort(FatalError);
if (errorOnNotFound)
{
FatalErrorIn
(
"void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::"
"setNextFaces"
"("
"label&, "
"label&, "
"label&, "
"const boolList&, "
"labelList&, "
"const DynamicList<label>&, "
"bool"
") const"
) << "Unable to set source and target faces" << abort(FatalError);
}
}
}
@ -447,77 +538,21 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
List<DynamicList<label> > tgtAddr(this->tgtPatch_.size());
List<DynamicList<scalar> > tgtWght(tgtAddr.size());
calcAddressing
(
srcAddr,
srcWght,
tgtAddr,
tgtWght,
srcFaceI,
tgtFaceI
);
// construct weights and addressing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nFacesRemaining = srcAddr.size();
// list of tgt face neighbour faces
DynamicList<label> nbrFaces(10);
// list of faces currently visited for srcFaceI to avoid multiple hits
DynamicList<label> visitedFaces(10);
// list to keep track of tgt faces used to seed src faces
labelList seedFaces(nFacesRemaining, -1);
seedFaces[srcFaceI] = tgtFaceI;
// list to keep track of whether src face can be mapped
boolList mapFlag(nFacesRemaining, true);
// reset starting seed
label startSeedI = 0;
DynamicList<label> nonOverlapFaces;
do
if (this->srcNonOverlap_.size() != 0)
{
// Do advancing front starting from srcFaceI,tgtFaceI
bool faceProcessed = processSourceFace
(
srcFaceI,
tgtFaceI,
nbrFaces,
visitedFaces,
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
mapFlag[srcFaceI] = false;
nFacesRemaining--;
if (!faceProcessed)
{
nonOverlapFaces.append(srcFaceI);
}
// choose new src face from current src face neighbour
if (nFacesRemaining > 0)
{
setNextFaces
(
startSeedI,
srcFaceI,
tgtFaceI,
mapFlag,
seedFaces,
visitedFaces
);
}
} while (nFacesRemaining > 0);
if (nonOverlapFaces.size() != 0)
{
Pout<< " AMI: " << nonOverlapFaces.size()
Pout<< " AMI: " << this->srcNonOverlap_.size()
<< " non-overlap faces identified"
<< endl;
this->srcNonOverlap_.transfer(nonOverlapFaces);
}

View File

@ -52,9 +52,9 @@ class faceAreaWeightAMI
public AMIMethod<SourcePatch, TargetPatch>
{
private:
protected:
// Private Member Functions
// Protected Member Functions
//- Disallow default bitwise copy construct
faceAreaWeightAMI(const faceAreaWeightAMI&);
@ -64,8 +64,19 @@ private:
// Marching front
//- Calculate addressing and weights using temporary storage
virtual void calcAddressing
(
List<DynamicList<label> >& srcAddress,
List<DynamicList<scalar> >& srcWeights,
List<DynamicList<label> >& tgtAddress,
List<DynamicList<scalar> >& tgtWeights,
label srcFaceI,
label tgtFaceI
);
//- Determine overlap contributions for source face srcFaceI
bool processSourceFace
virtual bool processSourceFace
(
const label srcFaceI,
const label tgtStartFaceI,
@ -78,7 +89,7 @@ private:
);
//- Attempt to re-evaluate source faces that have not been included
void restartUncoveredSourceFace
virtual void restartUncoveredSourceFace
(
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
@ -87,21 +98,22 @@ private:
);
//- Set the source and target seed faces
void setNextFaces
virtual void setNextFaces
(
label& startSeedI,
label& srcFaceI,
label& tgtFaceI,
const boolList& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces
const DynamicList<label>& visitedFaces,
bool errorOnNotFound = true
) const;
// Evaluation
//- Area of intersection between source and target faces
scalar interArea
virtual scalar interArea
(
const label srcFaceI,
const label tgtFaceI

View File

@ -0,0 +1,155 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "partialFaceAreaWeightAMI.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
(
label& startSeedI,
label& srcFaceI,
label& tgtFaceI,
const boolList& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces,
const bool errorOnNotFound
) const
{
faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
(
startSeedI,
srcFaceI,
tgtFaceI,
mapFlag,
seedFaces,
visitedFaces,
false // no error on not found
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::
partialFaceAreaWeightAMI
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget
)
:
faceAreaWeightAMI<SourcePatch, TargetPatch>
(
srcPatch,
tgtPatch,
srcMagSf,
tgtMagSf,
triMode,
reverseTarget
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::
~partialFaceAreaWeightAMI()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
bool Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::conformal() const
{
return false;
}
template<class SourcePatch, class TargetPatch>
void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFaceI,
label tgtFaceI
)
{
bool ok =
this->initialise
(
srcAddress,
srcWeights,
tgtAddress,
tgtWeights,
srcFaceI,
tgtFaceI
);
if (!ok)
{
return;
}
// temporary storage for addressing and weights
List<DynamicList<label> > srcAddr(this->srcPatch_.size());
List<DynamicList<scalar> > srcWght(srcAddr.size());
List<DynamicList<label> > tgtAddr(this->tgtPatch_.size());
List<DynamicList<scalar> > tgtWght(tgtAddr.size());
faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
(
srcAddr,
srcWght,
tgtAddr,
tgtWght,
srcFaceI,
tgtFaceI
);
// transfer data to persistent storage
forAll(srcAddr, i)
{
srcAddress[i].transfer(srcAddr[i]);
srcWeights[i].transfer(srcWght[i]);
}
forAll(tgtAddr, i)
{
tgtAddress[i].transfer(tgtAddr[i]);
tgtWeights[i].transfer(tgtWght[i]);
}
}
// ************************************************************************* //

View File

@ -0,0 +1,141 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::partialFaceAreaWeightAMI
Description
Partial face area weighted Arbitrary Mesh Interface (AMI) method
SourceFiles
partialFaceAreaWeightAMI.C
\*---------------------------------------------------------------------------*/
#ifndef partialFaceAreaWeightAMI_H
#define partialFaceAreaWeightAMI_H
#include "faceAreaWeightAMI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class partialFaceAreaWeightAMI Declaration
\*---------------------------------------------------------------------------*/
template<class SourcePatch, class TargetPatch>
class partialFaceAreaWeightAMI
:
public faceAreaWeightAMI<SourcePatch, TargetPatch>
{
private:
// Private Member Functions
//- Disallow default bitwise copy construct
partialFaceAreaWeightAMI(const partialFaceAreaWeightAMI&);
//- Disallow default bitwise assignment
void operator=(const partialFaceAreaWeightAMI&);
// Marching front
//- Set the source and target seed faces
virtual void setNextFaces
(
label& startSeedI,
label& srcFaceI,
label& tgtFaceI,
const boolList& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces,
bool errorOnNotFound = true
) const;
public:
//- Runtime type information
TypeName("partialFaceAreaWeightAMI");
// Constructors
//- Construct from components
partialFaceAreaWeightAMI
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget = false
);
//- Destructor
virtual ~partialFaceAreaWeightAMI();
// Member Functions
// Access
//- Flag to indicate that interpolation patches are conformal
virtual bool conformal() const;
// Manipulation
//- Update addressing and weights
virtual void calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFaceI = -1,
label tgtFaceI = -1
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "partialFaceAreaWeightAMI.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -28,6 +28,7 @@ License
#include "directAMI.H"
#include "mapNearestAMI.H"
#include "faceAreaWeightAMI.H"
#include "partialFaceAreaWeightAMI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -38,6 +39,7 @@ namespace Foam
makeAMIMethodType(AMIPatchToPatchInterpolation, directAMI);
makeAMIMethodType(AMIPatchToPatchInterpolation, mapNearestAMI);
makeAMIMethodType(AMIPatchToPatchInterpolation, faceAreaWeightAMI);
makeAMIMethodType(AMIPatchToPatchInterpolation, partialFaceAreaWeightAMI);
}

View File

@ -0,0 +1,129 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cyclicACMIGAMGInterfaceField.H"
#include "addToRunTimeSelectionTable.H"
#include "lduMatrix.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(cyclicACMIGAMGInterfaceField, 0);
addToRunTimeSelectionTable
(
GAMGInterfaceField,
cyclicACMIGAMGInterfaceField,
lduInterface
);
addToRunTimeSelectionTable
(
GAMGInterfaceField,
cyclicACMIGAMGInterfaceField,
lduInterfaceField
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cyclicACMIGAMGInterfaceField::cyclicACMIGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const lduInterfaceField& fineInterface
)
:
GAMGInterfaceField(GAMGCp, fineInterface),
cyclicACMIInterface_(refCast<const cyclicACMIGAMGInterface>(GAMGCp)),
doTransform_(false),
rank_(0)
{
const cyclicAMILduInterfaceField& p =
refCast<const cyclicAMILduInterfaceField>(fineInterface);
doTransform_ = p.doTransform();
rank_ = p.rank();
}
Foam::cyclicACMIGAMGInterfaceField::cyclicACMIGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
)
:
GAMGInterfaceField(GAMGCp, doTransform, rank),
cyclicACMIInterface_(refCast<const cyclicACMIGAMGInterface>(GAMGCp)),
doTransform_(doTransform),
rank_(rank)
{}
// * * * * * * * * * * * * * * * * Desstructor * * * * * * * * * * * * * * * //
Foam::cyclicACMIGAMGInterfaceField::~cyclicACMIGAMGInterfaceField()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::cyclicACMIGAMGInterfaceField::updateInterfaceMatrix
(
scalarField& result,
const scalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
) const
{
// Get neighbouring field
scalarField pnf
(
cyclicACMIInterface_.neighbPatch().interfaceInternalField(psiInternal)
);
// Transform according to the transformation tensors
transformCoupleField(pnf, cmpt);
if (cyclicACMIInterface_.owner())
{
pnf = cyclicACMIInterface_.AMI().interpolateToSource(pnf);
}
else
{
pnf = cyclicACMIInterface_.neighbPatch().AMI().interpolateToTarget(pnf);
}
const labelUList& faceCells = cyclicACMIInterface_.faceCells();
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
// ************************************************************************* //

Some files were not shown because too many files have changed in this diff Show More