zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'feature/cvMesh'

Browse Source
This commit is contained in:
laurence
2013-06-04 10:04:21 +01:00
parent 1ad81c9b19 a0a162948f
commit 00d0aeeb39
325 changed files with 112562 additions and 1402 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/OpenFOAM/Make/files
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

16
src/OpenFOAM/algorithms/indexedOctree/indexedOctree.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,

18
src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
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,

45
src/OpenFOAM/fields/Fields/quaternionField/quaternionField.C Normal file
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"
// ************************************************************************* //

64
src/OpenFOAM/fields/Fields/quaternionField/quaternionField.H Normal file
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
// ************************************************************************* //

43
src/OpenFOAM/fields/Fields/quaternionField/quaternionIOField.C Normal file
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
);
}
// ************************************************************************* //

49
src/OpenFOAM/fields/Fields/quaternionField/quaternionIOField.H Normal file
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
// ************************************************************************* //

14
src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
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
|| (

2
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
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.

13
src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
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();

22
src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatchNew.C
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);

3
src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.C
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"

12
src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H
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.

118
src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSortPoints.C Normal file
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;
}
// ************************************************************************* //

72
src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchPointAddressing.C
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;
// }
}

5
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
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;
}
}

10
src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
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"

8
src/OpenFOAM/primitives/Pair/Pair.H
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)
:

147
src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2D.H Normal file
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
// ************************************************************************* //

503
src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2DI.H Normal file
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
// ************************************************************************* //

85
src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.C Normal file
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
// ************************************************************************* //

63
src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.H Normal file
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
// ************************************************************************* //

9
src/OpenFOAM/primitives/Tensor2D/Tensor2D.H
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>&);
};

138
src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H
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> >
{

64
src/OpenFOAM/primitives/quaternion/quaternion.C
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);
}

17
src/OpenFOAM/primitives/quaternion/quaternion.H
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;}

77
src/OpenFOAM/primitives/quaternion/quaternionI.H
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)

138
src/OpenFOAM/primitives/triad/triad.C
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);
}
// ************************************************************************* //

5
src/OpenFOAM/primitives/triad/triad.H
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;}

428
src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
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
@ -1105,83 +1135,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 +1478,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;
}
}

3
src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
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);
};

206
src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
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()

5
src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
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;

68
src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
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

30
src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
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;

24
src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H
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_)
{

36
src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshTemplates.C
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_);

43
src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
View File

@ -403,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)

21
src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
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
@ -153,6 +153,25 @@ public:
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
);
// Member Functions

110
src/meshTools/indexedOctree/treeDataEdge.C
View File

@ -105,17 +105,6 @@ Foam::treeDataEdge::findNearestOp::findNearestOp
{}
Foam::treeDataEdge::findNearestOpSubset::findNearestOpSubset
(
const indexedOctree<treeDataEdge>& tree,
DynamicList<label>& shapeMask
)
:
tree_(tree),
shapeMask_(shapeMask)
{}
Foam::treeDataEdge::findIntersectOp::findIntersectOp
(
const indexedOctree<treeDataEdge>& tree
@ -281,105 +270,6 @@ void Foam::treeDataEdge::findNearestOp::operator()
}
void Foam::treeDataEdge::findNearestOpSubset::operator()
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& minIndex,
point& nearestPoint
) const
{
const treeDataEdge& shape = tree_.shapes();
forAll(indices, i)
{
const label index = indices[i];
const label edgeIndex = shape.edgeLabels()[index];
if (!shapeMask_.empty() && findIndex(shapeMask_, edgeIndex) != -1)
{
continue;
}
const edge& e = shape.edges()[edgeIndex];
pointHit nearHit = e.line(shape.points()).nearestDist(sample);
scalar distSqr = sqr(nearHit.distance());
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
minIndex = index;
nearestPoint = nearHit.rawPoint();
}
}
}
void Foam::treeDataEdge::findNearestOpSubset::operator()
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const
{
const treeDataEdge& shape = tree_.shapes();
// Best so far
scalar nearestDistSqr = magSqr(linePoint - nearestPoint);
forAll(indices, i)
{
const label index = indices[i];
const label edgeIndex = shape.edgeLabels()[index];
if (!shapeMask_.empty() && findIndex(shapeMask_, edgeIndex) != -1)
{
continue;
}
const edge& e = shape.edges()[edgeIndex];
// Note: could do bb test ? Worthwhile?
// Nearest point on line
point ePoint, lnPt;
scalar dist = e.line(shape.points()).nearestDist(ln, ePoint, lnPt);
scalar distSqr = sqr(dist);
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
minIndex = index;
linePoint = lnPt;
nearestPoint = ePoint;
{
point& minPt = tightest.min();
minPt = min(ln.start(), ln.end());
minPt.x() -= dist;
minPt.y() -= dist;
minPt.z() -= dist;
}
{
point& maxPt = tightest.max();
maxPt = max(ln.start(), ln.end());
maxPt.x() += dist;
maxPt.y() += dist;
maxPt.z() += dist;
}
}
}
}
bool Foam::treeDataEdge::findIntersectOp::operator()
(
const label index,

37
src/meshTools/indexedOctree/treeDataEdge.H
View File

@ -118,43 +118,6 @@ public:
) const;
};
class findNearestOpSubset
{
const indexedOctree<treeDataEdge>& tree_;
DynamicList<label>& shapeMask_;
public:
findNearestOpSubset
(
const indexedOctree<treeDataEdge>& tree,
DynamicList<label>& shapeMask
);
void operator()
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& minIndex,
point& nearestPoint
) const;
//- Calculates nearest (to line) point in shape.
// Returns point and distance (squared)
void operator()
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const;
};
class findIntersectOp
{

3
src/meshTools/indexedOctree/treeDataPrimitivePatch.C
View File

@ -116,7 +116,8 @@ bool Foam::treeDataPrimitivePatch<PatchType>::findIntersection
dir,
faceCentres[index],
points,
intersection::HALF_RAY
intersection::HALF_RAY,
shape.planarTol_
);
}

67
src/meshTools/meshTools/meshTools.C
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
@ -215,6 +215,20 @@ void Foam::meshTools::writeOBJ
}
void Foam::meshTools::writeOBJ
(
Ostream& os,
const triad& t,
const point& pt
)
{
forAll(t, dirI)
{
writeOBJ(os, pt, pt + t[dirI]);
}
}
void Foam::meshTools::writeOBJ
(
Ostream& os,
@ -248,57 +262,6 @@ void Foam::meshTools::writeOBJ
}
void Foam::meshTools::writeOBJ
(
Ostream& os,
const faceList& faces,
const pointField& points,
const labelList& faceLabels
)
{
Map<label> foamToObj(4*faceLabels.size());
label vertI = 0;
forAll(faceLabels, i)
{
const face& f = faces[faceLabels[i]];
forAll(f, fp)
{
if (foamToObj.insert(f[fp], vertI))
{
writeOBJ(os, points[f[fp]]);
vertI++;
}
}
os << 'l';
forAll(f, fp)
{
os << ' ' << foamToObj[f[fp]]+1;
}
os << ' ' << foamToObj[f[0]]+1 << endl;
}
}
void Foam::meshTools::writeOBJ
(
Ostream& os,
const faceList& faces,
const pointField& points
)
{
labelList allFaces(faces.size());
forAll(allFaces, i)
{
allFaces[i] = i;
}
writeOBJ(os, faces, points, allFaces);
}
void Foam::meshTools::writeOBJ
(
Ostream& os,

24
src/meshTools/meshTools/meshTools.H
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
@ -37,6 +37,7 @@ SourceFiles
#include "label.H"
#include "vector.H"
#include "triad.H"
#include "labelList.H"
#include "pointField.H"
#include "faceList.H"
@ -107,6 +108,15 @@ namespace meshTools
const point& pt
);
//- Write obj representation of a triad. Requires the location of the
// triad to be supplied
void writeOBJ
(
Ostream& os,
const triad& t,
const point& pt
);
//- Write obj representation of a line connecting two points
// Need to keep track of points that have been added. count starts at 0
void writeOBJ
@ -126,19 +136,21 @@ namespace meshTools
);
//- Write obj representation of faces subset
template<class FaceType>
void writeOBJ
(
Ostream& os,
const faceList&,
const UList<FaceType>&,
const pointField&,
const labelList& faceLabels
);
//- Write obj representation of faces
template<class FaceType>
void writeOBJ
(
Ostream& os,
const faceList&,
const UList<FaceType>&,
const pointField&
);
@ -332,6 +344,12 @@ namespace meshTools
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "meshToolsTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

79
src/meshTools/meshTools/meshToolsTemplates.C Normal file
View File

@ -0,0 +1,79 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
\*---------------------------------------------------------------------------*/
template<class FaceType>
void Foam::meshTools::writeOBJ
(
Ostream& os,
const UList<FaceType>& faces,
const pointField& points,
const labelList& faceLabels
)
{
Map<label> foamToObj(4*faceLabels.size());
label vertI = 0;
forAll(faceLabels, i)
{
const FaceType& f = faces[faceLabels[i]];
forAll(f, fp)
{
if (foamToObj.insert(f[fp], vertI))
{
writeOBJ(os, points[f[fp]]);
vertI++;
}
}
os << 'l';
forAll(f, fp)
{
os << ' ' << foamToObj[f[fp]]+1;
}
os << ' ' << foamToObj[f[0]]+1 << endl;
}
}
template<class FaceType>
void Foam::meshTools::writeOBJ
(
Ostream& os,
const UList<FaceType>& faces,
const pointField& points
)
{
labelList allFaces(faces.size());
forAll(allFaces, i)
{
allFaces[i] = i;
}
writeOBJ(os, faces, points, allFaces);
}
// ************************************************************************* //

1
src/meshTools/searchableSurface/searchableSurface.C
View File

@ -39,7 +39,6 @@ defineRunTimeSelectionTable(searchableSurface, dict);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::searchableSurface> Foam::searchableSurface::New
(
const word& searchableSurfaceType,

16
src/meshTools/searchableSurface/triSurfaceMesh.C
View File

@ -742,9 +742,19 @@ void Foam::triSurfaceMesh::getVolumeType
{
const point& pt = points[pointI];
// - use cached volume type per each tree node
// - cheat conversion since same values
volType[pointI] = tree().getVolumeType(pt);
if (!tree().bb().contains(pt))
{
// Have to calculate directly as outside the octree
volType[pointI] = tree().shapes().getVolumeType(tree(), pt);
}
else
{
// - use cached volume type per each tree node
volType[pointI] = tree().getVolumeType(pt);
}
// Info<< "octree : " << pt << " = "
// << volumeType::names[volType[pointI]] << endl;
}
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;

91
src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
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
@ -118,7 +118,11 @@ Foam::List<Foam::surfaceFeatures::edgeStatus> Foam::surfaceFeatures::toStatus()
// Set from value for every edge
void Foam::surfaceFeatures::setFromStatus(const List<edgeStatus>& edgeStat)
void Foam::surfaceFeatures::setFromStatus
(
const List<edgeStatus>& edgeStat,
const scalar includedAngle
)
{
// Count
@ -170,20 +174,24 @@ void Foam::surfaceFeatures::setFromStatus(const List<edgeStatus>& edgeStat)
}
}
// Calculate consistent feature points
calcFeatPoints(edgeStat);
const scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
calcFeatPoints(edgeStat, minCos);
}
//construct feature points where more than 2 feature edges meet
void Foam::surfaceFeatures::calcFeatPoints
(
const List<edgeStatus>& edgeStat
const List<edgeStatus>& edgeStat,
const scalar minCos
)
{
DynamicList<label> featurePoints(surf_.nPoints()/1000);
const labelListList& pointEdges = surf_.pointEdges();
const edgeList& edges = surf_.edges();
const pointField& localPoints = surf_.localPoints();
forAll(pointEdges, pointI)
{
@ -203,6 +211,27 @@ void Foam::surfaceFeatures::calcFeatPoints
{
featurePoints.append(pointI);
}
else if (nFeatEdges == 2)
{
// Check the angle between the two edges
DynamicList<vector> edgeVecs(2);
forAll(pEdges, i)
{
const label edgeI = pEdges[i];
if (edgeStat[edgeI] != NONE)
{
edgeVecs.append(edges[edgeI].vec(localPoints));
edgeVecs.last() /= mag(edgeVecs.last());
}
}
if (mag(edgeVecs[0] & edgeVecs[1]) < minCos)
{
featurePoints.append(pointI);
}
}
}
featurePoints_.transfer(featurePoints);
@ -213,7 +242,8 @@ void Foam::surfaceFeatures::classifyFeatureAngles
(
const labelListList& edgeFaces,
List<edgeStatus>& edgeStat,
const scalar minCos
const scalar minCos,
const bool geometricTestOnly
) const
{
const vectorField& faceNormals = surf_.faceNormals();
@ -229,14 +259,18 @@ void Foam::surfaceFeatures::classifyFeatureAngles
if (eFaces.size() != 2)
{
// Non-manifold. What to do here? Is region edge? external edge?
edgeStat[edgeI] = REGION;
edgeStat[edgeI] = NONE;
}
else
{
label face0 = eFaces[0];
label face1 = eFaces[1];
if (surf_[face0].region() != surf_[face1].region())
if
(
!geometricTestOnly
&& surf_[face0].region() != surf_[face1].region()
)
{
edgeStat[edgeI] = REGION;
}
@ -445,7 +479,8 @@ Foam::surfaceFeatures::surfaceFeatures
const triSurface& surf,
const scalar includedAngle,
const scalar minLen,
const label minElems
const label minElems,
const bool geometricTestOnly
)
:
surf_(surf),
@ -454,11 +489,11 @@ Foam::surfaceFeatures::surfaceFeatures
externalStart_(0),
internalStart_(0)
{
findFeatures(includedAngle);
findFeatures(includedAngle, geometricTestOnly);
if (minLen > 0 || minElems > 0)
{
trimFeatures(minLen, minElems);
trimFeatures(minLen, minElems, includedAngle);
}
}
@ -507,7 +542,8 @@ Foam::surfaceFeatures::surfaceFeatures
const triSurface& surf,
const pointField& points,
const edgeList& edges,
const scalar mergeTol
const scalar mergeTol,
const bool geometricTestOnly
)
:
surf_(surf),
@ -553,7 +589,13 @@ Foam::surfaceFeatures::surfaceFeatures
// Find whether an edge is external or internal
List<edgeStatus> edgeStat(dynFeatEdges.size(), NONE);
classifyFeatureAngles(dynFeatureEdgeFaces, edgeStat, GREAT);
classifyFeatureAngles
(
dynFeatureEdgeFaces,
edgeStat,
GREAT,
geometricTestOnly
);
// Transfer the edge status to a list encompassing all edges in the surface
// so that calcFeatPoints can be used.
@ -572,7 +614,7 @@ Foam::surfaceFeatures::surfaceFeatures
edgeStat.clear();
dynFeatEdges.clear();
setFromStatus(allEdgeStat);
setFromStatus(allEdgeStat, GREAT);
}
@ -632,16 +674,26 @@ Foam::labelList Foam::surfaceFeatures::selectFeatureEdges
}
void Foam::surfaceFeatures::findFeatures(const scalar includedAngle)
void Foam::surfaceFeatures::findFeatures
(
const scalar includedAngle,
const bool geometricTestOnly
)
{
scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
// Per edge whether is feature edge.
List<edgeStatus> edgeStat(surf_.nEdges(), NONE);
classifyFeatureAngles(surf_.edgeFaces(), edgeStat, minCos);
classifyFeatureAngles
(
surf_.edgeFaces(),
edgeStat,
minCos,
geometricTestOnly
);
setFromStatus(edgeStat);
setFromStatus(edgeStat, includedAngle);
}
@ -650,7 +702,8 @@ void Foam::surfaceFeatures::findFeatures(const scalar includedAngle)
Foam::labelList Foam::surfaceFeatures::trimFeatures
(
const scalar minLen,
const label minElems
const label minElems,
const scalar includedAngle
)
{
// Convert feature edge list to status per edge.
@ -774,7 +827,7 @@ Foam::labelList Foam::surfaceFeatures::trimFeatures
}
// Convert back to edge labels
setFromStatus(edgeStat);
setFromStatus(edgeStat, includedAngle);
return featLines;
}

40
src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
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
@ -133,14 +133,19 @@ private:
//- Construct feature points where more than 2 feature edges meet
void calcFeatPoints(const List<edgeStatus>&);
void calcFeatPoints
(
const List<edgeStatus>& edgeStat,
const scalar minCos
);
//- Classify the angles of the feature edges
void classifyFeatureAngles
(
const labelListList& edgeFaces,
List<edgeStatus>& edgeStat,
const scalar minCos
const scalar minCos,
const bool geometricTestOnly
) const;
//- Choose next unset feature edge.
@ -184,13 +189,16 @@ public:
);
//- Construct from surface, angle and min cumulative length and/or
// number of elements
// number of elements. If geometric test only is true, then region
// information is ignored and features are only assigned based on the
// geometric criteria
surfaceFeatures
(
const triSurface&,
const scalar includedAngle,
const scalar minLen = 0,
const label minElems = 0
const label minElems = 0,
const bool geometricTestOnly = false
);
//- Construct from dictionary
@ -205,7 +213,8 @@ public:
const triSurface&,
const pointField& points,
const edgeList& edges,
const scalar mergeTol = 1e-6
const scalar mergeTol = 1e-6,
const bool geometricTestOnly = false
);
//- Construct as copy
@ -275,17 +284,30 @@ public:
// Edit
//- Find feature edges using provided included angle
void findFeatures(const scalar includedAngle);
void findFeatures
(
const scalar includedAngle,
const bool geometricTestOnly
);
//- Delete small sets of edges. Edges are stringed up and any
// string of length < minLen (or nElems < minElems) is deleted.
labelList trimFeatures(const scalar minLen, const label minElems);
labelList trimFeatures
(
const scalar minLen,
const label minElems,
const scalar includedAngle
);
//- From member feature edges to status per edge.
List<edgeStatus> toStatus() const;
//- Set from status per edge
void setFromStatus(const List<edgeStatus>&);
void setFromStatus
(
const List<edgeStatus>& edgeStat,
const scalar includedAngle
);
// Find

4
src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
View File

@ -51,7 +51,7 @@ bool Foam::triSurfaceSearch::checkUniqueHit
nearLabel
);
if (nearType == 1)
if (nearType == triPointRef::POINT)
{
// near point
@ -78,7 +78,7 @@ bool Foam::triSurfaceSearch::checkUniqueHit
}
}
}
else if (nearType == 2)
else if (nearType == triPointRef::EDGE)
{
// near edge
// check if the other face of the edge is already hit