Merge branch 'master' into molecularDynamics

This commit is contained in:
graham
2008-08-19 13:57:22 +01:00
132 changed files with 4356 additions and 1588 deletions

View File

@ -1,9 +1,11 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
( cd OpenFOAM && wmakeLnInclude . )
( cd Pstream && ./Allwmake )
wmake libo OSspecific/$WM_OS
wmake libo OSspecific/$WM_OS
wmake libso OpenFOAM
wmake libso lagrangian/basic

View File

@ -0,0 +1,129 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::Keyed
Description
A container with an integer key attached to any item.
The key can useful for sorting.
SourceFiles
KeyedI.H
\*---------------------------------------------------------------------------*/
#ifndef Keyed_H
#define Keyed_H
#include "List.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of friend functions and operators
template<class T> class Keyed;
template<class T> Istream& operator>>(Istream&, Keyed<T>&);
template<class T> Ostream& operator<<(Ostream&, const Keyed<T>&);
/*---------------------------------------------------------------------------*\
Class Keyed Declaration
\*---------------------------------------------------------------------------*/
template<class T>
class Keyed
:
public T
{
// Private data
label key_;
public:
// Static Members
//- Add labels to a list of values
inline static List<Keyed<T> > createList
(
const List<T>&,
const label key=0
);
//- Add labels to a list of values
inline static List<Keyed<T> > createList
(
const List<T>&,
const List<label>& keys
);
// Constructors
//- Construct null
inline Keyed();
//- Construct as a copy of item, with a key
inline Keyed(const T& item, const label key=0);
//- Construct from Istream
inline Keyed(Istream&);
// Member Functions
// Access
//- Return const access to the integer key
inline label key() const;
//- Return non-const access to the integer key
inline label& key();
// IOstream Operators
friend Istream& operator>> <T>(Istream&, Keyed<T>&);
friend Ostream& operator<< <T>(Ostream&, const Keyed<T>&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "KeyedI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,146 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
//- Construct null
template<class T>
inline Foam::Keyed<T>::Keyed()
:
key_(-1)
{}
//- Construct from components
template<class T>
inline Foam::Keyed<T>::Keyed(const T& item, const label key)
:
T(item),
key_(key)
{}
template<class T>
inline Foam::Keyed<T>::Keyed(Istream& is)
{
is >> *this;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class T>
inline Foam::label Foam::Keyed<T>::key() const
{
return key_;
}
template<class T>
inline Foam::label& Foam::Keyed<T>::key()
{
return key_;
}
template<class T>
inline Foam::List<Foam::Keyed<T> >
Foam::Keyed<T>::createList(const List<T>& lst, const label key)
{
List<Keyed<T> > newList(lst.size());
forAll (lst, elemI)
{
newList[elemI] = Keyed(lst[elemI], key);
}
return newList;
}
template<class T>
inline Foam::List<Foam::Keyed<T> >
Foam::Keyed<T>::createList(const List<T>& lst, const List<label>& keys)
{
if (lst.size() != keys.size())
{
FatalErrorIn
(
"Foam::Keyed<T>::createList(const List<T>&, const List<label>&)"
)
<< "size mismatch adding keys to a list:" << nl
<< "List has size " << lst.size()
<< " and keys has size " << keys.size() << nl
<< abort(FatalError);
}
List<Keyed<T> > newList(lst.size());
forAll (lst, elemI)
{
newList[elemI] = Keyed(lst[elemI], keys[elemI]);
}
return newList;
}
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template<class T>
inline Foam::Istream& Foam::operator>>(Istream& is, Keyed<T>& item)
{
// Read beginning of Keyed item/key pair
is.readBegin("Keyed<T>");
is >> static_cast<T&>(item);
is >> item.key_;
// Read end of Keyed item/key pair
is.readEnd("Keyed<T>");
// Check state of Ostream
is.check("Istream& operator>>(Istream&, Keyed&)");
return is;
}
template<class T>
inline Foam::Ostream& Foam::operator<<(Ostream& os, const Keyed<T>& item)
{
os << token::BEGIN_LIST
<< static_cast<const T&>(item)
<< token::SPACE << item.key_
<< token::END_LIST;
return os;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -53,14 +53,9 @@ inline Foam::label Foam::globalIndex::toGlobal(const label i) const
//- Is on local processor
inline bool Foam::globalIndex::isLocal(const label i) const
{
label localI =
(
Pstream::myProcNo() == 0
? i
: i - offsets_[Pstream::myProcNo()-1]
);
return localI >= 0 && localI < offsets_[Pstream::myProcNo()];
return
(i < offsets_[Pstream::myProcNo()])
&& (i >= (Pstream::myProcNo() == 0 ? 0 : offsets_[Pstream::myProcNo()-1]));
}
@ -69,7 +64,7 @@ inline Foam::label Foam::globalIndex::toLocal(const label procI, const label i)
{
label localI = (procI == 0 ? i : i - offsets_[procI-1]);
if (localI < 0 || localI >= offsets_[procI])
if (localI < 0 || i >= offsets_[procI])
{
FatalErrorIn("globalIndex::toLocal(const label, const label)")
<< "Global " << i << " does not belong on processor "

View File

@ -0,0 +1,173 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "PrimitivePatchExtra.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
PrimitivePatchExtra
(
const ListType<Face>& faces,
const pointField& points
)
:
PrimitivePatch<Face, ListType, PointField, PointType>(faces, points),
sortedEdgeFacesPtr_(NULL),
edgeOwnerPtr_(NULL)
{}
// Construct as copy
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
PrimitivePatchExtra
(
const PrimitivePatchExtra<Face, ListType, PointField, PointType>& pp
)
:
PrimitivePatch<Face, ListType, PointField, PointType>(pp),
sortedEdgeFacesPtr_(NULL),
edgeOwnerPtr_(NULL)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
~PrimitivePatchExtra()
{
clearOut();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
clearOut()
{
PrimitivePatch<Face, ListType, PointField, PointType>::clearOut();
clearTopology();
}
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
clearTopology()
{
PrimitivePatch<Face, ListType, PointField, PointType>::clearTopology();
deleteDemandDrivenData(sortedEdgeFacesPtr_);
deleteDemandDrivenData(edgeOwnerPtr_);
}
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
const Foam::labelListList&
Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
sortedEdgeFaces() const
{
if (!sortedEdgeFacesPtr_)
{
calcSortedEdgeFaces();
}
return *sortedEdgeFacesPtr_;
}
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
const Foam::labelList&
Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
edgeOwner() const
{
if (!edgeOwnerPtr_)
{
calcEdgeOwner();
}
return *edgeOwnerPtr_;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "PrimitivePatchExtraAddressing.C"
#include "PrimitivePatchExtraCleanup.C"
#include "PrimitivePatchExtraSearch.C"
// ************************************************************************* //

View File

@ -0,0 +1,208 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::PrimitivePatchExtra
Description
PrimitivePatch with some extra functionality.
SourceFiles
PrimitivePatchExtra.C
\*---------------------------------------------------------------------------*/
#ifndef PrimitivePatchExtra_H
#define PrimitivePatchExtra_H
#include "PrimitivePatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class PrimitivePatchExtra Declaration
\*---------------------------------------------------------------------------*/
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType=point
>
class PrimitivePatchExtra
:
public PrimitivePatch<Face, ListType, PointField, PointType>
{
public:
// Public typedefs
typedef Face FaceType;
typedef ListType<Face> FaceListType;
typedef PointField PointFieldType;
private:
// Private typedefs
typedef PrimitivePatch<Face, ListType, PointField, PointType> TemplateType;
// Private data
// Demand driven private data
//- Edge-face addressing (sorted)
mutable labelListList* sortedEdgeFacesPtr_;
//- Label of face that 'owns' edge
// i.e. e.vec() is righthanded walk along face
mutable labelList* edgeOwnerPtr_;
// Private Member Functions
//- Calculate sorted edgeFaces
void calcSortedEdgeFaces() const;
//- Calculate owner
void calcEdgeOwner() const;
public:
// Constructors
//- Construct from components
PrimitivePatchExtra
(
const ListType<Face>& faces,
const pointField& points
);
//- Construct as copy
PrimitivePatchExtra
(
const PrimitivePatchExtra<Face, ListType, PointField, PointType>&
);
// Destructor
virtual ~PrimitivePatchExtra();
void clearOut();
void clearTopology();
// Member Functions
// Access functions for demand driven data
// Topological data; no mesh required.
//- Return edge-face addressing sorted
// (for edges with more than 2 faces) according to the
// angle around the edge.
// Orientation is anticlockwise looking from
// edge.vec(localPoints())
const labelListList& sortedEdgeFaces() const;
//- If 2 face neighbours: label of face where ordering of edge
// is consistent with righthand walk.
// If 1 neighbour: label of only face.
// If >2 neighbours: undetermined.
const labelList& edgeOwner() const;
// Addressing into mesh
// Other patch operations
// Check
//- Check multiply-connected edges.
void checkEdges(const bool verbose) const;
//- Check orientation (normals) and normals of neighbouring faces
boolList checkOrientation(const bool verbose) const;
// Edit
//- Fill faceZone with currentZone for every face reachable
// from faceI without crossing edge marked in borderEdge.
// Note: faceZone has to be sized nFaces before calling.
void markZone
(
const boolList& borderEdge,
const label faceI,
const label currentZone,
labelList& faceZone
) const;
//- (size and) fills faceZone with zone of face.
// Zone is area reachable by edge crossing without crossing borderEdge
// (bool for every edge in surface). Returns number of zones.
label markZones
(
const boolList& borderEdge,
labelList& faceZone
) const;
//- Determine the mapping for a sub-mesh.
// Only include faces for which boolList entry is true
// @param[out] pointMap mapping new to old localPoints
// @param[out] faceMap mapping new to old faces
void subsetMap
(
const boolList& include,
labelList& pointMap,
labelList& faceMap
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "PrimitivePatchExtra.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,210 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Contains fix for PrimitivePatch addressing (which doesn't work if surface
is non-manifold). Should be moved into PrimitivePatch.
\*---------------------------------------------------------------------------*/
#include "PrimitivePatchExtra.H"
#include "HashTable.H"
#include "SortableList.H"
#include "transform.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
calcSortedEdgeFaces() const
{
if (sortedEdgeFacesPtr_)
{
FatalErrorIn("PrimitivePatchExtra<Face, ListType, PointField>::calcSortedEdgeFaces()")
<< "sortedEdgeFacesPtr_ already set"
<< abort(FatalError);
}
const labelListList& eFaces = TemplateType::edgeFaces();
const edgeList& edgeLst = TemplateType::edges();
const pointField& locPointLst = TemplateType::localPoints();
const List<FaceType>& locFaceLst = TemplateType::localFaces();
// create the lists for the various results. (resized on completion)
sortedEdgeFacesPtr_ = new labelListList(eFaces.size());
labelListList& sortedEdgeFaces = *sortedEdgeFacesPtr_;
forAll(eFaces, edgeI)
{
const labelList& myFaceNbs = eFaces[edgeI];
if (myFaceNbs.size() > 2)
{
// Get point on edge and normalized direction of edge (= e2 base
// of our coordinate system)
const edge& e = edgeLst[edgeI];
const point& edgePt = locPointLst[e.start()];
vector e2 = e.vec(locPointLst);
e2 /= mag(e2) + VSMALL;
// Get opposite vertex for 0th face
const Face& f = locFaceLst[myFaceNbs[0]];
label fp0 = findIndex(f, e[0]);
label fp1 = f.fcIndex(fp0);
label vertI = (f[fp1] != e[1] ? f[fp1] : f.fcIndex(fp1));
// Get vector normal both to e2 and to edge from opposite vertex
// to edge (will be x-axis of our coordinate system)
vector e0 = e2 ^ (locPointLst[vertI] - edgePt);
e0 /= mag(e0) + VSMALL;
// Get y-axis of coordinate system
vector e1 = e2 ^ e0;
SortableList<scalar> faceAngles(myFaceNbs.size());
// e0 is reference so angle is 0
faceAngles[0] = 0;
for (label nbI = 1; nbI < myFaceNbs.size(); nbI++)
{
// Get opposite vertex
const FaceType& f = locFaceLst[myFaceNbs[nbI]];
label fp0 = findIndex(f, e[0]);
label fp1 = f.fcIndex(fp0);
label vertI = (f[fp1] != e[1] ? f[fp1] : f.fcIndex(fp1));
vector vec = e2 ^ (locPointLst[vertI] - edgePt);
vec /= mag(vec) + VSMALL;
faceAngles[nbI] = pseudoAngle
(
e0,
e1,
vec
);
}
faceAngles.sort();
sortedEdgeFaces[edgeI] = IndirectList<label>
(
myFaceNbs,
faceAngles.indices()
);
}
else
{
// No need to sort. Just copy.
sortedEdgeFaces[edgeI] = myFaceNbs;
}
}
}
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
calcEdgeOwner() const
{
if (edgeOwnerPtr_)
{
FatalErrorIn("PrimitivePatchExtra<Face, ListType, PointField>::calcEdgeOwner()")
<< "edgeOwnerPtr_ already set"
<< abort(FatalError);
}
// create the owner list
edgeOwnerPtr_ = new labelList
(
TemplateType::nEdges()
);
labelList& edgeOwner = *edgeOwnerPtr_;
const edgeList& edgeLst = TemplateType::edges();
const labelListList& eFaces = TemplateType::edgeFaces();
const List<FaceType>& locFaceLst = TemplateType::localFaces();
forAll(edgeLst, edgeI)
{
const edge& e = edgeLst[edgeI];
const labelList& myFaces = eFaces[edgeI];
if (myFaces.size() == 1)
{
edgeOwner[edgeI] = myFaces[0];
}
else
{
// Find the first face whose vertices are aligned with the edge.
// (in case of multiply connected edge the best we can do)
edgeOwner[edgeI] = -1;
forAll(myFaces, i)
{
const FaceType& f = locFaceLst[myFaces[i]];
if (f.findEdge(e) > 0)
{
edgeOwner[edgeI] = myFaces[i];
break;
}
}
if (edgeOwner[edgeI] == -1)
{
FatalErrorIn("PrimitivePatchExtra<Face, ListType, PointField>::calcEdgeOwner()")
<< "Edge " << edgeI << " vertices:" << e
<< " is used by faces " << myFaces
<< " vertices:"
<< IndirectList<FaceType>(locFaceLst, myFaces)()
<< " none of which use the edge vertices in the same order"
<< nl << "I give up" << abort(FatalError);
}
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -0,0 +1,214 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "PrimitivePatchExtra.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Check/fix edges with more than two faces
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
checkEdges
(
const bool verbose
) const
{
const labelListList& eFaces = TemplateType::edgeFaces();
const edgeList& edgeLst = TemplateType::edges();
forAll(eFaces, edgeI)
{
const labelList& myFaces = eFaces[edgeI];
// boundary edges have one face
// interior edges have two faces
if (myFaces.size() == 0)
{
FatalErrorIn("PrimitivePatchExtra::checkEdges(bool verbose)")
<< "Edge " << edgeI << " with vertices " << edgeLst[edgeI]
<< " has no edgeFaces"
<< exit(FatalError);
}
else if (myFaces.size() > 2)
{
WarningIn
(
"PrimitivePatchExtra::checkEdges(bool verbose)"
) << "Edge " << edgeI << " with vertices " << edgeLst[edgeI]
<< " has more than 2 faces connected to it : " << myFaces
<< endl;
}
}
}
// Check normals and orientation
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
Foam::boolList
Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
checkOrientation
(
const bool verbose
) const
{
const ListType<FaceType>& faceLst = *this;
const edgeList& edgeLst = TemplateType::edges();
const labelListList& faceEs = TemplateType::faceEdges();
const label numEdges = TemplateType::nEdges();
const pointField& pointLst = TemplateType::points();
const vectorField& normLst = TemplateType::faceNormals();
// Check edge normals, face normals, point normals.
forAll(faceEs, faceI)
{
const labelList& edgeLabels = faceEs[faceI];
if (edgeLabels.size() < 3)
{
FatalErrorIn("PrimitivePatchExtra::checkOrientation(bool)")
<< "face " << faceLst[faceI]
<< " has fewer than 3 edges. Edges:" << edgeLabels
<< exit(FatalError);
}
bool valid = true;
forAll(edgeLabels, i)
{
if (edgeLabels[i] < 0 || edgeLabels[i] >= numEdges)
{
WarningIn
(
"PrimitivePatchExtra::checkOrientation(bool)"
) << "edge number " << edgeLabels[i] << " on face " << faceI
<< " out-of-range\n"
<< "This usually means that the input surface has "
<< "edges with more than 2 faces connected.\n"
<< endl;
valid = false;
}
}
if (!valid)
{
continue;
}
//
//- Compute normal from 3 points, use the first as the origin
// minor warpage should not be a problem
const FaceType& f = faceLst[faceI];
const point p0(pointLst[f[0]]);
const point p1(pointLst[f[1]]);
const point p2(pointLst[f[f.size()-1]]);
const vector pointNormal((p1 - p0) ^ (p2 - p0));
if ((pointNormal & normLst[faceI]) < 0)
{
FatalErrorIn("PrimitivePatchExtra::checkOrientation(bool)")
<< "Normal calculated from points not consistent with"
" faceNormal" << endl
<< "face: " << f << endl
<< "points: " << p0 << ' ' << p1 << ' ' << p2 << endl
<< "pointNormal:" << pointNormal << endl
<< "faceNormal:" << normLst[faceI]
<< exit(FatalError);
}
}
const labelListList& eFaces = TemplateType::edgeFaces();
const pointField& locPointsLst = TemplateType::localPoints();
// Storage for holding status of edge. True if normal flips across this
// edge
boolList borderEdge(numEdges, false);
forAll(edgeLst, edgeI)
{
const edge& e = edgeLst[edgeI];
const labelList& neighbours = eFaces[edgeI];
if (neighbours.size() == 2)
{
const FaceType& faceA = faceLst[neighbours[0]];
const FaceType& faceB = faceLst[neighbours[1]];
// The edge cannot be going in the same direction if both faces
// are oriented counterclockwise.
// Thus the next face point *must* different between the faces.
if
(
faceA[faceA.fcIndex(findIndex(faceA, e.start()))]
== faceB[faceB.fcIndex(findIndex(faceB, e.start()))]
)
{
borderEdge[edgeI] = true;
if (verbose)
{
WarningIn("PrimitivePatchExtra::checkOrientation(bool)")
<< "face orientation incorrect." << nl
<< "edge[" << edgeI << "] " << e
<< " between faces " << neighbours << ":" << nl
<< "face[" << neighbours[0] << "] " << faceA << nl
<< "face[" << neighbours[1] << "] " << faceB << endl;
}
}
}
else if (neighbours.size() != 1)
{
if (verbose)
{
WarningIn("PrimitivePatchExtra::checkOrientation(bool)")
<< "Wrong number of edge neighbours." << endl
<< "edge[" << edgeI << "] " << e
<< "with points:" << locPointsLst[e.start()]
<< ' ' << locPointsLst[e.end()]
<< " has neighbours:" << neighbours << endl;
}
borderEdge[edgeI] = true;
}
}
return borderEdge;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -0,0 +1,233 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "PrimitivePatchExtra.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Finds area, starting at faceI, delimited by borderEdge. Marks all visited
// faces (from face-edge-face walk) with currentZone.
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::markZone
(
const boolList& borderEdge,
const label faceI,
const label currentZone,
labelList& faceZone
) const
{
// List of faces whose faceZone has been set.
labelList changedFaces(1, faceI);
const labelListList& faceEs = TemplateType::faceEdges();
const labelListList& eFaces = TemplateType::edgeFaces();
while (true)
{
// Pick up neighbours of changedFaces
DynamicList<label> newChangedFaces(2*changedFaces.size());
forAll(changedFaces, i)
{
label faceI = changedFaces[i];
const labelList& fEdges = faceEs[faceI];
forAllfEdges, i)
{
label edgeI = fEdges[i];
if (!borderEdge[edgeI])
{
const labelList& eFaceLst = eFaces[edgeI];
forAll(eFaceLst, j)
{
label nbrFaceI = eFaceLst[j];
if (faceZone[nbrFaceI] == -1)
{
faceZone[nbrFaceI] = currentZone;
newChangedFaces.append(nbrFaceI);
}
else if (faceZone[nbrFaceI] != currentZone)
{
FatalErrorIn
(
"PrimitivePatchExtra<Face, ListType, PointField>::markZone"
"(const boolList&, const label, const label, labelList&) const"
)
<< "Zones " << faceZone[nbrFaceI]
<< " at face " << nbrFaceI
<< " connects to zone " << currentZone
<< " at face " << faceI
<< abort(FatalError);
}
}
}
}
}
if (newChangedFaces.size() == 0)
{
break;
}
changedFaces.transfer(newChangedFaces.shrink());
newChangedFaces.clear();
}
}
// Finds areas delimited by borderEdge (or 'real' edges).
// Fills faceZone accordingly
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
Foam::label Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
markZones
(
const boolList& borderEdge,
labelList& faceZone
) const
{
const label numEdges = TemplateType::nEdges();
const label numFaces = TemplateType::size();
faceZone.setSize(numFaces);
faceZone = -1;
if (borderEdge.size() != numEdges)
{
FatalErrorIn
(
"PrimitivePatchExtra<Face, ListType, PointField>::markZones"
"(const boolList&, labelList&)"
)
<< "borderEdge boolList not same size as number of edges" << endl
<< "borderEdge:" << borderEdge.size() << endl
<< "nEdges :" << numEdges
<< exit(FatalError);
}
label zoneI = 0;
label startFaceI = 0;
for (;;zoneI++)
{
// Find first non-coloured face
for (; startFaceI < numFaces; startFaceI++)
{
if (faceZone[startFaceI] == -1)
{
break;
}
}
if (startFaceI >= numFaces)
{
break;
}
faceZone[startFaceI] = zoneI;
markZone(borderEdge, startFaceI, zoneI, faceZone);
}
return zoneI;
}
// Finds areas delimited by borderEdge (or 'real' edges).
// Fills faceZone accordingly
template
<
class Face,
template<class> class ListType,
class PointField,
class PointType
>
void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::
subsetMap
(
const boolList& include,
labelList& pointMap,
labelList& faceMap
) const
{
const List<FaceType>& locFaces = TemplateType::localFaces();
const label numPoints = TemplateType::nPoints();
label faceI = 0;
label pointI = 0;
faceMap.setSize(locFaces.size());
pointMap.setSize(numPoints);
boolList pointHad(numPoints, false);
forAll(include, oldFaceI)
{
if (include[oldFaceI])
{
// Store new faces compact
faceMap[faceI++] = oldFaceI;
// Renumber labels for face
const FaceType& f = locFaces[oldFaceI];
forAll(f, fp)
{
const label ptLabel = f[fp];
if (!pointHad[ptLabel])
{
pointHad[ptLabel] = true;
pointMap[pointI++] = ptLabel;
}
}
}
}
// Trim
faceMap.setSize(faceI);
pointMap.setSize(pointI);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -33,19 +33,15 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#if INT_MAX != 2147483647
# error "INT_MAX != 2147483647"
# error "The random number generator may not work!"
#endif
#ifdef USE_RANDOM
# include <climits>
# if INT_MAX != 2147483647
# error "INT_MAX != 2147483647"
# error "The random number generator random() may not work!"
# endif
#else
# include <cstdlib>
#endif
@ -77,7 +73,7 @@ int Random::bit()
# ifdef USE_RANDOM
if (random() > INT_MAX/2)
# else
if (mrand48() > 0)
if (lrand48() > INT_MAX/2)
# endif
{
return 1;

View File

@ -533,7 +533,7 @@ Foam::labelList Foam::meshRefinement::getRefineCandidateFaces
const labelList& refineCell
) const
{
labelList testFaces(mesh_.nCells());
labelList testFaces(mesh_.nFaces());
label nTest = 0;

View File

@ -2,7 +2,7 @@
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.openfoam.org |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
@ -12,7 +12,6 @@ FoamFile
class dictionary;
object dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//dynamicFvMeshLib "libtopoChangerFvMesh.so";
@ -38,7 +37,6 @@ mixerFvMeshCoeffs
}
}
// Refinement
dynamicRefineFvMeshCoeffs
{

View File

@ -46,6 +46,13 @@ void Foam::Cloud<ParticleType>::initCloud(const bool checkClass)
readFields();
}
}
else
{
WarningIn("Cloud<ParticleType>::initCloud(const bool checkClass)")
<< "Cannot read particle positions file " << nl
<< " " << ioP.path() << nl
<< " assuming the initial cloud contains 0 particles." << endl;
}
}

View File

@ -158,7 +158,7 @@ public:
//- Return velocity
inline const vector& U() const;
//- The nearest distance to a wall that
//- The nearest distance to a wall that
// the particle can be in the n direction
inline scalar wallImpactDistance(const vector& n) const;
@ -226,6 +226,14 @@ public:
const vector& separation
);
// I-O
static void readFields(Cloud<solidParticle>& c);
static void writeFields(const Cloud<solidParticle>& c);
// Ostream Operator
friend Ostream& operator<<(Ostream&, const solidParticle&);
@ -238,12 +246,6 @@ inline bool contiguous<solidParticle>()
return true;
}
template<>
void Cloud<solidParticle>::readFields();
template<>
void Cloud<solidParticle>::writeFields() const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -40,9 +40,13 @@ namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::solidParticleCloud::solidParticleCloud(const fvMesh& mesh)
Foam::solidParticleCloud::solidParticleCloud
(
const fvMesh& mesh,
const word& cloudName
)
:
Cloud<solidParticle>(mesh),
Cloud<solidParticle>(mesh, cloudName, false),
mesh_(mesh),
particleProperties_
(
@ -58,7 +62,9 @@ Foam::solidParticleCloud::solidParticleCloud(const fvMesh& mesh)
rhop_(dimensionedScalar(particleProperties_.lookup("rhop")).value()),
e_(dimensionedScalar(particleProperties_.lookup("e")).value()),
mu_(dimensionedScalar(particleProperties_.lookup("mu")).value())
{}
{
solidParticle::readFields(*this);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -82,4 +88,10 @@ void Foam::solidParticleCloud::move(const dimensionedVector& g)
}
void Foam::solidParticleCloud::writeFields() const
{
solidParticle::writeFields(*this);
}
// ************************************************************************* //

View File

@ -83,7 +83,11 @@ public:
// Constructors
//- Construct given mesh
solidParticleCloud(const fvMesh&);
solidParticleCloud
(
const fvMesh&,
const word& cloudName = "defaultCloud"
);
// Member Functions
@ -97,13 +101,17 @@ public:
inline scalar mu() const;
// Check
// Edit
//- Move the particles under the influence of the given
// gravitational acceleration
void move(const dimensionedVector& g);
//- Move the particles under the influence of the given
// gravitational acceleration
void move(const dimensionedVector& g);
// I-O
//- Write fields
void writeFields() const;
};

View File

@ -36,7 +36,7 @@ Foam::solidParticle::solidParticle
bool readFields
)
:
Particle<solidParticle>(cloud, is)
Particle<solidParticle>(cloud, is, readFields)
{
if (readFields)
{
@ -60,20 +60,20 @@ Foam::solidParticle::solidParticle
}
namespace Foam
void Foam::solidParticle::readFields(Cloud<solidParticle>& c)
{
if (!c.size())
{
return;
}
IOField<scalar> d(c.fieldIOobject("d"));
c.checkFieldIOobject(c, d);
template<>
void Cloud<solidParticle>::readFields()
{
IOField<scalar> d(fieldIOobject("d"));
checkFieldIOobject(*this, d);
IOField<vector> U(fieldIOobject("U"));
checkFieldIOobject(*this, U);
IOField<vector> U(c.fieldIOobject("U"));
c.checkFieldIOobject(c, U);
label i = 0;
forAllIter(Cloud<solidParticle>::iterator, *this, iter)
forAllIter(Cloud<solidParticle>::iterator, c, iter)
{
solidParticle& p = iter();
@ -84,16 +84,17 @@ void Cloud<solidParticle>::readFields()
}
template<>
void Cloud<solidParticle>::writeFields() const
void Foam::solidParticle::writeFields(const Cloud<solidParticle>& c)
{
label np = size();
Particle<solidParticle>::writeFields(c);
IOField<scalar> d(fieldIOobject("d"), np);
IOField<vector> U(fieldIOobject("U"), np);
label np = c.size();
IOField<scalar> d(c.fieldIOobject("d"), np);
IOField<vector> U(c.fieldIOobject("U"), np);
label i = 0;
forAllConstIter(Cloud<solidParticle>, *this, iter)
forAllConstIter(Cloud<solidParticle>, c, iter)
{
const solidParticle& p = iter();
@ -106,8 +107,6 @@ void Cloud<solidParticle>::writeFields() const
U.write();
}
};
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //

View File

@ -1,103 +0,0 @@
ALL FILES:
* construct null
* allow assignment of name/origin
* MAJOR CHANGE
- a point is a vector, but not vice versa
- previous methods toGlobal/toLocal include the origin translation
which is incorrect for vectors
- new methods with explicit names:
globalPosition, globalVector, relativePosition, relativeVector
* change affects:
finiteVolume/fields/fvPatchFields/derivedFvPatchFields/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
utilities/mesh/generation/blockMesh/curvedEdges/arcEdge.C
------------------------
coordinateRotation.C
coordinateRotation.H
* coordinateRotation is a tensor with restricted contructors
* null contructor yields the global coordinate system
* a null dictionary corresponds to a null constructor
* new axis/direction constructor with functionality taken from
porousZone and previous coordinateSystem implementation
* allow any combination (e1/e2), (e2/e3), (e3,e1) of local vectors
in addition to the axis/direction specification.
- Added check for co-linear axes.
- Could possibly eliminate the non-orthogonality check.
* allow assigment from dictionary, automatically descend in a
'coordinateRotation' sub-dictionary
* add hook in New() for type 'coordinateRotation' with alias 'axes'
------------------------
coordinateSystem.C
coordinateSystem.H
* remove pure virtual restriction - coordinateSystem can be used
directly and has the same properties as a cartesianCS
* null contructor yields the global coordinate system
* a null dictionary entry corresponds to the global coordinate system
* dictionary constructor w/o name uses the typeName
* the coordinateSystem now has a coordinateRotation,
use coordinateRotation constructors instead of calculating any
rotations ourselves
* allow assigment from dictionary, automatically descend into a
'coordinateSystem' sub-dictionary
This simplifies the addition to other dictionary constructors
(eg, porousZone) without requiring a flat structure or special
parsing within the constructor.
eg,
Foam::porousZone::porousZone(const fvMesh& mesh, Istream& is)
:
mesh_(mesh),
name_(is),
dict_(is),
cellZoneID_(mesh_.cellZones().findZoneID(name_)),
csys_(dict_),
C0_(0),
C1_(0),
D_("D", dimensionSet(0, -2, 0, 0, 0), tensor::zero),
F_("F", dimensionSet(0, -1, 0, 0, 0), tensor::zero)
{
...
}
* could consider eliminating Rtr_ member and just use R_.T() instead
* add hook in New() for type 'coordinateSystem'
------------------------
cartesianCS.C
cartesianCS.H
* eliminate redundant virtual functions
------------------------
cylindricalCS.H
* include coordinateSystem.H and not cartesianCS.H
------------------------
coordinateSystems.C
coordinateSystems.H
* provide a means to access an invariant set of coordinate systems

View File

@ -210,25 +210,19 @@ void Foam::coordinateRotation::operator=(const dictionary& rhs)
);
vector axis1, axis2;
axisOrder order = e3e1;
axisOrder order(e3e1);
if (dict.found("e1") && dict.found("e2"))
if (dict.readIfPresent("e1", axis1) && dict.readIfPresent("e2", axis2))
{
order = e1e2;
dict.lookup("e1") >> axis1;
dict.lookup("e2") >> axis2;
}
else if (dict.found("e2") && dict.found("e3"))
else if (dict.readIfPresent("e2", axis1) && dict.readIfPresent("e3", axis2))
{
order = e2e3;
dict.lookup("e2") >> axis1;
dict.lookup("e3") >> axis2;
}
else if (dict.found("e3") && dict.found("e1"))
else if (dict.readIfPresent("e3", axis1) && dict.readIfPresent("e1", axis2))
{
order = e3e1;
dict.lookup("e3") >> axis1;
dict.lookup("e1") >> axis2;
}
else if (dict.found("axis") || dict.found("direction"))
{

View File

@ -267,14 +267,8 @@ void Foam::coordinateSystem::operator=(const dictionary& rhs)
);
// unspecified origin is (0 0 0)
if (dict.found("origin"))
{
dict.lookup("origin") >> origin_;
}
else
{
origin_ = vector::zero;
}
origin_ = vector::zero;
dict.readIfPresent("origin", origin_);
// specify via coordinateRotation
if (dict.found("coordinateRotation"))

View File

@ -117,10 +117,7 @@ Foam::autoPtr<Foam::coordinateSystem> Foam::coordinateSystem::New
// default type is self
word coordType(typeName_());
if (dict.found("type"))
{
dict.lookup("type") >> coordType;
}
dict.readIfPresent("type", coordType);
// can (must) construct base class directly
if (coordType == typeName_())

View File

@ -34,12 +34,6 @@ Description
#include "scalar.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
//- Not VSMALL - but smaller than SMALL
const Foam::scalar Foam::triangleFuncs::NVSMALL = Foam::sqrt(Foam::VSMALL);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::triangleFuncs::setIntersection
@ -120,7 +114,8 @@ bool Foam::triangleFuncs::intersectAxesBundle
scalar det = v2*u1 - u2*v1;
if (Foam::mag(det) < SMALL)
// Fix for V0:(-31.71428 0 -15.10714) V10:(-1.285715 8.99165e-16 -1.142858) V20:(0 0 -1.678573) i0:0
if (Foam::mag(det)/max(max(mag(V10),mag(V20)),1) < SMALL)
{
// Triangle parallel to dir
return false;
@ -522,7 +517,7 @@ bool Foam::triangleFuncs::classify
bool hit = false;
if (Foam::mag(u1) < NVSMALL)
if (Foam::mag(u1) < ROOTVSMALL)
{
beta = u0/u2;

View File

@ -71,11 +71,6 @@ public:
private:
// Private data
//- Not VSMALL - but smaller than SMALL
static const scalar NVSMALL;
// Private Member Functions
//- Helper function for intersect. Sets pt to be anywhere on the edge

View File

@ -2,7 +2,7 @@
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.openfoam.org |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
@ -12,7 +12,6 @@ FoamFile
class dictionary;
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application oodles;

View File

@ -1,22 +1,19 @@
/*-------------------------------*- C++ -*---------------------------------*\
| ========= |
| \\ / OpenFOAM 1.4.1 |
| \\ / |
| \\ / The Open Source CFD Toolbox |
| \\/ http://www.OpenFOAM.org |
\*-------------------------------------------------------------------------*/
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location system;
object probesDict;
version 2.0;
format ascii;
class dictionary;
object probesDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Fields to be probed. runTime modifiable!
fields
(

View File

@ -161,7 +161,7 @@ void triSurface::printTriangle
const pointField& points
)
{
os
os
<< pre.c_str() << "vertex numbers:"
<< f[0] << ' ' << f[1] << ' ' << f[2] << endl
<< pre.c_str() << "vertex coords :"
@ -274,13 +274,13 @@ void triSurface::checkTriangles(const bool verbose)
"triSurface::checkTriangles(bool verbose)"
) << "triangles share the same vertices:\n"
<< " face 1 :" << faceI << endl;
printTriangle(Warning, " ", f, points());
printTriangle(Warning, " ", f, points());
Warning
<< endl
<< " face 2 :"
<< neighbours[neighbourI] << endl;
printTriangle(Warning, " ", n, points());
printTriangle(Warning, " ", n, points());
}
break;
@ -375,8 +375,8 @@ boolList triSurface::checkOrientation(const bool verbose)
(
"triSurface::checkOrientation(bool)"
) << "edge number " << edgeLabels[i] << " on face " << facei
<< " out of range"
<< "\nThis usually means that the input surface has "
<< " out-of-range\n"
<< "This usually means that the input surface has "
<< "edges with more than 2 triangles connected.\n"
<< endl;
valid = false;
@ -413,63 +413,39 @@ boolList triSurface::checkOrientation(const bool verbose)
const labelListList& eFaces = edgeFaces();
// Storage for holding status of edge. True if normal flips across this
// edge
boolList borderEdge(nEdges(), false);
forAll(es, edgei)
forAll(es, edgeI)
{
const labelList& neighbours = eFaces[edgei];
const edge& e = es[edgeI];
const labelList& neighbours = eFaces[edgeI];
if (neighbours.size() == 2)
{
// Two triangles, A and B. Check if edge orientation is
// anticlockwise on both.
const labelList& fEdgesA = faceEdges()[neighbours[0]];
const labelledTri& faceA = (*this)[neighbours[0]];
const labelledTri& faceB = (*this)[neighbours[1]];
// Get next edge after edgei
label nextEdgeA = fEdgesA.fcIndex(findIndex(fEdgesA, edgei));
const labelList& fEdgesB = faceEdges()[neighbours[1]];
label nextEdgeB = fEdgesB.fcIndex(findIndex(fEdgesB, edgei));
// Now check if nextEdgeA and nextEdgeB have any common points
// The edge cannot be going in the same direction if both faces
// are oriented counterclockwise.
// Thus the next face point *must* different between the faces.
if
(
(es[nextEdgeA].start() == es[nextEdgeB].start())
|| (es[nextEdgeA].start() == es[nextEdgeB].end())
|| (es[nextEdgeA].end() == es[nextEdgeB].start())
|| (es[nextEdgeA].end() == es[nextEdgeB].end())
faceA[faceA.fcIndex(findIndex(faceA, e.start()))]
== faceB[faceB.fcIndex(findIndex(faceB, e.start()))]
)
{
borderEdge[edgei] = true;
borderEdge[edgeI] = true;
if (verbose)
{
WarningIn("triSurface::checkOrientation(bool)")
<< "Triangle orientation incorrect." << endl
<< "edge neighbours:" << neighbours << endl
<< "triangle " << neighbours[0] << " has edges "
<< fEdgesA << endl
<< " with points " << endl
<< " " << es[fEdgesA[0]].start() << ' '
<< es[fEdgesA[0]].end() << endl
<< " " << es[fEdgesA[1]].start() << ' '
<< es[fEdgesA[1]].end() << endl
<< " " << es[fEdgesA[2]].start() << ' '
<< es[fEdgesA[2]].end() << endl
<< "triangle " << neighbours[1] << " has edges "
<< fEdgesB << endl
<< " with points " << endl
<< " " << es[fEdgesB[0]].start() << ' '
<< es[fEdgesB[0]].end() << endl
<< " " << es[fEdgesB[1]].start() << ' '
<< es[fEdgesB[1]].end() << endl
<< " " << es[fEdgesB[2]].start() << ' '
<< es[fEdgesB[2]].end() << endl
<< endl;
WarningIn("PrimitivePatchExtra::checkOrientation(bool)")
<< "face orientation incorrect." << nl
<< "edge[" << edgeI << "] " << e
<< " between faces " << neighbours << ":" << nl
<< "face[" << neighbours[0] << "] " << faceA << nl
<< "face[" << neighbours[1] << "] " << faceB << endl;
}
}
}
@ -477,15 +453,14 @@ boolList triSurface::checkOrientation(const bool verbose)
{
if (verbose)
{
const edge& e = es[edgei];
WarningIn("triSurface::checkOrientation(bool)")
<< "Wrong number of edge neighbours." << endl
<< "Edge:" << e
<< "edge[" << edgeI << "] " << e
<< "with points:" << localPoints()[e.start()]
<< ' ' << localPoints()[e.end()]
<< " has neighbours:" << neighbours << endl;
}
borderEdge[edgei] = true;
borderEdge[edgeI] = true;
}
}
@ -1109,7 +1084,7 @@ void triSurface::subsetMeshMap
if (!pointHad[a])
{
pointHad[a] = true;
pointMap[pointI++] = a;
pointMap[pointI++] = a;
}
label b = tri[1];
@ -1147,7 +1122,7 @@ triSurface triSurface::subsetMesh
// Fill pointMap, faceMap
subsetMeshMap(include, pointMap, faceMap);
// Create compact coordinate list and forward mapping array
pointField newPoints(pointMap.size());