Commit 097bba7b authored by Franjo's avatar Franjo

Merge branch 'defect-release-v1.1' into developmentPublicRepo

parents c4e209ff 4c63d7cb
......@@ -32,6 +32,8 @@ Description
#include "gzstream.h"
#include "triSurface.H"
#include "helperFunctions.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
......@@ -163,6 +165,192 @@ void triSurf::writeToFMS(const fileName& fName) const
fStream << subsets;
}
void triSurf::topologyCheck()
{
const pointField& pts = this->points();
const LongList<labelledTri>& trias = this->facets();
//- check for inf and nan points
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 100)
# endif
forAll(pts, pointI)
{
const point& p = pts[pointI];
if( help::isnan(p) || help::isinf(p) )
{
# ifdef USE_OMP
# pragma omp critical
# endif
{
FatalErrorIn
(
"void triSurf::topologyCheck()"
) << "Point " << pointI << " has invalid coordinates "
<< p << exit(FatalError);
}
}
}
//- check whether the nodes are within the scope
//- report duplicate nodes in the same triangle
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 100)
# endif
forAll(trias, triI)
{
const labelledTri& ltri = trias[triI];
forAll(ltri, pI)
{
if( ltri[pI] < 0 || (ltri[pI] >= pts.size()) )
{
# ifdef USE_OMP
# pragma omp critical
# endif
FatalErrorIn
(
"void triSurf::topologyCheck()"
) << "Point " << ltri[pI] << " in triangle " << ltri
<< " is out of scope 0 " << pts.size() << exit(FatalError);
}
if( ltri[pI] == ltri[(pI+1)%3] || ltri[pI] == ltri[(pI+2)%3] )
{
# ifdef USE_OMP
# pragma omp critical
# endif
WarningIn
(
"void triSurf::topologyCheck()"
) << "Triangle " << ltri << " has duplicated points. "
<< "This may cause problems in the meshing process!" << endl;
}
}
}
//- check feature edges
const edgeLongList& featureEdges = this->featureEdges();
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 100)
# endif
forAll(featureEdges, eI)
{
const edge& fe = featureEdges[eI];
forAll(fe, pI)
{
if( fe[pI] < 0 || (fe[pI] >= pts.size()) )
{
# ifdef USE_OMP
# pragma omp critical
# endif
FatalErrorIn
(
"void triSurf::topologyCheck()"
) << "Feature edge " << fe << " point " << fe[pI]
<< " is out of scope 0 " << pts.size() << exit(FatalError);
}
}
if( fe.start() == fe.end() )
{
# ifdef USE_OMP
# pragma omp critical
# endif
WarningIn
(
"void triSurf::topologyCheck()"
) << "Feature edge " << fe << " has duplicated points. "
<< "This may cause problems in the meshing process!" << endl;
}
}
//- check point subsets
DynList<label> subsetIds;
this->pointSubsetIndices(subsetIds);
forAll(subsetIds, i)
{
labelLongList elmts;
this->pointsInSubset(subsetIds[i], elmts);
forAll(elmts, elmtI)
{
const label elI = elmts[elmtI];
if( elI < 0 || elI >= pts.size() )
{
# ifdef USE_OMP
# pragma omp critical
# endif
FatalErrorIn
(
"void triSurf::topologyCheck()"
) << "Point " << elI << " in point subset "
<< this->pointSubsetName(subsetIds[i])
<< " is out of scope 0 " << pts.size() << exit(FatalError);
}
}
}
//- check face subsets
subsetIds.clear();
this->facetSubsetIndices(subsetIds);
forAll(subsetIds, i)
{
labelLongList elmts;
this->facetsInSubset(subsetIds[i], elmts);
forAll(elmts, elmtI)
{
const label elI = elmts[elmtI];
if( elI < 0 || elI >= trias.size() )
{
# ifdef USE_OMP
# pragma omp critical
# endif
FatalErrorIn
(
"void triSurf::topologyCheck()"
) << "Triangle " << elI << " in facet subset "
<< this->facetSubsetName(subsetIds[i])
<< " is out of scope 0 " << trias.size() << exit(FatalError);
}
}
}
//- check feature edge subsets
subsetIds.clear();
this->edgeSubsetIndices(subsetIds);
forAll(subsetIds, i)
{
labelLongList elmts;
this->edgesInSubset(subsetIds[i], elmts);
forAll(elmts, elmtI)
{
const label elI = elmts[elmtI];
if( elI < 0 || elI >= featureEdges.size() )
{
# ifdef USE_OMP
# pragma omp critical
# endif
FatalErrorIn
(
"void triSurf::topologyCheck()"
) << "Feature edge " << elI << " in edge subset "
<< this->edgeSubsetName(subsetIds[i])
<< " is out of scope 0 " << featureEdges.size()
<< exit(FatalError);
}
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
triSurf::triSurf()
......@@ -186,7 +374,9 @@ triSurf::triSurf
triSurfFacets(triangles, patches),
triSurfFeatureEdges(featureEdges),
triSurfAddressing(triSurfPoints::points_, triSurfFacets::triangles_)
{}
{
topologyCheck();
}
//- Read from file
triSurf::triSurf(const fileName& fName)
......@@ -197,6 +387,8 @@ triSurf::triSurf(const fileName& fName)
triSurfAddressing(triSurfPoints::points_, triSurfFacets::triangles_)
{
readSurface(fName);
topologyCheck();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
......
......@@ -74,6 +74,8 @@ class triSurf
inline LongList<labelledTri>& accessToFacets();
inline geometricSurfacePatchList& accessToPatches();
void topologyCheck();
//- Disallow default bitwise assignment
void operator=(const triSurf&);
......
......@@ -29,6 +29,8 @@ Description
#include "VRWGraphSMPModifier.H"
#include "demandDrivenData.H"
#include <set>
# ifdef USE_OMP
#include <omp.h>
# endif
......@@ -70,49 +72,42 @@ void triSurfAddressing::calculateEdges() const
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(facets_, fI)
forAll(pFacets, pI)
{
const labelledTri& tri = facets_[fI];
std::set<std::pair<label, label> > edgesAtPoint;
forAll(tri, pI)
forAllRow(pFacets, pI, pfI)
{
const edge fe(tri[pI], tri[(pI+1)%3]);
const label s = fe.start();
const label e = fe.end();
DynList<label> edgeFaces;
bool store(true);
const label triI = pFacets(pI, pfI);
const labelledTri& tri = facets_[triI];
//- find all faces attached to this edge
//- store the edge in case the face faceI is the face
//- with the smallest label
forAllRow(pFacets, s, pfI)
label pos(-1);
forAll(tri, i)
{
const label ofI = pFacets(s, pfI);
const labelledTri& of = facets_[ofI];
label pos(-1);
for(label i=0;i<3;++i)
if( of[i] == e )
{
pos = i;
break;
}
if( pos < 0 )
continue;
if( ofI < fI )
if( tri[i] == pI )
{
store = false;
break;
pos = i;
if( tri[(pos+1)%3] >= pI )
edgesAtPoint.insert
(
std::make_pair(pI, tri[(pos+1)%3])
);
if( tri[(pos+2)%3] >= pI )
edgesAtPoint.insert
(
std::make_pair(pI, tri[(pos+2)%3])
);
}
edgeFaces.append(ofI);
}
if( store )
edgesHelper.append(fe);
}
std::set<std::pair<label, label> >::const_iterator it;
for(it=edgesAtPoint.begin();it!=edgesAtPoint.end();++it)
edgesHelper.append(edge(it->first, it->second));
}
//- this enables other threads to see the number of edges
......@@ -158,56 +153,58 @@ void triSurfAddressing::calculateFacetEdges() const
const edgeLongList& edges = this->edges();
const VRWGraph& pointFaces = this->pointFacets();
facetEdgesPtr_ = new VRWGraph(facets_.size());
facetEdgesPtr_ = new VRWGraph(facets_.size(), 3, -1);
VRWGraph& faceEdges = *facetEdgesPtr_;
labelList nfe(facets_.size());
# ifdef USE_OMP
const label nThreads = 3 * omp_get_num_procs();
# pragma omp parallel num_threads(nThreads)
# pragma omp parallel for schedule(dynamic, 100)
# endif
forAll(edges, edgeI)
{
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(facets_, fI)
nfe[fI] = 3;
# ifdef USE_OMP
# pragma omp barrier
# pragma omp master
# endif
VRWGraphSMPModifier(faceEdges).setSizeAndRowSize(nfe);
# ifdef USE_OMP
# pragma omp barrier
const edge ee = edges[edgeI];
const label pI = ee.start();
# pragma omp for schedule(static)
# endif
forAll(edges, edgeI)
forAllRow(pointFaces, pI, pfI)
{
const edge ee = edges[edgeI];
const label pI = ee.start();
const label fI = pointFaces(pI, pfI);
forAllRow(pointFaces, pI, pfI)
const labelledTri& tri = facets_[fI];
forAll(tri, eI)
{
const label fI = pointFaces(pI, pfI);
const edge e(tri[eI], tri[(eI+1)%3]);
const labelledTri& tri = facets_[fI];
forAll(tri, eI)
if( e == ee )
{
const edge e(tri[eI], tri[(eI+1)%3]);
if( e == ee )
{
faceEdges(fI, eI) = edgeI;
break;
}
faceEdges(fI, eI) = edgeI;
}
}
}
}
# ifdef DEBUGTriSurfAddressing
forAll(faceEdges, triI)
{
forAllRow(faceEdges, triI, feI)
{
if( facets_[triI][feI] < 0 || facets_[triI][feI] >= points_.size() )
FatalErrorIn
(
"void triSurfAddressing::calculateFacetEdges() const"
) << "Invalid entry in triangle " << triI
<< " " << facets_[triI] << abort(FatalError);
const label edgeI = faceEdges(triI, feI);
if( edgeI < 0 || edgeI >= edges.size() )
FatalErrorIn
(
"void triSurfAddressing::calculateFacetEdges() const"
) << "Invalid entry in face " << triI << " "
<< facets_[triI] << " edges "
<< faceEdges[triI] << abort(FatalError);
}
}
# endif
}
void triSurfAddressing::calculateEdgeFacets() const
......
......@@ -39,6 +39,10 @@ Description
//#define DEBUGVrt
# ifdef DEBUGVrt
#include "OFstream.H"
# endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
......@@ -46,6 +50,7 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
# ifdef DEBUGVrt
void writeOctreeToVTK
(
const fileName& fName,
......@@ -113,6 +118,7 @@ void writeOctreeToVTK
file << nl;
}
# endif
void meshOctreeAddressing::createOctreePoints() const
{
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment