songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
Henry 2012-10-04 10:08:55 +01:00
commit 80b2446564
16 changed files with 533 additions and 277 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

144
applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
View File

@ -7,6 +7,7 @@
#include "pointSet.H"
#include "IOmanip.H"
#include "emptyPolyPatch.H"
#include "processorPolyPatch.H"
Foam::label Foam::checkTopology
(
@ -289,10 +290,18 @@ Foam::label Foam::checkTopology
}
}
if (!Pstream::parRun())
{
Pout<< "\nChecking patch topology for multiply connected surfaces ..."
<< endl;
if (!Pstream::parRun())
{
Info<< "\nChecking patch topology for multiply connected"
<< " surfaces..." << endl;
}
else
{
Info<< "\nChecking basic patch addressing..." << endl;
}
const polyBoundaryMesh& patches = mesh.boundaryMesh();
@ -301,90 +310,101 @@ Foam::label Foam::checkTopology
(
mesh,
"nonManifoldPoints",
mesh.nPoints()/100
mesh.nPoints()/1000
);
Pout.setf(ios_base::left);
Pout<< " "
Info<< " "
<< setw(20) << "Patch"
<< setw(9) << "Faces"
<< setw(9) << "Points"
<< setw(34) << "Surface topology";
<< setw(9) << "Points";
if (!Pstream::parRun())
{
Info<< setw(34) << "Surface topology";
}
if (allGeometry)
{
Pout<< " Bounding box";
Info<< " Bounding box";
}
Pout<< endl;
Info<< endl;
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
Pout<< " "
if (!isA<processorPolyPatch>(pp))
{
Info<< " "
<< setw(20) << pp.name()
<< setw(9) << pp.size()
<< setw(9) << pp.nPoints();
<< setw(9) << returnReduce(pp.size(), sumOp<label>())
<< setw(9) << returnReduce(pp.nPoints(), sumOp<label>());
primitivePatch::surfaceTopo pTyp = pp.surfaceType();
if (pp.empty())
{
Pout<< setw(34) << "ok (empty)";
}
else if (pTyp == primitivePatch::MANIFOLD)
{
if (pp.checkPointManifold(true, &points))
if (!Pstream::parRun())
{
Pout<< setw(34) << "multiply connected (shared point)";
}
else
{
Pout<< setw(34) << "ok (closed singly connected)";
}
primitivePatch::surfaceTopo pTyp = pp.surfaceType();
// Add points on non-manifold edges to make set complete
pp.checkTopology(false, &points);
}
else
{
pp.checkTopology(false, &points);
if (pTyp == primitivePatch::OPEN)
{
Pout<< setw(34) << "ok (non-closed singly connected)";
}
else
{
Pout<< setw(34) << "multiply connected (shared edge)";
}
}
if (allGeometry)
{
const pointField& pts = pp.points();
const labelList& mp = pp.meshPoints();
if (returnReduce(mp.size(), sumOp<label>()) > 0)
{
boundBox bb(point::max, point::min);
forAll (mp, i)
if (pp.empty())
{
bb.min() = min(bb.min(), pts[mp[i]]);
bb.max() = max(bb.max(), pts[mp[i]]);
Info<< setw(34) << "ok (empty)";
}
else if (pTyp == primitivePatch::MANIFOLD)
{
if (pp.checkPointManifold(true, &points))
{
Info<< setw(34)
<< "multiply connected (shared point)";
}
else
{
Info<< setw(34) << "ok (closed singly connected)";
}
// Add points on non-manifold edges to make set complete
pp.checkTopology(false, &points);
}
else
{
pp.checkTopology(false, &points);
if (pTyp == primitivePatch::OPEN)
{
Info<< setw(34)
<< "ok (non-closed singly connected)";
}
else
{
Info<< setw(34)
<< "multiply connected (shared edge)";
}
}
reduce(bb.min(), minOp<vector>());
reduce(bb.max(), maxOp<vector>());
Pout<< ' ' << bb;
}
if (allGeometry)
{
const pointField& pts = pp.points();
const labelList& mp = pp.meshPoints();
if (returnReduce(mp.size(), sumOp<label>()) > 0)
{
boundBox bb(point::max, point::min);
forAll (mp, i)
{
bb.min() = min(bb.min(), pts[mp[i]]);
bb.max() = max(bb.max(), pts[mp[i]]);
}
reduce(bb.min(), minOp<vector>());
reduce(bb.max(), maxOp<vector>());
Info<< ' ' << bb;
}
}
Info<< endl;
}
Pout<< endl;
}
if (points.size())
{
Pout<< " <<Writing " << points.size()
Info<< " <<Writing " << returnReduce(points.size(), sumOp<label>())
<< " conflicting points to set "
<< points.name() << endl;
@ -392,7 +412,7 @@ Foam::label Foam::checkTopology
points.write();
}
//Pout.setf(ios_base::right);
//Info.setf(ios_base::right);
}
// Force creation of all addressing if requested.

11
applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -115,6 +115,7 @@ int main(int argc, char *argv[])
runTimeMaster
)
);
const word oldInstance = masterMesh.pointsInstance();
Info<< "Reading mesh to add for time = " << runTimeToAdd.timeName() << nl;
@ -139,7 +140,13 @@ int main(int argc, char *argv[])
masterMesh.addMesh(meshToAdd);
masterMesh.merge();
masterMesh.polyMesh::write();
if (overwrite)
{
masterMesh.setInstance(oldInstance);
}
masterMesh.write();
Info<< "\nEnd\n" << endl;

22
src/OSspecific/POSIX/printStack.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -167,6 +167,26 @@ void getSymbolForRaw
}
void error::safePrintStack(std::ostream& os)
{
// Get raw stack symbols
void *array[100];
size_t size = backtrace(array, 100);
char **strings = backtrace_symbols(array, size);
// See if they contain function between () e.g. "(__libc_start_main+0xd0)"
// and see if cplus_demangle can make sense of part before +
for (size_t i = 0; i < size; i++)
{
string msg(strings[i]);
fileName programFile;
word address;
os << '#' << label(i) << '\t' << msg << std::endl;
}
}
void error::printStack(Ostream& os)
{
// Reads the starting addresses for the dynamically linked libraries

5
src/OpenFOAM/db/IOstreams/Sstreams/ISstream.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,8 +29,7 @@ Description
SourceFiles
ISstreamI.H
ISread.C
ISreadToken.C
ISstream.C
\*---------------------------------------------------------------------------*/

6
src/OpenFOAM/db/error/error.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -162,6 +162,10 @@ public:
operator dictionary() const;
//- Helper function to print a stack (if OpenFOAM IO not yet
// initialised)
static void safePrintStack(std::ostream&);
//- Helper function to print a stack
static void printStack(Ostream&);

13
src/OpenFOAM/fields/pointPatchFields/basic/fixedValue/fixedValuePointPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -121,6 +121,17 @@ public:
}
// Member functions
// Access
//- Return true if this patch field fixes a value
virtual bool fixesValue() const
{
return true;
}
// Member operators
// Disable assignment operators

8
src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -278,6 +278,12 @@ public:
return internalField_;
}
//- Return true if this patch field fixes a value
virtual bool fixesValue() const
{
return false;
}
//- Return true if this patch field is coupled
virtual bool coupled() const
{

43
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
View File

@ -710,7 +710,8 @@ void Foam::autoLayerDriver::setNumLayers
const indirectPrimitivePatch& pp,
pointField& patchDisp,
labelList& patchNLayers,
List<extrudeMode>& extrudeStatus
List<extrudeMode>& extrudeStatus,
label& nAddedCells
) const
{
const fvMesh& mesh = meshRefiner_.mesh();
@ -797,6 +798,24 @@ void Foam::autoLayerDriver::setNumLayers
}
}
// Calculate number of cells to create
nAddedCells = 0;
forAll(pp.localFaces(), faceI)
{
const face& f = pp.localFaces()[faceI];
// Get max of extrusion per point
label nCells = 0;
forAll(f, fp)
{
nCells = max(nCells, patchNLayers[f[fp]]);
}
nAddedCells += nCells;
}
reduce(nAddedCells, sumOp<label>());
//reduce(nConflicts, sumOp<label>());
//
//Info<< "Set displacement to zero for " << nConflicts
@ -2491,9 +2510,13 @@ void Foam::autoLayerDriver::addLayers
// extrudeStatus = EXTRUDE.
labelList patchNLayers(pp().nPoints(), 0);
// Ideal number of cells added
label nIdealAddedCells = 0;
// Whether to add edge for all pp.localPoints.
List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
{
// Get number of layer per point from number of layers per patch
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -2506,7 +2529,8 @@ void Foam::autoLayerDriver::addLayers
patchDisp,
patchNLayers,
extrudeStatus
extrudeStatus,
nIdealAddedCells
);
// Precalculate mesh edge labels for patch edges
@ -2776,6 +2800,7 @@ void Foam::autoLayerDriver::addLayers
layerParams.nSmoothNormals(),
layerParams.nSmoothSurfaceNormals(),
layerParams.minMedianAxisAngleCos(),
layerParams.featureAngle(),
dispVec,
medialRatio,
@ -3101,10 +3126,24 @@ void Foam::autoLayerDriver::addLayers
label nExtruded = countExtrusion(pp, extrudeStatus);
label nTotFaces = returnReduce(pp().size(), sumOp<label>());
label nAddedCells = 0;
{
forAll(flaggedCells, cellI)
{
if (flaggedCells[cellI])
{
nAddedCells++;
}
}
reduce(nAddedCells, sumOp<label>());
}
Info<< "Extruding " << nExtruded
<< " out of " << nTotFaces
<< " faces (" << 100.0*nExtruded/nTotFaces << "%)."
<< " Removed extrusion at " << nTotChanged << " faces."
<< endl
<< "Added " << nAddedCells << " out of " << nIdealAddedCells
<< " cells (" << 100.0*nAddedCells/nIdealAddedCells << "%)."
<< endl;
if (nTotChanged == 0)

4
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
View File

@ -192,7 +192,8 @@ class autoLayerDriver
const indirectPrimitivePatch& pp,
pointField& patchDisp,
labelList& patchNLayers,
List<extrudeMode>& extrudeStatus
List<extrudeMode>& extrudeStatus,
label& nIdealAddedCells
) const;
//- Helper function to make a pointVectorField with correct
@ -469,6 +470,7 @@ class autoLayerDriver
const label nSmoothNormals,
const label nSmoothSurfaceNormals,
const scalar minMedianAxisAngleCos,
const scalar featureAngle,
pointVectorField& dispVec,
pointScalarField& medialRatio,

143
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
View File

@ -33,9 +33,10 @@ Description
#include "motionSmoother.H"
#include "pointData.H"
#include "PointEdgeWave.H"
#include "OFstream.H"
#include "OBJstream.H"
#include "meshTools.H"
#include "PatchTools.H"
#include "unitConversion.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -775,6 +776,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
const label nSmoothNormals,
const label nSmoothSurfaceNormals,
const scalar minMedianAxisAngleCos,
const scalar featureAngle,
pointVectorField& dispVec,
pointScalarField& medialRatio,
@ -904,7 +906,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
const edge& e = edges[edgeI];
// Approximate medial axis location on edge.
//const point medialAxisPt = e.centre(points);
vector eVec = e.vec(mesh.points());
vector eVec = e.vec(points);
scalar eMag = mag(eVec);
if (eMag > VSMALL)
{
@ -961,42 +963,111 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
labelHashSet adaptPatches(meshMover.adaptPatchIDs());
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
//const pointPatchVectorField& pvf =
// meshMover.displacement().boundaryField()[patchI];
const pointPatchVectorField& pvf =
meshMover.displacement().boundaryField()[patchI];
if
(
!pp.coupled()
&& !isA<emptyPolyPatch>(pp)
//&& pvf.constraintType() != word::null //exclude fixedValue
&& !adaptPatches.found(patchI)
)
{
Info<< "Inserting points on patch " << pp.name() << endl;
const labelList& meshPoints = pp.meshPoints();
forAll(meshPoints, i)
// Check the type of the patchField. The types are
// - fixedValue (0 or more layers) but the >0 layers have
// already been handled in the adaptPatches loop
// - constraint (but not coupled) types, e.g. symmetryPlane,
// slip.
if (pvf.fixesValue())
{
label pointI = meshPoints[i];
// Disable all movement on fixedValue patchFields
Info<< "Inserting all points on patch " << pp.name()
<< endl;
if (!pointMedialDist[pointI].valid(dummyTrackData))
forAll(meshPoints, i)
{
maxPoints.append(pointI);
maxInfo.append
(
pointData
label pointI = meshPoints[i];
if (!pointMedialDist[pointI].valid(dummyTrackData))
{
maxPoints.append(pointI);
maxInfo.append
(
points[pointI],
0.0,
pointI, // passive data
vector::zero // passive data
)
);
pointMedialDist[pointI] = maxInfo.last();
pointData
(
points[pointI],
0.0,
pointI, // passive data
vector::zero // passive data
)
);
pointMedialDist[pointI] = maxInfo.last();
}
}
}
else
{
// Based on geometry: analyse angle w.r.t. nearest moving
// point. In the pointWallDist we transported the
// normal as the passive vector. Note that this points
// out of the originating wall so inside of the domain
// on this patch.
Info<< "Inserting points on patch " << pp.name()
<< " if angle to nearest layer patch > "
<< featureAngle << " degrees." << endl;
scalar featureAngleCos = Foam::cos(degToRad(featureAngle));
pointField pointNormals
(
PatchTools::pointNormals
(
mesh,
pp,
identity(pp.size())+pp.start()
)
);
forAll(meshPoints, i)
{
label pointI = meshPoints[i];
if (!pointMedialDist[pointI].valid(dummyTrackData))
{
// Check if angle not too large.
scalar cosAngle =
(
-pointWallDist[pointI].v()
& pointNormals[i]
);
if (cosAngle > featureAngleCos)
{
// Extrusion direction practically perpendicular
// to the patch. Disable movement at the patch.
maxPoints.append(pointI);
maxInfo.append
(
pointData
(
points[pointI],
0.0,
pointI, // passive data
vector::zero // passive data
)
);
pointMedialDist[pointI] = maxInfo.last();
}
else
{
// Extrusion direction makes angle with patch
// so allow sliding - don't insert zero points
}
}
}
}
}
@ -1141,13 +1212,12 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
// reduce thickness where thickness/medial axis distance large
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
autoPtr<OFstream> str;
label vertI = 0;
autoPtr<OBJstream> str;
if (debug)
{
str.reset
(
new OFstream
new OBJstream
(
mesh.time().path()
/ "thicknessRatioExcludePoints_"
@ -1159,13 +1229,12 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
<< str().name() << endl;
}
autoPtr<OFstream> medialVecStr;
label medialVertI = 0;
autoPtr<OBJstream> medialVecStr;
if (debug)
{
medialVecStr.reset
(
new OFstream
new OBJstream
(
mesh.time().path()
/ "thicknessRatioExcludeMedialVec_"
@ -1227,21 +1296,19 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
if (str.valid())
{
const point& pt = mesh.points()[pointI];
meshTools::writeOBJ(str(), pt);
vertI++;
meshTools::writeOBJ(str(), pt+patchDisp[patchPointI]);
vertI++;
str()<< "l " << vertI-1 << ' ' << vertI << nl;
str().write(linePointRef(pt, pt+patchDisp[patchPointI]));
}
if (medialVecStr.valid())
{
const point& pt = mesh.points()[pointI];
meshTools::writeOBJ(medialVecStr(), pt);
medialVertI++;
meshTools::writeOBJ(medialVecStr(), medialVec[pointI]);
medialVertI++;
medialVecStr()<< "l " << medialVertI-1
<< ' ' << medialVertI << nl;
medialVecStr().write
(
linePointRef
(
pt,
medialVec[pointI]
)
);
}
}
}

18
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
View File

@ -178,6 +178,24 @@ class autoSnapDriver
vector& edgeOffset // offset from pt to point on edge
) const;
//- For any reverse (so from feature back to mesh) attraction:
// add attraction if diagonal points on face attracted
void stringFeatureEdges
(
const label iter,
const scalar featureCos,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
const vectorField& rawPatchAttraction,
const List<pointConstraint>& rawPatchConstraints,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const;
//- Return hit if on multiple points
pointIndexHit findMultiPatchPoint
(

378
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
View File

@ -848,17 +848,17 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
vector d = r.refPoint()-pt;
d -= (d&n)*n;
// Correct for attraction to non-dominant face
correctAttraction
(
surfacePoints,
surfaceCount,
r.refPoint(),
n, // normalised normal
pt,
d // perpendicular offset vector
);
//// Correct for attraction to non-dominant face
//correctAttraction
//(
// surfacePoints,
// surfaceCount,
// r.refPoint(),
// n, // normalised normal
// pt,
//
// d // perpendicular offset vector
//);
// Trim to snap distance
if (magSqr(d) > sqr(snapDist[pointI]))
@ -893,6 +893,15 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
patchConstraint.applyConstraint(surfaceNormals[0]);
patchConstraint.applyConstraint(surfaceNormals[1]);
patchConstraint.applyConstraint(surfaceNormals[2]);
//Pout<< "# Feature point " << pt << nl;
//meshTools::writeOBJ(Pout, pt);
//meshTools::writeOBJ(Pout, surfacePoints[0]);
//meshTools::writeOBJ(Pout, surfacePoints[1]);
//meshTools::writeOBJ(Pout, surfacePoints[2]);
//Pout<< "l 1 2" << nl
// << "l 1 3" << nl
// << "l 1 4" << nl;
}
}
@ -1001,6 +1010,195 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
}
void Foam::autoSnapDriver::stringFeatureEdges
(
const label iter,
const scalar featureCos,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
const vectorField& rawPatchAttraction,
const List<pointConstraint>& rawPatchConstraints,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const
{
// Snap edges to feature edges
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Walk existing edges and snap remaining ones (that are marked as
// feature edges in rawPatchConstraints)
// What this does is fill in any faces where not all points
// on the face are being attracted:
/*
+
/ \
/ \
---+ +---
\ /
\ /
+
*/
// so the top and bottom will never get attracted since the nearest
// back from the feature edge will always be one of the left or right
// points since the face is diamond like. So here we walk the feature edges
// and add any non-attracted points.
while (true)
{
label nChanged = 0;
const labelListList& pointEdges = pp.pointEdges();
forAll(pointEdges, pointI)
{
if (patchConstraints[pointI].first() == 2)
{
const point& pt = pp.localPoints()[pointI];
const labelList& pEdges = pointEdges[pointI];
const vector& featVec = patchConstraints[pointI].second();
// Detect whether there are edges in both directions.
// (direction along the feature edge that is)
bool hasPos = false;
bool hasNeg = false;
forAll(pEdges, pEdgeI)
{
const edge& e = pp.edges()[pEdges[pEdgeI]];
label nbrPointI = e.otherVertex(pointI);
if (patchConstraints[nbrPointI].first() > 1)
{
const point& nbrPt = pp.localPoints()[nbrPointI];
const point featPt =
nbrPt + patchAttraction[nbrPointI];
const scalar cosAngle = (featVec & (featPt-pt));
if (cosAngle > 0)
{
hasPos = true;
}
else
{
hasNeg = true;
}
}
}
if (!hasPos || !hasNeg)
{
//Pout<< "**Detected feature string end at "
// << pp.localPoints()[pointI] << endl;
// No string. Assign best choice on either side
label bestPosPointI = -1;
scalar minPosDistSqr = GREAT;
label bestNegPointI = -1;
scalar minNegDistSqr = GREAT;
forAll(pEdges, pEdgeI)
{
const edge& e = pp.edges()[pEdges[pEdgeI]];
label nbrPointI = e.otherVertex(pointI);
if
(
patchConstraints[nbrPointI].first() <= 1
&& rawPatchConstraints[nbrPointI].first() > 1
)
{
const vector& nbrFeatVec =
rawPatchConstraints[pointI].second();
if (mag(featVec&nbrFeatVec) > featureCos)
{
// nbrPointI attracted to sameish feature
// Note: also check on position.
scalar d2 = magSqr
(
rawPatchAttraction[nbrPointI]
);
const point featPt =
pp.localPoints()[nbrPointI]
+ rawPatchAttraction[nbrPointI];
const scalar cosAngle =
(featVec & (featPt-pt));
if (cosAngle > 0)
{
if (!hasPos && d2 < minPosDistSqr)
{
minPosDistSqr = d2;
bestPosPointI = nbrPointI;
}
}
else
{
if (!hasNeg && d2 < minNegDistSqr)
{
minNegDistSqr = d2;
bestNegPointI = nbrPointI;
}
}
}
}
}
if (bestPosPointI != -1)
{
// Use reconstructed-feature attraction. Use only
// part of it since not sure...
//const point& bestPt =
// pp.localPoints()[bestPosPointI];
//Pout<< "**Overriding point " << bestPt
// << " on reconstructed feature edge at "
// << rawPatchAttraction[bestPosPointI]+bestPt
// << " to attracted-to-feature-edge." << endl;
patchAttraction[bestPosPointI] =
0.5*rawPatchAttraction[bestPosPointI];
patchConstraints[bestPosPointI] =
rawPatchConstraints[bestPosPointI];
nChanged++;
}
if (bestNegPointI != -1)
{
// Use reconstructed-feature attraction. Use only
// part of it since not sure...
//const point& bestPt =
// pp.localPoints()[bestNegPointI];
//Pout<< "**Overriding point " << bestPt
// << " on reconstructed feature edge at "
// << rawPatchAttraction[bestNegPointI]+bestPt
// << " to attracted-to-feature-edge." << endl;
patchAttraction[bestNegPointI] =
0.5*rawPatchAttraction[bestNegPointI];
patchConstraints[bestNegPointI] =
rawPatchConstraints[bestNegPointI];
nChanged++;
}
}
}
}
reduce(nChanged, sumOp<label>());
Info<< "Stringing feature edges : changed " << nChanged << " points"
<< endl;
if (nChanged == 0)
{
break;
}
}
}
Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
(
const indirectPrimitivePatch& pp,
@ -2113,156 +2311,20 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
// Walk existing edges and snap remaining ones (that are marked as
// feature edges in allPatchConstraints)
while (true)
{
label nChanged = 0;
stringFeatureEdges
(
iter,
featureCos,
const labelListList& pointEdges = pp.pointEdges();
forAll(pointEdges, pointI)
{
if (patchConstraints[pointI].first() == 2)
{
const point& pt = pp.localPoints()[pointI];
const labelList& pEdges = pointEdges[pointI];
const vector& featVec = patchConstraints[pointI].second();
pp,
snapDist,
// Detect whether there are edges in both directions.
// (direction along the feature edge that is)
bool hasPos = false;
bool hasNeg = false;
forAll(pEdges, pEdgeI)
{
const edge& e = pp.edges()[pEdges[pEdgeI]];
label nbrPointI = e.otherVertex(pointI);
if (patchConstraints[nbrPointI].first() > 1)
{
const point& nbrPt = pp.localPoints()[nbrPointI];
const point featPt =
nbrPt + patchAttraction[nbrPointI];
const scalar cosAngle = (featVec & (featPt-pt));
if (cosAngle > 0)
{
hasPos = true;
}
else
{
hasNeg = true;
}
}
}
if (!hasPos || !hasNeg)
{
//Pout<< "**Detected feature string end at "
// << pp.localPoints()[pointI] << endl;
// No string. Assign best choice on either side
label bestPosPointI = -1;
scalar minPosDistSqr = GREAT;
label bestNegPointI = -1;
scalar minNegDistSqr = GREAT;
forAll(pEdges, pEdgeI)
{
const edge& e = pp.edges()[pEdges[pEdgeI]];
label nbrPointI = e.otherVertex(pointI);
if
(
patchConstraints[nbrPointI].first() <= 1
&& allPatchConstraints[nbrPointI].first() > 1
)
{
const vector& nbrFeatVec =
allPatchConstraints[pointI].second();
if (mag(featVec&nbrFeatVec) > featureCos)
{
// nbrPointI attracted to sameish feature
// Note: also check on position.
scalar d2 = magSqr
(
allPatchAttraction[nbrPointI]
);
const point featPt =
pp.localPoints()[nbrPointI]
+ allPatchAttraction[nbrPointI];
const scalar cosAngle =
(featVec & (featPt-pt));
if (cosAngle > 0)
{
if (!hasPos && d2 < minPosDistSqr)
{
minPosDistSqr = d2;
bestPosPointI = nbrPointI;
}
}
else
{
if (!hasNeg && d2 < minNegDistSqr)
{
minNegDistSqr = d2;
bestNegPointI = nbrPointI;
}
}
}
}
}
if (bestPosPointI != -1)
{
// Use reconstructed-feature attraction. Use only
// part of it since not sure...
//const point& bestPt =
// pp.localPoints()[bestPosPointI];
//Pout<< "**Overriding point " << bestPt
// << " on reconstructed feature edge at "
// << allPatchAttraction[bestPosPointI]+bestPt
// << " to attracted-to-feature-edge." << endl;
patchAttraction[bestPosPointI] =
0.5*allPatchAttraction[bestPosPointI];
patchConstraints[bestPosPointI] =
allPatchConstraints[bestPosPointI];
nChanged++;
}
if (bestNegPointI != -1)
{
// Use reconstructed-feature attraction. Use only
// part of it since not sure...
//const point& bestPt =
// pp.localPoints()[bestNegPointI];
//Pout<< "**Overriding point " << bestPt
// << " on reconstructed feature edge at "
// << allPatchAttraction[bestNegPointI]+bestPt
// << " to attracted-to-feature-edge." << endl;
patchAttraction[bestNegPointI] =
0.5*allPatchAttraction[bestNegPointI];
patchConstraints[bestNegPointI] =
allPatchConstraints[bestNegPointI];
nChanged++;
}
}
}
}
reduce(nChanged, sumOp<label>());
Info<< "Stringing feature edges : changed " << nChanged << " points"
<< endl;
if (nChanged == 0)
{
break;
}
}
allPatchAttraction,
allPatchConstraints,
patchAttraction,
patchConstraints
);
// Avoid diagonal attraction

3
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
View File

@ -312,7 +312,8 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
{
Pout<< "Adding particle from point:" << pointI
<< " coord:" << featureMesh.points()[pointI]
<< " pEdges:" << pointEdges[pointI]
<< " since number of emanating edges:"
<< pointEdges[pointI].size()
<< endl;
}

2
src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermo.C
View File

@ -31,7 +31,7 @@ License
namespace Foam
{
defineTypeNameAndDebug(solidReactionThermo, 0);
defineRunTimeSelectionTable(solidReactionThermo, mesh);
defineRunTimeSelectionTable(solidReactionThermo, fvMesh);
defineRunTimeSelectionTable(solidReactionThermo, dictionary);
}

2
src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermo.H
View File

@ -66,7 +66,7 @@ public:
(
autoPtr,
solidReactionThermo,
mesh,
fvMesh,
(const fvMesh& mesh),
(mesh)
);

8
src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermoNew.C
View File

@ -53,16 +53,16 @@ Foam::autoPtr<Foam::solidReactionThermo> Foam::solidReactionThermo::New
Info<< "Selecting thermodynamics package " << modelType << endl;
meshConstructorTable::iterator cstrIter =
meshConstructorTablePtr_->find(modelType);
fvMeshConstructorTable::iterator cstrIter =
fvMeshConstructorTablePtr_->find(modelType);
if (cstrIter == meshConstructorTablePtr_->end())
if (cstrIter == fvMeshConstructorTablePtr_->end())
{
FatalErrorIn("solidReactionThermo::New(const fvMesh&)")
<< "Unknown solidReactionThermo type "
<< modelType << nl << nl
<< "Valid solidReactionThermo types:" << nl
<< meshConstructorTablePtr_->sortedToc() << nl
<< fvMeshConstructorTablePtr_->sortedToc() << nl
<< exit(FatalError);
}