Commit 060cbdf9 authored by franjo.juretic@c-fields.com's avatar franjo.juretic@c-fields.com
Browse files

Improvements on feature capturing


git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@60 fdcce57e-7e00-11e2-b579-49867b4cea03
parent b7cab04b
......@@ -6,5 +6,6 @@
EXE_INC = \
$(OMP_FLAGS) \
-g -ggdb -DFULLDEBUG -O0 \
-I$(LIB_SRC)/triSurface/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2005-2007 Franjo Juretic
\\/ 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Class
refLabelledPointScalar
Description
A class containing a label and labelledPointScalar. It is used for
exchanging data over processors
SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef refLabelledPointScalar_H
#define refLabelledPointScalar_H
#include "labelledPointScalar.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class refLabelledPointScalar Declaration
\*---------------------------------------------------------------------------*/
class refLabelledPointScalar
{
// Private data
//- point label
label pLabel_;
//- labelledPointScalar
labelledPointScalar lps_;
public:
// Constructors
//- Null construct
refLabelledPointScalar()
:
pLabel_(-1),
lps_()
{}
//- Construct from label and labelledPointScalar
refLabelledPointScalar(const label pl, const labelledPointScalar& lps)
:
pLabel_(pl),
lps_(lps)
{}
// Destructor
~refLabelledPointScalar()
{}
// Member functions
//- return object label
inline label objectLabel() const
{
return pLabel_;
}
//- return labelledPointScalar
inline const labelledPointScalar& lps() const
{
return lps_;
}
// Member operators
inline void operator=(const refLabelledPointScalar& lps)
{
pLabel_ = lps.pLabel_;
lps_ = lps.lps_;
}
inline bool operator==(const refLabelledPointScalar& lps) const
{
if( pLabel_ == lps.pLabel_ )
return true;
return false;
}
inline bool operator!=(const refLabelledPointScalar& lps) const
{
return !this->operator==(lps);
}
// Friend operators
friend Ostream& operator<<(Ostream& os, const refLabelledPointScalar& lps)
{
os << token::BEGIN_LIST;
os << lps.pLabel_ << token::SPACE;
os << lps.lps_ << token::END_LIST;
// Check state of Ostream
os.check("operator<<(Ostream&, const labelledPointScalarS&");
return os;
}
friend Istream& operator>>(Istream& is, refLabelledPointScalar& lps)
{
// Read beginning of refLabelledPointScalar
is.readBegin("refLabelledPointScalar");
is >> lps.pLabel_;
is >> lps.lps_;
// Read end of refLabelledPointScalar
is.readEnd("refLabelledPointScalar");
// Check state of Istream
is.check("operator>>(Istream&, refLabelledPointScalar");
return is;
}
};
//- Specify data associated with refLabelledPointScalar type is contiguous
template<>
inline bool contiguous<refLabelledPointScalar>() {return true;}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
......@@ -211,9 +211,10 @@ void triSurfAddressing::calculateFacetEdges() const
void triSurfAddressing::calculateEdgeFacets() const
{
const edgeLongList& edges = this->edges();
const VRWGraph& faceEdges = this->facetEdges();
edgeFacetsPtr_ = new VRWGraph();
edgeFacetsPtr_ = new VRWGraph(edges.size());
VRWGraphSMPModifier(*edgeFacetsPtr_).reverseAddressing(faceEdges);
}
......@@ -222,9 +223,9 @@ void triSurfAddressing::calculatePointEdges() const
{
const edgeLongList& edges = this->edges();
pointEdgesPtr_ = new VRWGraph();
pointEdgesPtr_ = new VRWGraph(points_.size());
VRWGraphSMPModifier(*pointEdgesPtr_).reverseAddressing(edges);
pointEdgesPtr_->reverseAddressing(edges);
}
void triSurfAddressing::calculateFacetFacetsEdges() const
......
......@@ -229,9 +229,9 @@ bool meshOctreeAutomaticRefinement::refineBasedOnContainedPartitions
const triSurfacePartitioner& sPart = this->partitioner();
//- find leaves which contains corner nodes
const List<labelHashSet>& pPart = sPart.partitionPartitions();
const labelList& edgePartition = sPart.edgePartitions();
const List<labelHashSet>& ePart = sPart.edgePartitionEdgePartitions();
const List<labelHashSet>& pPatches = sPart.patchPatches();
const labelList& edgeGroups = sPart.edgeGroups();
const List<labelHashSet>& eNeiGroups = sPart.edgeGroupEdgeGroups();
# ifdef DEBUGAutoRef
Info << "pPart " << pPart << endl;
......@@ -243,10 +243,10 @@ bool meshOctreeAutomaticRefinement::refineBasedOnContainedPartitions
const triSurf& surf = octree_.surface();
const LongList<meshOctreeCube*>& leaves = octreeModifier.leavesAccess();
DynList<label> patches, ePartitions, helper;
DynList<label> patches, eGroups, helper;
# ifdef USE_OMP
# pragma omp parallel for if( refCandidates.size() > 1000 ) \
private(patches, ePartitions, helper) \
private(patches, eGroups, helper) \
reduction(+ : nMarked) schedule(dynamic, 20)
# endif
forAll(refCandidates, refI)
......@@ -271,9 +271,9 @@ bool meshOctreeAutomaticRefinement::refineBasedOnContainedPartitions
//- find edge partitions contained in this box
helper.clear();
octree_.findEdgesInBox(bb, helper);
ePartitions.clear();
eGroups.clear();
forAll(helper, i)
ePartitions.appendIfNotIn(edgePartition[helper[i]]);
eGroups.appendIfNotIn(edgeGroups[helper[i]]);
# ifdef DEBUGAutoRef
Info << "patches for leaf " << leafI << " are " << patches << endl;
......@@ -283,7 +283,7 @@ bool meshOctreeAutomaticRefinement::refineBasedOnContainedPartitions
forAll(patches, patchI)
{
for(label patchJ=(patchI+1);patchJ<patches.size();++patchJ)
if( !pPart[patches[patchI]].found(patches[patchJ]) )
if( !pPatches[patches[patchI]].found(patches[patchJ]) )
{
# ifdef DEBUGAutoRef
Info << "2.Here" << endl;
......@@ -294,10 +294,10 @@ bool meshOctreeAutomaticRefinement::refineBasedOnContainedPartitions
}
}
forAll(ePartitions, ePartI)
forAll(eGroups, egI)
{
for(label ePartJ=ePartI+1;ePartJ<ePartitions.size();++ePartJ)
if( !ePart[ePartitions[ePartI]].found(ePartitions[ePartJ]) )
for(label egJ=egI+1;egJ<eGroups.size();++egJ)
if( !eNeiGroups[eGroups[egI]].found(eGroups[egJ]) )
{
refine = true;
break;
......
......@@ -933,7 +933,7 @@ bool edgeExtractor::distributeBoundaryFacesNormalAlignment()
const triSurf& surf = meshOctree_.surface();
const pointField& sPoints = surf.points();
label nCorrected;
label nCorrected, nIterations(0);
Map<label> otherProcNewPatch;
do
......@@ -1068,7 +1068,7 @@ bool edgeExtractor::distributeBoundaryFacesNormalAlignment()
//- transfer the new patches back
facePatch_.transfer(newBoundaryPatches);
}
} while( nCorrected != 0 );
} while( (nCorrected != 0) && (++nIterations < 5) );
return changed;
}
......@@ -1207,9 +1207,9 @@ void edgeExtractor::findEdgeCandidates()
}
//- start post-processing gathered data
const labelList& edgePartition = partitioner.edgePartitions();
const labelList& edgeGroup = partitioner.edgeGroups();
List<List<labelledScalar> > edgePartitionAndWeights(edges.size());
List<List<labelledScalar> > edgeGroupAndWeights(edges.size());
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 40) \
......@@ -1231,11 +1231,11 @@ void edgeExtractor::findEdgeCandidates()
DynList<labelledScalar> weights;
forAll(sc, i)
{
const label sPart = edgePartition[sc[i].scalarLabel()];
const label sPart = edgeGroup[sc[i].scalarLabel()];
forAll(ec, j)
{
const label ePart = edgePartition[ec[j].scalarLabel()];
const label ePart = edgeGroup[ec[j].scalarLabel()];
if( (sPart >= 0) && (sPart == ePart) )
{
......@@ -1252,16 +1252,16 @@ void edgeExtractor::findEdgeCandidates()
}
//- store the data
edgePartitionAndWeights[edgeI].setSize(weights.size());
forAll(edgePartitionAndWeights[edgeI], epI)
edgePartitionAndWeights[edgeI][epI] = weights[epI];
edgeGroupAndWeights[edgeI].setSize(weights.size());
forAll(edgeGroupAndWeights[edgeI], epI)
edgeGroupAndWeights[edgeI][epI] = weights[epI];
//- sort the data according to the weights
stableSort(edgePartitionAndWeights[edgeI]);
stableSort(edgeGroupAndWeights[edgeI]);
}
}
Info << "Edge partitions and weights " << edgePartitionAndWeights << endl;
Info << "Edge partitions and weights " << edgeGroupAndWeights << endl;
}
bool edgeExtractor::checkConcaveEdgeCells()
......@@ -1759,11 +1759,6 @@ bool edgeExtractor::checkFacePatchesTopology()
if( (newPatch < 0) || (newPatch == facePatch_[bfI]) )
continue;
//- do not swap if the number of neighbours in the patch
//- does not dominate
if( nNeiEdges <= (bf.size() / 2) )
continue;
//- check whether the edges shared ith the neighbour patch form
//- a singly linked chain
DynList<bool> sharedEdge;
......@@ -1790,14 +1785,8 @@ bool edgeExtractor::checkFacePatchesTopology()
//- transfer the new patches back
facePatch_.transfer(newBoundaryPatches);
const triSurf* surfPtr = surfaceWithPatches();
surfPtr->writeSurface("checkFacePatches_"+help::scalarToText(nIter)+".stl");
deleteDemandDrivenData(surfPtr);
}
//break;
} while( nCorrected != 0 && (nIter < 3) );
return changed;
......@@ -2773,7 +2762,23 @@ void edgeExtractor::extractEdges()
Info << "No topological adjustment was needed" << endl;
}
return;
Info << "Checking quality of corner points" << endl;
if( checkCorners() )
{
Info << "Finished checking quality of corner points" << endl;
# ifdef DEBUGEdgeExtractor
Info << "Changes due to face patches" << endl;
fileName sName("checkPoints"+help::scalarToText(nIter)+".stl");
sPtr = surfaceWithPatches();
sPtr->writeSurface(sName);
deleteDemandDrivenData(sPtr);
# endif
}
else
{
Info << "No adjustment of corner points was needed" << endl;
}
//- project the selected feature vertices onto the corresponding
//- objects on the surface mesh and improve the distribution of patches
......@@ -2787,7 +2792,7 @@ void edgeExtractor::extractEdges()
//- project feature vertices onto the surface mesh
projectDeterminedFeatureVertices();
if( checkFacePatchesGeometry() )
if( checkCorners() )
{
# ifdef DEBUGEdgeExtractor
Info << "Changes due to untangling" << endl;
......@@ -2856,6 +2861,59 @@ const triSurf* edgeExtractor::surfaceWithPatches() const
return surfPtr;
}
const triSurf* edgeExtractor::surfaceWithPatches(const label bpI) const
{
//- allocate the memory for the surface mesh
triSurf* surfPtr = new triSurf();
//- surface of the volume mesh
const meshSurfaceEngine& mse = surfaceEngine();
const faceList::subList& bFaces = mse.boundaryFaces();
const VRWGraph& pFaces = mse.pointFaces();
const labelList& bp = mse.bp();
const pointFieldPMG& points = mesh_.points();
//- modifier of the new surface mesh
triSurfModifier surfModifier(*surfPtr);
surfModifier.patchesAccess() = meshOctree_.surface().patches();
pointField& sPts = surfModifier.pointsAccess();
//- create the triangulation of the volume mesh surface
labelLongList newPointLabel(points.size(), -1);
label nPoints(0);
forAllRow(pFaces, bpI, pfI)
{
const label bfI = pFaces(bpI, pfI);
const face& bf = bFaces[bfI];
forAll(bf, pI)
if( newPointLabel[bf[pI]] == -1 )
newPointLabel[bf[pI]] = nPoints++;
labelledTri tri;
tri.region() = facePatch_[bfI];
tri[0] = newPointLabel[bf[0]];
for(label i=bf.size()-2;i>0;--i)
{
tri[1] = newPointLabel[bf[i]];
tri[2] = newPointLabel[bf[i+1]];
surfPtr->appendTriangle(tri);
}
}
//- copy points
sPts.setSize(nPoints);
forAll(newPointLabel, pointI)
if( newPointLabel[pointI] != -1 )
{
sPts[newPointLabel[pointI]] = points[pointI];
}
return surfPtr;
}
void edgeExtractor::updateMeshPatches()
{
const triSurf& surface = meshOctree_.surface();
......
......@@ -230,10 +230,14 @@ public:
//- assemble the above functionality into a workflow
void extractEdges();
//- generate a surface mesh based and store the created patches
//- generate a surface mesh and store the created patches
//- this is mainly intended for debugging purposes
const triSurf* surfaceWithPatches() const;
//- generate a surface mesh constin of facets adjacent to the requested
//- surface point
const triSurf* surfaceWithPatches(const label bpI) const;
//- update mesh with selected patches
void updateMeshPatches();
};
......
......@@ -52,6 +52,7 @@ void meshSurfaceEdgeExtractorNonTopo::distributeBoundaryFaces()
{
edgeExtractor extractor(mesh_, meshOctree_);
Info << "Extracting edges" << endl;
extractor.extractEdges();
extractor.updateMeshPatches();
......
......@@ -31,6 +31,7 @@ Description
#include "meshSurfaceMapper.H"
#include "meshOctree.H"
#include "refLabelledPoint.H"
#include "refLabelledPointScalar.H"
#include "helperFunctionsPar.H"
# ifdef USE_OMP
......@@ -54,8 +55,27 @@ void meshSurfaceMapper::preMapVertices(const label nIterations)
const pointFieldPMG& points = surfaceEngine_.points();
const vectorField& faceCentres = surfaceEngine_.faceCentres();
const VRWGraph& pointFaces = surfaceEngine_.pointFaces();
const VRWGraph& pointInFace = surfaceEngine_.pointInFaces();
const faceList::subList& bFaces = surfaceEngine_.boundaryFaces();
List<labelledPoint> preMapPositions(boundaryPoints.size());
List<labelledPointScalar> preMapPositions(boundaryPoints.size());
List<DynList<scalar, 6> > faceCentreDistances(bFaces.size());
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 20)
# endif
forAll(bFaces, bfI)
{
const point& c = faceCentres[bfI];
const face& bf = bFaces[bfI];
faceCentreDistances[bfI].setSize(bf.size());
forAll(bf, pI)
{
faceCentreDistances[bfI][pI] = magSqr(points[bf[pI]] - c);
}
}
for(label iterI=0;iterI<nIterations;++iterI)
{
......@@ -65,12 +85,21 @@ void meshSurfaceMapper::preMapVertices(const label nIterations)
# endif
forAll(pointFaces, bpI)
{
labelledPoint lp(0, vector::zero);
labelledPointScalar lp(bpI, vector::zero, 0.0);
forAllRow(pointFaces, bpI, bfI)
const point& p = points[boundaryPoints[bpI]];
forAllRow(pointFaces, bpI, pfI)
{
++lp.pointLabel();
lp.coordinates() += faceCentres[pointFaces(bpI, bfI)];
const label bfI = pointFaces(bpI, pfI);
const point& fc = faceCentres[pointFaces(bpI, pfI)];
const label pos = pointInFace(bpI, pfI);
const scalar w
(
max(magSqr(p - fc) / faceCentreDistances[bfI][pos], SMALL)
);
lp.coordinates() += w * faceCentres[bfI];
lp.scalarValue() += w;
}
preMapPositions[bpI] = lp;
......@@ -89,14 +118,14 @@ void meshSurfaceMapper::preMapVertices(const label nIterations)
surfaceEngine_.globalToLocalBndPointAddressing();
//- collect data to be sent to other processors
std::map<label, LongList<refLabelledPoint> > exchangeData;
std::map<label, LongList<labelledPointScalar> > exchangeData;
forAll(surfaceEngine_.bpNeiProcs(), i)
exchangeData.insert
(
std::make_pair
(
surfaceEngine_.bpNeiProcs()[i],
LongList<refLabelledPoint>()
LongList<labelledPointScalar>()
)
);
......@@ -113,30 +142,30 @@ void meshSurfaceMapper::preMapVertices(const label nIterations)
exchangeData[neiProc].append
(
refLabelledPoint
labelledPointScalar
(
globalPointLabel[bpI],
preMapPositions[bpI]
preMapPositions[bpI].coordinates(),
preMapPositions[bpI].scalarValue()
)
);
}
}
//- exchange data with other processors
LongList<refLabelledPoint> receivedData;
LongList<labelledPointScalar> receivedData;
help::exchangeMap(exchangeData, receivedData);
//- combine collected data with the available data
forAll(receivedData, i)
{
const refLabelledPoint& rlp = receivedData[i];
const labelledPoint& lps = rlp.lPoint();
const labelledPointScalar& lps = receivedData[i];
const label bpI = globalToLocal[rlp.objectLabel()];
const label bpI = globalToLocal[lps.pointLabel()];
labelledPoint& lp = preMapPositions[bpI];