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 'cvm' of ssh://hunt//home/noisy3/OpenFOAM/OpenFOAM-dev into cvm

Browse Source
This commit is contained in:
graham
2011-06-20 21:28:06 +01:00
parent 6fc2e685f4 3d3175930d
commit 9ccf31bc65
17 changed files with 1091 additions and 245 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
applications/utilities/mesh/generation/cvMesh/cvMeshDict
View File

@ -210,16 +210,15 @@ initialPoints
// to the surface. Is fraction of local target cell size (see below)
minimumSurfaceDistanceCoeff 0.55;
initialPointsMethod hierarchicalDensityWeightedStochastic;
initialPointsMethod autoDensity;
// initialPointsMethod uniformGrid;
// initialPointsMethod bodyCentredCubic;
// initialPointsMethod pointFile;
// initialPointsMethod densityWeightedStochastic;
// Take boundbox of all geometry. Sample with this box. If too much
// samples in box (due to target cell size) split box.
hierarchicalDensityWeightedStochasticDetails
autoDensityDetails
{
// Initial number of refinement levels. Needs to be enough to pick
// up features due to size ratio. If not enough it will take longer
@ -233,12 +232,6 @@ initialPoints
surfaceSampleResolution 5;
}
densityWeightedStochasticDetails
{
// Volume of space to be meshed. Use e.g. surfaceInertia to find this.
totalVolume 1.56e-05;
}
uniformGridDetails
{
// Absolute cell size.

1
applications/utilities/mesh/manipulation/topoSet/Make/files
View File

@ -1,3 +1,4 @@
timeSelector.C
topoSet.C
EXE = $(FOAM_APPBIN)/topoSet

289
applications/utilities/mesh/manipulation/topoSet/timeSelector.C Normal file
View File

@ -0,0 +1,289 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2011 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 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 "timeSelector.H"
#include "ListOps.H"
#include "argList.H"
#include "Time.H"
#include "IStringStream.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::timeSelector::timeSelector()
:
scalarRanges()
{}
Foam::timeSelector::timeSelector(Istream& is)
:
scalarRanges(is)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::timeSelector::selected(const instant& value) const
{
return scalarRanges::selected(value.value());
}
Foam::List<bool> Foam::timeSelector::selected(const instantList& Times) const
{
List<bool> lst(Times.size(), false);
// check ranges, avoid false positive on constant/
forAll(Times, timeI)
{
if (Times[timeI].name() != "constant" && selected(Times[timeI]))
{
lst[timeI] = true;
}
}
// check specific values
forAll(*this, rangeI)
{
if (operator[](rangeI).isExact())
{
scalar target = operator[](rangeI).value();
int nearestIndex = -1;
scalar nearestDiff = Foam::GREAT;
forAll(Times, timeI)
{
if (Times[timeI].name() == "constant") continue;
scalar diff = fabs(Times[timeI].value() - target);
if (diff < nearestDiff)
{
nearestDiff = diff;
nearestIndex = timeI;
}
}
if (nearestIndex >= 0)
{
lst[nearestIndex] = true;
}
}
}
return lst;
}
Foam::List<Foam::instant> Foam::timeSelector::select
(
const instantList& Times
) const
{
return subset(selected(Times), Times);
}
void Foam::timeSelector::inplaceSelect(instantList& Times) const
{
inplaceSubset(selected(Times), Times);
}
void Foam::timeSelector::addOptions
(
const bool constant,
const bool zeroTime
)
{
if (constant)
{
argList::addBoolOption
(
"constant",
"include the 'constant/' dir in the times list"
);
}
if (zeroTime)
{
argList::addBoolOption
(
"zeroTime",
"include the '0/' dir in the times list"
);
}
argList::addBoolOption
(
"noZero",
"exclude the '0/' dir from the times list, "
"has precedence over the -zeroTime option"
);
argList::addBoolOption
(
"latestTime",
"select the latest time"
);
argList::addOption
(
"time",
"ranges",
"comma-separated time ranges - eg, ':10,20,40-70,1000:'"
);
}
Foam::List<Foam::instant> Foam::timeSelector::select
(
const instantList& timeDirs,
const argList& args
)
{
if (timeDirs.size())
{
List<bool> selectTimes(timeDirs.size(), true);
// determine locations of constant/ and 0/ directories
label constantIdx = -1;
label zeroIdx = -1;
forAll(timeDirs, timeI)
{
if (timeDirs[timeI].name() == "constant")
{
constantIdx = timeI;
}
else if (timeDirs[timeI].value() == 0)
{
zeroIdx = timeI;
}
if (constantIdx >= 0 && zeroIdx >= 0)
{
break;
}
}
// determine latestTime selection (if any)
// this must appear before the -time option processing
label latestIdx = -1;
if (args.optionFound("latestTime"))
{
selectTimes = false;
latestIdx = timeDirs.size() - 1;
// avoid false match on constant/
if (latestIdx == constantIdx)
{
latestIdx = -1;
}
}
if (args.optionFound("time"))
{
// can match 0/, but can never match constant/
selectTimes = timeSelector
(
args.optionLookup("time")()
).selected(timeDirs);
}
// add in latestTime (if selected)
if (latestIdx >= 0)
{
selectTimes[latestIdx] = true;
}
if (constantIdx >= 0)
{
// only add constant/ if specifically requested
selectTimes[constantIdx] = args.optionFound("constant");
}
// special treatment for 0/
if (zeroIdx >= 0)
{
if (args.optionFound("noZero"))
{
// exclude 0/ if specifically requested
selectTimes[zeroIdx] = false;
}
else if (argList::validOptions.found("zeroTime"))
{
// with -zeroTime enabled, drop 0/ unless specifically requested
selectTimes[zeroIdx] = args.optionFound("zeroTime");
}
}
return subset(selectTimes, timeDirs);
}
else
{
return timeDirs;
}
}
Foam::List<Foam::instant> Foam::timeSelector::select0
(
Time& runTime,
const argList& args,
const bool useOptionsOnly
)
{
if
(
useOptionsOnly
|| (
args.optionFound("latestTime")
|| args.optionFound("time")
|| args.optionFound("constant")
|| args.optionFound("noZero")
|| args.optionFound("zeroTime")
)
)
{
instantList timeDirs = timeSelector::select(runTime.times(), args);
if (timeDirs.empty())
{
FatalErrorIn(args.executable())
<< "No times selected"
<< exit(FatalError);
}
runTime.setTime(timeDirs[0], 0);
return timeDirs;
}
else
{
// No timeSelector option specified. Do not change runTime.
return instantList(1, instant(runTime.value(), runTime.timeName()));
}
}
// ************************************************************************* //

163
applications/utilities/mesh/manipulation/topoSet/timeSelector.H Normal file
View File

@ -0,0 +1,163 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2011 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 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::timeSelector
Description
A List of scalarRange for selecting times.
The timeSelector provides a convenient means of selecting multiple
times. A typical use would be the following:
\verbatim
timeSelector::addOptions();
// add other options
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
...
forAll(timeDirs, timeI)
{
...
}
\endverbatim
The result program would receive \b -time, @b -latestTime, @b -constant
and \b -noZero options. The @b -constant option explicitly includes the
\c constant/ directory in the time list and the \b -noZero option
explicitly excludes the \c 0/ directory from the time list.
There may however also be many cases in which neither the \c constant/
directory nor the \c 0/ directory contain particularly relevant
information. This might occur, for example, when post-processing
results. In this case, addOptions is called with optional boolean
arguments.
\verbatim
timeSelector::addOptions(false, true);
\endverbatim
The first argument avoids adding the \b -constant option. The second
argument adds an additional \b -zeroTime option and also prevents the
\c 0/ directory from being included in the default time range and in the
\b -latestTime selection.
SourceFiles
timeSelector.C
\*---------------------------------------------------------------------------*/
#ifndef timeSelector_H
#define timeSelector_H
#include "scalarRanges.H"
#include "instantList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class argList;
class Time;
/*---------------------------------------------------------------------------*\
Class timeSelector Declaration
\*---------------------------------------------------------------------------*/
class timeSelector
:
public scalarRanges
{
public:
// Constructors
//- Construct Null
timeSelector();
//- Construct from Istream
timeSelector(Istream&);
// Member Functions
//- Return true if the given instant is within the ranges
bool selected(const instant&) const;
//- Return the set of selected instants in the given list that are
// within the ranges
List<bool> selected(const List<instant>&) const;
//- Select a list of Time values that are within the ranges
instantList select(const List<instant>&) const;
//- Select a list of Time values that are within the ranges
void inplaceSelect(List<instant>&) const;
//- Add the options handled by timeSelector to argList::validOptions
//
// \param constant
// Add the \b -constant option to include the \c constant/ directory
//
// \param zeroTime
// Enable the \b -zeroTime option and alter the normal time selection
// behaviour (and \b -latestTime behaviour) to exclude the \c 0/
// directory. The \c 0/ directory will only be included when
// \b -zeroTime is specified.
// The \b -noZero option has precedence over the @b -zeroTime option.
static void addOptions
(
const bool constant=true,
const bool zeroTime=false
);
//- Return the set of times selected based on the argList options
static instantList select
(
const instantList&,
const argList& args
);
//- Return the set of times selected based on the argList options
// also set the runTime to the first instance
static instantList select0
(
Time& runTime,
const argList& args,
const bool useOptionsOnly
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

284
applications/utilities/mesh/manipulation/topoSet/topoSet.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,6 +26,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "timeSelector.H"
#include "argList.H"
#include "Time.H"
#include "polyMesh.H"
@ -33,15 +34,71 @@ Description
#include "cellSet.H"
#include "faceSet.H"
#include "pointSet.H"
#include "globalMeshData.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void printMesh(const Time& runTime, const polyMesh& mesh)
{
Info<< "Time:" << runTime.timeName()
<< " cells:" << mesh.globalData().nTotalCells()
<< " faces:" << mesh.globalData().nTotalFaces()
<< " points:" << mesh.globalData().nTotalPoints()
<< " patches:" << mesh.boundaryMesh().size()
<< " bb:" << mesh.bounds() << nl;
}
polyMesh::readUpdateState meshReadUpdate(polyMesh& mesh)
{
polyMesh::readUpdateState stat = mesh.readUpdate();
switch(stat)
{
case polyMesh::UNCHANGED:
{
Info<< " mesh not changed." << endl;
break;
}
case polyMesh::POINTS_MOVED:
{
Info<< " points moved; topology unchanged." << endl;
break;
}
case polyMesh::TOPO_CHANGE:
{
Info<< " topology changed; patches unchanged." << nl
<< " ";
printMesh(mesh.time(), mesh);
break;
}
case polyMesh::TOPO_PATCH_CHANGE:
{
Info<< " topology changed and patches changed." << nl
<< " ";
printMesh(mesh.time(), mesh);
break;
}
default:
{
FatalErrorIn("meshReadUpdate(polyMesh&)")
<< "Illegal mesh update state "
<< stat << abort(FatalError);
break;
}
}
return stat;
}
// Main program:
int main(int argc, char *argv[])
{
timeSelector::addOptions(true, false);
argList::addOption
(
"dict",
@ -57,6 +114,14 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
instantList timeDirs = timeSelector::select0
(
runTime,
args,
false // not override runTime if no time options
);
# include "createNamedPolyMesh.H"
const bool noSync = args.optionFound("noSync");
@ -99,128 +164,143 @@ int main(int argc, char *argv[])
// Read set construct info from dictionary
PtrList<dictionary> patchSources(topoSetDict.lookup("actions"));
PtrList<dictionary> actions(topoSetDict.lookup("actions"));
forAll(patchSources, i)
forAll(timeDirs, timeI)
{
const dictionary& dict = patchSources[i];
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Time = " << runTime.timeName() << endl;
const word setName(dict.lookup("name"));
const word actionName(dict.lookup("action"));
const word setType(dict.lookup("type"));
// Optionally re-read mesh
meshReadUpdate(mesh);
topoSetSource::setAction action = topoSetSource::toAction(actionName);
autoPtr<topoSet> currentSet;
if
(
(action == topoSetSource::NEW)
|| (action == topoSetSource::CLEAR)
)
// Execute all actions
forAll(actions, i)
{
currentSet = topoSet::New(setType, mesh, setName, 10000);
Info<< "Created set " << setName << endl;
}
else if (action == topoSetSource::REMOVE)
{
//?
}
else
{
currentSet = topoSet::New
const dictionary& dict = actions[i];
const word setName(dict.lookup("name"));
const word actionName(dict.lookup("action"));
const word setType(dict.lookup("type"));
topoSetSource::setAction action = topoSetSource::toAction
(
setType,
mesh,
setName,
IOobject::MUST_READ
actionName
);
Info<< "Read set " << setName << " with size "
<< currentSet().size() << endl;
}
// Handle special actions (clear, invert) locally, rest through sources.
switch (action)
{
case topoSetSource::NEW:
case topoSetSource::ADD:
case topoSetSource::DELETE:
autoPtr<topoSet> currentSet;
if
(
(action == topoSetSource::NEW)
|| (action == topoSetSource::CLEAR)
)
{
Info<< " Applying source " << word(dict.lookup("source"))
<< endl;
autoPtr<topoSetSource> source = topoSetSource::New
(
dict.lookup("source"),
mesh,
dict.subDict("sourceInfo")
);
source().applyToSet(action, currentSet());
// Synchronize for coupled patches.
if (!noSync) currentSet().sync(mesh);
currentSet().write();
currentSet = topoSet::New(setType, mesh, setName, 10000);
Info<< "Created set " << setName << endl;
}
break;
case topoSetSource::SUBSET:
else if (action == topoSetSource::REMOVE)
{
Info<< " Applying source " << word(dict.lookup("source"))
<< endl;
autoPtr<topoSetSource> source = topoSetSource::New
//?
}
else
{
currentSet = topoSet::New
(
dict.lookup("source"),
setType,
mesh,
dict.subDict("sourceInfo")
setName,
IOobject::MUST_READ
);
Info<< "Read set " << setName << " with size "
<< currentSet().size() << endl;
}
// Backup current set.
autoPtr<topoSet> oldSet
(
topoSet::New
// Handle special actions (clear, invert) locally, rest through
// sources.
switch (action)
{
case topoSetSource::NEW:
case topoSetSource::ADD:
case topoSetSource::DELETE:
{
Info<< " Applying source " << word(dict.lookup("source"))
<< endl;
autoPtr<topoSetSource> source = topoSetSource::New
(
setType,
dict.lookup("source"),
mesh,
currentSet().name() + "_old2",
currentSet()
)
);
dict.subDict("sourceInfo")
);
currentSet().clear();
source().applyToSet(topoSetSource::NEW, currentSet());
source().applyToSet(action, currentSet());
// Synchronize for coupled patches.
if (!noSync) currentSet().sync(mesh);
currentSet().write();
}
break;
// Combine new value of currentSet with old one.
currentSet().subset(oldSet());
// Synchronize for coupled patches.
if (!noSync) currentSet().sync(mesh);
currentSet().write();
case topoSetSource::SUBSET:
{
Info<< " Applying source " << word(dict.lookup("source"))
<< endl;
autoPtr<topoSetSource> source = topoSetSource::New
(
dict.lookup("source"),
mesh,
dict.subDict("sourceInfo")
);
// Backup current set.
autoPtr<topoSet> oldSet
(
topoSet::New
(
setType,
mesh,
currentSet().name() + "_old2",
currentSet()
)
);
currentSet().clear();
source().applyToSet(topoSetSource::NEW, currentSet());
// Combine new value of currentSet with old one.
currentSet().subset(oldSet());
// Synchronize for coupled patches.
if (!noSync) currentSet().sync(mesh);
currentSet().write();
}
break;
case topoSetSource::CLEAR:
Info<< " Clearing set" << endl;
currentSet().clear();
currentSet().write();
break;
case topoSetSource::INVERT:
Info<< " Inverting set" << endl;
currentSet().invert(currentSet().maxSize(mesh));
currentSet().write();
break;
default:
WarningIn(args.executable())
<< "Unhandled action " << action << endl;
break;
}
break;
case topoSetSource::CLEAR:
Info<< " Clearing set" << endl;
currentSet().clear();
currentSet().write();
break;
case topoSetSource::INVERT:
Info<< " Inverting set" << endl;
currentSet().invert(currentSet().maxSize(mesh));
currentSet().write();
break;
default:
WarningIn(args.executable())
<< "Unhandled action " << action << endl;
break;
}
if (currentSet.valid())
{
Info<< " Set " << currentSet().name()
<< " now size " << currentSet().size()
<< endl;
if (currentSet.valid())
{
Info<< " Set " << currentSet().name()
<< " now size " << currentSet().size()
<< endl;
}
}
}

230
applications/utilities/surface/surfaceClean/collapseBase.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -66,8 +66,8 @@ static void writeRegionOBJ
triSurface regionSurf(surf.subsetMesh(include, pointMap, faceMap));
//Pout<< "Region " << regionI << " surface:" << nl;
//regionSurf.writeStats(Pout);
Pout<< "Region " << regionI << " surface:" << nl;
regionSurf.writeStats(Pout);
regionSurf.write(regionName);
@ -97,7 +97,7 @@ static void splitTri
DynamicList<labelledTri>& tris
)
{
label oldNTris = tris.size();
//label oldNTris = tris.size();
label fp = findIndex(f, e[0]);
label fp1 = f.fcIndex(fp);
@ -175,7 +175,7 @@ static void splitTri
}
else
{
FatalErrorIn("splitTri")
FatalErrorIn("splitTri(..)")
<< "Edge " << e << " not part of triangle " << f
<< " fp:" << fp
<< " fp1:" << fp1
@ -183,13 +183,13 @@ static void splitTri
<< abort(FatalError);
}
Pout<< "Split face " << f << " along edge " << e
<< " into triangles:" << endl;
for (label i = oldNTris; i < tris.size(); i++)
{
Pout<< " " << tris[i] << nl;
}
//Pout<< "Split face " << f << " along edge " << e
// << " into triangles:" << endl;
//
//for (label i = oldNTris; i < tris.size(); i++)
//{
// Pout<< " " << tris[i] << nl;
//}
}
@ -206,14 +206,14 @@ static bool insertSorted
{
if (findIndex(sortedVerts, vertI) != -1)
{
FatalErrorIn("insertSorted") << "Trying to insert vertex " << vertI
FatalErrorIn("insertSorted(..)") << "Trying to insert vertex " << vertI
<< " which is already in list of sorted vertices "
<< sortedVerts << abort(FatalError);
}
if (weight <= 0 || weight >= 1)
{
FatalErrorIn("insertSorted") << "Trying to insert vertex " << vertI
FatalErrorIn("insertSorted(..)") << "Trying to insert vertex " << vertI
<< " with illegal weight " << weight
<< " into list of sorted vertices "
<< sortedVerts << abort(FatalError);
@ -228,7 +228,7 @@ static bool insertSorted
if (mag(w - weight) < SMALL)
{
WarningIn("insertSorted")
WarningIn("insertSorted(..)")
<< "Trying to insert weight " << weight << " which is close to"
<< " existing weight " << w << " in " << sortedWeights
<< endl;
@ -263,64 +263,103 @@ static bool insertSorted
}
// Is triangle candidate for collapse? Small height or small quality
bool isSliver
(
const triSurface& surf,
const scalar minLen,
const scalar minQuality,
const label faceI,
const label edgeI
)
{
const pointField& localPoints = surf.localPoints();
// Check
// - opposite vertex projects onto base edge
// - normal distance is small
// - or triangle quality is small
label opposite0 =
triSurfaceTools::oppositeVertex
(
surf,
faceI,
edgeI
);
const edge& e = surf.edges()[edgeI];
const labelledTri& f = surf[faceI];
pointHit pHit =
e.line(localPoints).nearestDist
(
localPoints[opposite0]
);
if
(
pHit.hit()
&& (
pHit.distance() < minLen
|| f.tri(surf.points()).quality() < minQuality
)
)
{
// Remove faceI and split all other faces using this
// edge. This is done by 'replacing' the edgeI with the
// opposite0 vertex
//Pout<< "Splitting face " << faceI << " since distance "
// << pHit.distance()
// << " from vertex " << opposite0
// << " to edge " << edgeI
// << " points "
// << localPoints[e[0]]
// << localPoints[e[1]]
// << " is too small or triangle quality "
// << f.tri(surf.points()).quality()
// << " too small." << endl;
return true;
}
else
{
return false;
}
}
// Mark all faces that are going to be collapsed.
// faceToEdge: per face -1 or the base edge of the face.
static void markCollapsedFaces
(
const triSurface& surf,
const scalar minLen,
const scalar minQuality,
labelList& faceToEdge
)
{
faceToEdge.setSize(surf.size());
faceToEdge = -1;
const pointField& localPoints = surf.localPoints();
const labelListList& edgeFaces = surf.edgeFaces();
forAll(edgeFaces, edgeI)
{
const edge& e = surf.edges()[edgeI];
const labelList& eFaces = surf.edgeFaces()[edgeI];
forAll(eFaces, i)
{
label faceI = eFaces[i];
// Check distance of vertex to edge.
label opposite0 =
triSurfaceTools::oppositeVertex
(
surf,
faceI,
edgeI
);
bool isCandidate = isSliver(surf, minLen, minQuality, faceI, edgeI);
pointHit pHit =
e.line(localPoints).nearestDist
(
localPoints[opposite0]
);
if (pHit.hit() && pHit.distance() < minLen)
if (isCandidate)
{
// Remove faceI and split all other faces using this
// edge. This is done by 'replacing' the edgeI with the
// opposite0 vertex
Pout<< "Splitting face " << faceI << " since distance "
<< pHit.distance()
<< " from vertex " << opposite0
<< " to edge " << edgeI
<< " points "
<< localPoints[e[0]]
<< localPoints[e[1]]
<< " is too small" << endl;
// Mark face as being collapsed
if (faceToEdge[faceI] != -1)
{
FatalErrorIn("markCollapsedFaces")
FatalErrorIn("markCollapsedFaces(..)")
<< "Cannot collapse face " << faceI << " since "
<< " is marked to be collapsed both to edge "
<< faceToEdge[faceI] << " and " << edgeI
@ -347,7 +386,7 @@ static void markRegion
{
if (faceToEdge[faceI] == -1 || collapseRegion[faceI] != -1)
{
FatalErrorIn("markRegion")
FatalErrorIn("markRegion(..)")
<< "Problem : crossed into uncollapsed/regionized face"
<< abort(FatalError);
}
@ -383,7 +422,7 @@ static void markRegion
}
else if (collapseRegion[nbrFaceI] != regionI)
{
FatalErrorIn("markRegion")
FatalErrorIn("markRegion(..)")
<< "Edge:" << edgeI << " between face " << faceI
<< " with region " << regionI
<< " and face " << nbrFaceI
@ -411,8 +450,8 @@ static label markRegions
{
if (collapseRegion[faceI] == -1 && faceToEdge[faceI] != -1)
{
Pout<< "markRegions : Marking region:" << regionI
<< " starting from face " << faceI << endl;
//Pout<< "markRegions : Marking region:" << regionI
// << " starting from face " << faceI << endl;
// Collapsed face. Mark connected region with current region number
markRegion(surf, faceToEdge, regionI++, faceI, collapseRegion);
@ -728,7 +767,12 @@ static void getSplitVerts
}
label collapseBase(triSurface& surf, const scalar minLen)
label collapseBase
(
triSurface& surf,
const scalar minLen,
const scalar minQuality
)
{
label nTotalSplit = 0;
@ -743,7 +787,7 @@ label collapseBase(triSurface& surf, const scalar minLen)
labelList faceToEdge(surf.size(), -1);
// Calculate faceToEdge (face collapses)
markCollapsedFaces(surf, minLen, faceToEdge);
markCollapsedFaces(surf, minLen, minQuality, faceToEdge);
// Find regions of connected collapsed faces
@ -754,8 +798,8 @@ label collapseBase(triSurface& surf, const scalar minLen)
label nRegions = markRegions(surf, faceToEdge, collapseRegion);
Pout<< "Detected " << nRegions << " regions of faces to be collapsed"
<< nl << endl;
//Pout<< "Detected " << nRegions << " regions of faces to be collapsed"
// << nl << endl;
// Pick up all vertices on outside of region
labelListList outsideVerts
@ -772,10 +816,10 @@ label collapseBase(triSurface& surf, const scalar minLen)
{
spanPoints[regionI] = getSpanPoints(surf, outsideVerts[regionI]);
Pout<< "For region " << regionI << " found extrema at points "
<< surf.localPoints()[spanPoints[regionI][0]]
<< surf.localPoints()[spanPoints[regionI][1]]
<< endl;
//Pout<< "For region " << regionI << " found extrema at points "
// << surf.localPoints()[spanPoints[regionI][0]]
// << surf.localPoints()[spanPoints[regionI][1]]
// << endl;
// Project all non-span points onto the span edge.
projectNonSpanPoints
@ -787,21 +831,21 @@ label collapseBase(triSurface& surf, const scalar minLen)
orderedWeights[regionI]
);
Pout<< "For region:" << regionI
<< " span:" << spanPoints[regionI]
<< " orderedVerts:" << orderedVertices[regionI]
<< " orderedWeights:" << orderedWeights[regionI]
<< endl;
//Pout<< "For region:" << regionI
// << " span:" << spanPoints[regionI]
// << " orderedVerts:" << orderedVertices[regionI]
// << " orderedWeights:" << orderedWeights[regionI]
// << endl;
writeRegionOBJ
(
surf,
regionI,
collapseRegion,
outsideVerts[regionI]
);
//writeRegionOBJ
//(
// surf,
// regionI,
// collapseRegion,
// outsideVerts[regionI]
//);
Pout<< endl;
//Pout<< endl;
}
@ -864,20 +908,19 @@ label collapseBase(triSurface& surf, const scalar minLen)
// Split edge using splitVerts. All non-collapsed triangles
// using edge will get split.
{
const pointField& localPoints = surf.localPoints();
Pout<< "edge " << edgeI << ' ' << e
<< " points "
<< localPoints[e[0]] << ' ' << localPoints[e[1]]
<< " split into edges with extra points:"
<< endl;
forAll(splitVerts, i)
{
Pout<< " " << splitVerts[i] << " weight "
<< splitWeights[i] << nl;
}
}
//{
// const pointField& localPoints = surf.localPoints();
// Pout<< "edge " << edgeI << ' ' << e
// << " points "
// << localPoints[e[0]] << ' ' << localPoints[e[1]]
// << " split into edges with extra points:"
// << endl;
// forAll(splitVerts, i)
// {
// Pout<< " " << splitVerts[i] << " weight "
// << splitWeights[i] << nl;
// }
//}
const labelList& eFaces = surf.edgeFaces()[edgeI];
@ -914,7 +957,8 @@ label collapseBase(triSurface& surf, const scalar minLen)
}
}
Pout<< "collapseBase : splitting " << nSplit << " triangles"
Info<< "collapseBase : collapsing " << nSplit
<< " triangles by splitting their base edge."
<< endl;
nTotalSplit += nSplit;
@ -927,15 +971,15 @@ label collapseBase(triSurface& surf, const scalar minLen)
// Pack the triangles
newTris.shrink();
Pout<< "Resetting surface from " << surf.size() << " to "
<< newTris.size() << " triangles" << endl;
//Pout<< "Resetting surface from " << surf.size() << " to "
// << newTris.size() << " triangles" << endl;
surf = triSurface(newTris, surf.patches(), surf.localPoints());
{
fileName fName("bla" + name(iter) + ".obj");
Pout<< "Writing surf to " << fName << endl;
surf.write(fName);
}
//{
// fileName fName("bla" + name(iter) + ".obj");
// Pout<< "Writing surf to " << fName << endl;
// surf.write(fName);
//}
iter++;
}

11
applications/utilities/surface/surfaceClean/collapseBase.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -41,9 +41,14 @@ using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Keep collapsing all triangles whose height is < minLen.
//- Keep collapsing all triangles whose height is < minLen or quality < minQ.
// Returns number of triangles collapsed.
label collapseBase(triSurface& surf, const scalar minLen);
label collapseBase
(
triSurface& surf,
const scalar minLen,
const scalar minQuality
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

9
applications/utilities/surface/surfaceClean/collapseEdge.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -121,8 +121,8 @@ label collapseEdge(triSurface& surf, const scalar minLen)
pointMap[v1] = v;
newPoints[v] = 0.5*(localPoints[v1] + localPoints[v]);
Pout<< "Collapsing triange " << faceI << " to edge mid "
<< newPoints[v] << endl;
//Pout<< "Collapsing triange " << faceI
// << " to edge mid " << newPoints[v] << endl;
nCollapsed++;
okToCollapse[faceI] = false;
@ -136,7 +136,8 @@ label collapseEdge(triSurface& surf, const scalar minLen)
}
}
Pout<< "collapseEdge : collapsing " << nCollapsed << " triangles"
Info<< "collapseEdge : collapsing " << nCollapsed
<< " triangles to a single edge."
<< endl;
nTotalCollapsed += nCollapsed;

10
applications/utilities/surface/surfaceClean/surfaceClean.C
View File

@ -50,6 +50,7 @@ int main(int argc, char *argv[])
argList::noParallel();
argList::validArgs.append("surfaceFile");
argList::validArgs.append("min length");
argList::validArgs.append("min quality");
argList::validArgs.append("output surfaceFile");
argList::addBoolOption
(
@ -60,10 +61,13 @@ int main(int argc, char *argv[])
const fileName inFileName = args[1];
const scalar minLen = args.argRead<scalar>(2);
const fileName outFileName = args[3];
const scalar minQuality = args.argRead<scalar>(3);
const fileName outFileName = args[4];
Info<< "Reading surface " << inFileName << nl
<< "Collapsing all triangles with edges or heights < " << minLen << nl
<< "Collapsing all triangles with" << nl
<< " edges or heights < " << minLen << nl
<< " quality < " << minQuality << nl
<< "Writing result to " << outFileName << nl << endl;
@ -90,7 +94,7 @@ int main(int argc, char *argv[])
}
while (true)
{
label nSplitEdge = collapseBase(surf, minLen);
label nSplitEdge = collapseBase(surf, minLen, minQuality);
if (nSplitEdge == 0)
{

4
applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C
View File

@ -51,8 +51,8 @@ void Foam::buildCGALPolyhedron::operator()
HalfedgeDS& hds
)
{
typedef typename HalfedgeDS::Traits Traits;
typedef typename Traits::Point_3 Point;
typedef HalfedgeDS::Traits Traits;
typedef Traits::Point_3 Point;
// Postcondition: `hds' is a valid polyhedral surface.
CGAL::Polyhedron_incremental_builder_3<HalfedgeDS> B(hds, false);

40
differencesWrtDev.txt
View File

@ -1,7 +1,7 @@
src/
mesh/conformalVoronoiMesh
- all the meshing.
- all the meshing. See separate section below.
polyMeshGeometry:
- disabled tetquality check on face-center decomp tet.
@ -93,5 +93,43 @@ surfaceSplitByTopology:
conformalVoronoiMesh/
--------------------
- Make sure the surface does not have any sliver triangles. These are
hard to get the surface normal correct for so might cause bleeding.
- Use surfaceCheck to find out triangle quality and size of smallest edges
- Use surfaceClean .. 1e-5 .. to get rid of any edges < 1e-5m.
- If you get bleeding you might see in face filtering:
...
cells with with zero or one non-boundary face : 1
...
Initial check before face collapse, found 48 bad quality faces
but this was real - the cell that got created inside the cone and sphere by the
bad triangle was actually attached to a lot of faces. This screwed up the
subsequent filtering as it stopped too early.
I ran:
surfaceClean coneAndSphere.obj 1e-5 coneAndSphere_clean.obj
and re-ran with that surface and got
...
cells with with zero or one non-boundary face : 0
...
Initial check before face collapse, found 0 bad quality faces
and the bad cells inside are gone.
That group of cells would be picked up at the end of the meshing as the
cvMesh_remainingProtrusions
cellSet which can be deleted.

166
doc/cvMesh.org Normal file
View File

@ -0,0 +1,166 @@
# -*- mode: org; -*-
#
#+TITLE: cvMesh usage
#+AUTHOR: OpenCFD Ltd.
#+DATE: June 2011
#+LINK: http://www.OpenFOAM.com
#+OPTIONS: author:nil ^:{}
#+STARTUP: hidestars
#+STARTUP: odd
* Method
+ initial point seeding. These are the background points to
- start the triangulation from (a subset)
- do cellsize and alignment queries on
+ mesh movement based on cell size and alignment
+ surface conformation by introducing duplicate points
+ do final simplification by face filtering
* Typical work flow
+ clean surfaces (surfaceClean, surfaceOrient)
+ set in system/controlDict
- startTime to latestTime
- endTime to number of meshing iterations
- writeInterval to large value
+ select defaultCellSize [m]
+ select refinement patterns based on e.g. surface
(cellSizeControlGeometry)
+ select autoDensity initial point seeding
+ force output of unfiltered mesh (to save time) by specifying unobtainable
mesh quality:
continueFilteringOnBadInitialPolyMesh false;
maxNonOrtho 1;
+ check visually
- features resolved?
- cell size ok? Load in targetCellSize volScalarField.
- if any surface protusions that are unacceptable
- lower targetCellSize
- lower maxSurfaceProtrusionCoeff (usually 0.1 of targetCellSize)
+ restart. Repeat until surface good enough.
+ produce final mesh
- set initialPointsMethod to pointFile so it reuses the Delaunay vertices
- re-enable mesh filtering
continueFilteringOnBadInitialPolyMesh true;
maxNonOrtho 65; //or whatever;
* Debugging
+ run surfaceFeatureExtract with -writeObj to dump all features
+ set initialPointsMethod to pointFile
+ set writeInterval to 1
+ run e.g. topoSet to create slices
* Annotated output
#+BEGIN_EXAMPLE
Create time
Reading geometryToConformTo
Reading cellSizeControlGeometry
Selecting initialPointsMethod autoDensity
Selecting relaxationModel adaptiveLinear
Selecting faceAreaWeightModel piecewiseLinearRamp
Conforming to feature points
Inserted 32 feature vertices
#+END_EXAMPLE
In this particular case there are 8 straight corners so these can be recreated
with a triplet each.
#+BEGIN_EXAMPLE
Inserting initial points
autoDensity
6922 points placed
9139 locations queried
0.7574132837 success rate
#+END_EXAMPLE
The initial points get created by the =autoDensity= method. It will accept only
part (75% in this cae) of the initial points to create randomly distributed locations.
#+BEGIN_EXAMPLE
6922 points to insert...
6922 points inserted
Build coarse surface conformation
Initial conformation
Number of vertices 6962
Number of surface hits 2245
Number of edge hits 238
#+END_EXAMPLE
The initial points get inserted into the triangulation and it detects the number
of surface pierces (surface hits) and the number of XXX (edge hits).
#+BEGIN_EXAMPLE
Stopping iterations when:
total number of hits drops below 0.001 of initial hits (2)
or
maximum number of iterations (15) is reached
#+END_EXAMPLE
It will introduce additional duplicate vertices on either side of the surface
pierce to conform to the geometry. This needs to repeated until the number of
pierces is low enough.
#+BEGIN_EXAMPLE
Conformation iteration 0
Number of vertices 12166
Number of surface hits 234
Number of edge hits 38
..
Conformation iteration 4
Number of vertices 12908
Number of surface hits 0
Number of edge hits 0
--- [ cpuTime 5.46 s, delta 0.8 s, "Conformation iteration 4" ] ---
--- [ 5.46 s, mem size 114228 kB, mem peak 114228 kB, mem rss 28712 kB ] ---
Total hits (0) less than limit (2), stopping iterations
Storing surface conformation
Stored 5946 vertices
Store size and alignment
#+END_EXAMPLE
The vertices needed to create the conformation are stored together with the size
and alignment.
#+BEGIN_EXAMPLE
Time = 1
Relaxation = 1
Looking up target cell alignment and size
Determining vertex displacements
Reinserting stored feature points
Reinserted 32 vertices
Inserting displaced tessellation
4328 points to insert...
4328 points inserted
Reinserting stored surface conformation
Reinserted 5946 vertices
Total displacement = (0.03629314022 0.1324815449 -0.2774017193)
Total distance = 159.3482907
ExecutionTime = 8.16 s ClockTime = 8 s
Time = 2
Relaxation = 0.9875
#+END_EXAMPLE

1
src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H
View File

@ -37,7 +37,6 @@ SourceFiles
#include "searchableSurfaces.H"
#include "searchableSurfacesQueries.H"
#include "triSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

24
src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
View File

@ -52,7 +52,6 @@ SourceFiles
#include "plane.H"
#include "SortableList.H"
#include "meshTools.H"
#include "triSurfaceTools.H"
#include "indexedOctree.H"
#include "treeDataPoint.H"
#include "unitConversion.H"
@ -683,17 +682,18 @@ private:
bool filterFaces
);
// //- Tet mesh calculation
// void calcTetMesh
// (
// pointField& points,
// faceList& faces,
// labelList& owner,
// labelList& neighbour,
// wordList& patchNames,
// labelList& patchSizes,
// labelList& patchStarts
// );
//- Tet mesh calculation
void calcTetMesh
(
pointField& points,
faceList& faces,
labelList& owner,
labelList& neighbour,
wordList& patchTypes,
wordList& patchNames,
labelList& patchSizes,
labelList& patchStarts
);
//- Determines if the dual face constructed by the Delaunay
// edge is a boundary face

2
src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
View File

@ -38,7 +38,7 @@ SourceFiles
#include "searchableSurfaces.H"
#include "searchableSurfacesQueries.H"
#include "extendedFeatureEdgeMesh.H"
#include "triSurface.H"
#include "boolList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@ -331,6 +331,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
),
triSurface(s),
tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
minQuality_(-1),
maxTreeDepth_(10),
surfaceClosed_(-1)
{
@ -377,6 +378,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
)
),
tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
minQuality_(-1),
maxTreeDepth_(10),
surfaceClosed_(-1)
{
@ -426,6 +428,7 @@ Foam::triSurfaceMesh::triSurfaceMesh
)
),
tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
minQuality_(-1),
maxTreeDepth_(10),
surfaceClosed_(-1)
{
@ -449,6 +452,13 @@ Foam::triSurfaceMesh::triSurfaceMesh
<< tolerance_ << endl;
}
// Have optional minimum quality for normal calculation
if (dict.readIfPresent("minQuality", minQuality_) && minQuality_ > 0)
{
Info<< searchableSurface::name()
<< " : ignoring triangles with quality < "
<< minQuality_ << " for normals calculation." << endl;
}
// Have optional non-standard tree-depth to limit storage.
if (dict.readIfPresent("maxTreeDepth", maxTreeDepth_) && maxTreeDepth_ > 0)
@ -820,22 +830,70 @@ void Foam::triSurfaceMesh::getNormal
{
normal.setSize(info.size());
forAll(info, i)
if (minQuality_ >= 0)
{
if (info[i].hit())
{
label faceI = info[i].index();
//- Cached:
//normal[i] = faceNormals()[faceI];
// Make sure we don't use triangles with low quality since
// normal is not reliable.
//- Uncached
normal[i] = triSurface::operator[](faceI).normal(points());
normal[i] /= mag(normal[i]) + VSMALL;
}
else
const triSurface& s = static_cast<const triSurface&>(*this);
const labelListList& faceFaces = s.faceFaces();
forAll(info, i)
{
// Set to what?
normal[i] = vector::zero;
if (info[i].hit())
{
label faceI = info[i].index();
scalar qual = s[faceI].tri(points()).quality();
if (qual < minQuality_)
{
// Search neighbouring triangles
const labelList& fFaces = faceFaces[faceI];
forAll(fFaces, j)
{
label nbrI = fFaces[j];
scalar nbrQual = s[nbrI].tri(points()).quality();
if (nbrQual > qual)
{
qual = nbrQual;
normal[i] = s[nbrI].normal(points());
}
}
}
else
{
normal[i] = s[faceI].normal(points());
}
normal[i] /= mag(normal[i]);
}
else
{
// Set to what?
normal[i] = vector::zero;
}
}
}
else
{
forAll(info, i)
{
if (info[i].hit())
{
label faceI = info[i].index();
//- Cached:
//normal[i] = faceNormals()[faceI];
//- Uncached
normal[i] = triSurface::operator[](faceI).normal(points());
normal[i] /= mag(normal[i]) + VSMALL;
}
else
{
// Set to what?
normal[i] = vector::zero;
}
}
}
}

7
src/meshTools/searchableSurface/triSurfaceMesh.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -31,6 +31,7 @@ Description
- scale : scaling factor.
- tolerance : relative tolerance for doing intersections
(see triangle::intersection)
- minQuality: discard triangles with low quality when getting normal
SourceFiles
triSurfaceMesh.C
@ -70,6 +71,10 @@ private:
//- Optional tolerance to use in searches
scalar tolerance_;
//- Optional min triangle quality. Triangles below this get
// ignored for normal calculation
scalar minQuality_;
//- Optional max tree depth of octree
label maxTreeDepth_;