mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
Working on feature point handling by properly constructing planes
This commit is contained in:
graham
@ -0,0 +1,178 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software; you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by the
|
||||
Free Software Foundation; either version 2 of the License, or (at your
|
||||
option) any later version.
|
||||
|
||||
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
||||
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM; if not, write to the Free Software Foundation,
|
||||
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
\*----------------------------------------------------------------------------*/
|
||||
|
||||
#include "CV3D.H"
|
||||
#include "plane.H"
|
||||
#include "triSurfaceTools.H"
|
||||
#include "mathematicalConstants.H"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
void Foam::CV3D::insertFeaturePoints()
|
||||
{
|
||||
const edgeList& edges = qSurf_.edges();
|
||||
const pointField& localPts = qSurf_.localPoints();
|
||||
|
||||
labelList featPoints(0);
|
||||
labelListList featPointFeatEdges(0);
|
||||
|
||||
qSurf_.extractFeatures3D
|
||||
(
|
||||
controls_.featAngle,
|
||||
featPoints,
|
||||
featPointFeatEdges
|
||||
);
|
||||
|
||||
forAll(featPoints, i)
|
||||
{
|
||||
label ptI = featPoints[i];
|
||||
const point& featPt = localPts[ptI];
|
||||
|
||||
const labelList& featEdges = featPointFeatEdges[i];
|
||||
|
||||
forAll(featEdges, fE)
|
||||
{
|
||||
label edgeI = featEdges[fE];
|
||||
const edge& featEdge = edges[edgeI];
|
||||
|
||||
// Check direction of edge, if the feature point is at the end()
|
||||
// the reverse direction.
|
||||
|
||||
scalar edgeDirection = 1.0;
|
||||
|
||||
if (ptI == featEdge.end())
|
||||
{
|
||||
edgeDirection = -1.0;
|
||||
}
|
||||
|
||||
point edgeLocalFeatPt = featPt
|
||||
+ 2.0*tols_.ppDist*edgeDirection
|
||||
* featEdge.vec(localPts)/featEdge.mag(localPts);
|
||||
|
||||
// Pick up the two faces adjacent to the feature edge
|
||||
const labelList& eFaces = qSurf_.edgeFaces()[edgeI];
|
||||
|
||||
label faceA = eFaces[0];
|
||||
vector nA = qSurf_.faceNormals()[faceA];
|
||||
|
||||
label faceB = eFaces[1];
|
||||
vector nB = qSurf_.faceNormals()[faceB];
|
||||
|
||||
// Intersect planes parallel to faceA and faceB offset by ppDist
|
||||
// and the plane defined by edgeLocalFeatPt and the edge vector.
|
||||
plane planeA(edgeLocalFeatPt - tols_.ppDist*nA, nA);
|
||||
plane planeB(edgeLocalFeatPt - tols_.ppDist*nB, nB);
|
||||
|
||||
plane planeF(edgeLocalFeatPt, featEdge.vec(localPts));
|
||||
|
||||
point refPt = planeF.planePlaneIntersect(planeA,planeB);
|
||||
|
||||
point faceAVert =
|
||||
localPts[triSurfaceTools::oppositeVertex(qSurf_, faceA, edgeI)];
|
||||
|
||||
// Determine convex or concave angle
|
||||
if (((faceAVert - edgeLocalFeatPt) & nB) < 0)
|
||||
{
|
||||
// Convex. So refPt will be inside domain and hence a master point
|
||||
|
||||
// Insert the master point refering the the first slave
|
||||
label masterPtIndex = insertPoint(refPt, number_of_vertices() + 1);
|
||||
|
||||
// Insert the slave points by reflecting refPt in both faces.
|
||||
// with each slave refering to the master
|
||||
|
||||
point reflectedA = refPt + 2*((edgeLocalFeatPt - refPt) & nA)*nA;
|
||||
insertPoint(reflectedA, masterPtIndex);
|
||||
|
||||
point reflectedB = refPt + 2*((edgeLocalFeatPt - refPt) & nB)*nB;
|
||||
insertPoint(reflectedB, masterPtIndex);
|
||||
}
|
||||
else
|
||||
{
|
||||
// Concave. master and reflected points inside the domain.
|
||||
// Generate reflected master to be outside.
|
||||
point reflMasterPt = refPt + 2*(edgeLocalFeatPt - refPt);
|
||||
|
||||
// Reflect refPt in both faces.
|
||||
point reflectedA =
|
||||
reflMasterPt + 2*((edgeLocalFeatPt - reflMasterPt) & nA)*nA;
|
||||
|
||||
point reflectedB =
|
||||
reflMasterPt + 2*((edgeLocalFeatPt - reflMasterPt) & nB)*nB;
|
||||
|
||||
scalar totalAngle =
|
||||
180*(mathematicalConstant::pi + acos(mag(nA & nB)))
|
||||
/mathematicalConstant::pi;
|
||||
|
||||
// Number of quadrants the angle should be split into
|
||||
int nQuads = int(totalAngle/controls_.maxQuadAngle) + 1;
|
||||
|
||||
// The number of additional master points needed to obtain the
|
||||
// required number of quadrants.
|
||||
int nAddPoints = min(max(nQuads - 2, 0), 2);
|
||||
|
||||
// index of reflMaster
|
||||
label reflectedMaster = number_of_vertices() + 2 + nAddPoints;
|
||||
|
||||
// Master A is inside.
|
||||
label reflectedAI = insertPoint(reflectedA, reflectedMaster);
|
||||
|
||||
// Master B is inside.
|
||||
insertPoint(reflectedB, reflectedMaster);
|
||||
|
||||
if (nAddPoints == 1)
|
||||
{
|
||||
// One additinal point is the reflection of the slave point,
|
||||
// i.e. the original reference point
|
||||
insertPoint(refPt, reflectedMaster);
|
||||
}
|
||||
else if (nAddPoints == 2)
|
||||
{
|
||||
point reflectedAa = refPt
|
||||
- ((edgeLocalFeatPt - reflMasterPt) & nB)*nB;
|
||||
insertPoint(reflectedAa, reflectedMaster);
|
||||
|
||||
point reflectedBb = refPt
|
||||
- ((edgeLocalFeatPt - reflMasterPt) & nA)*nA;
|
||||
insertPoint(reflectedBb, reflectedMaster);
|
||||
}
|
||||
|
||||
// Slave is outside.
|
||||
insertPoint(reflMasterPt, reflectedAI);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (controls_.writeFeatureTriangulation)
|
||||
{
|
||||
writePoints("feat_allPoints.obj", false);
|
||||
writePoints("feat_points.obj", true);
|
||||
writeTriangles("feat_triangles.obj", true);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,296 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software; you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by the
|
||||
Free Software Foundation; either version 2 of the License, or (at your
|
||||
option) any later version.
|
||||
|
||||
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
||||
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM; if not, write to the Free Software Foundation,
|
||||
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
\*----------------------------------------------------------------------------*/
|
||||
|
||||
#include "querySurface.H"
|
||||
#include "mathematicalConstants.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::querySurface::querySurface(const fileName& surfaceFileName)
|
||||
:
|
||||
triSurface(surfaceFileName),
|
||||
rndGen_(12345),
|
||||
bb_(localPoints()),
|
||||
tree_
|
||||
(
|
||||
treeDataTriSurface(*this),
|
||||
bb_.extend(rndGen_, 1e-3), // slightly randomize bb
|
||||
8, // maxLevel
|
||||
4, //10, // leafsize
|
||||
10.0 //3.0 // duplicity
|
||||
)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::querySurface::~querySurface()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::labelList Foam::querySurface::extractFeatures2D
|
||||
(
|
||||
const scalar featAngle
|
||||
) const
|
||||
{
|
||||
scalar featCos = cos(mathematicalConstant::pi*featAngle/180.0);
|
||||
|
||||
const labelListList& edgeFaces = this->edgeFaces();
|
||||
const pointField& localPoints = this->localPoints();
|
||||
const edgeList& edges = this->edges();
|
||||
const vectorField& faceNormals = this->faceNormals();
|
||||
|
||||
DynamicList<label> featEdges(edges.size());
|
||||
|
||||
forAll(edgeFaces, edgeI)
|
||||
{
|
||||
const edge& e = edges[edgeI];
|
||||
|
||||
if (magSqr(e.vec(localPoints) & vector(1,1,0)) < SMALL)
|
||||
{
|
||||
const labelList& eFaces = edgeFaces[edgeI];
|
||||
|
||||
if
|
||||
(
|
||||
eFaces.size() == 2
|
||||
&& (faceNormals[eFaces[0]] & faceNormals[eFaces[1]]) < featCos
|
||||
)
|
||||
{
|
||||
featEdges.append(edgeI);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return featEdges.shrink();
|
||||
}
|
||||
|
||||
void Foam::querySurface::extractFeatures3D
|
||||
(
|
||||
const scalar featAngle,
|
||||
labelList& featPoints,
|
||||
labelListList& featPointFeatEdges
|
||||
) const
|
||||
{
|
||||
scalar featCos = cos(mathematicalConstant::pi*featAngle/180.0);
|
||||
|
||||
const labelListList& edgeFaces = this->edgeFaces();
|
||||
const labelListList& pointEdges = this->pointEdges();
|
||||
const pointField& localPoints = this->localPoints();
|
||||
// const edgeList& edges = this->edges();
|
||||
const vectorField& faceNormals = this->faceNormals();
|
||||
|
||||
DynamicList<label> tempFeatPoints(localPoints.size());
|
||||
|
||||
DynamicList<DynamicList<label> > tempFeatPointFeatEdges(localPoints.size());
|
||||
|
||||
forAll(pointEdges, ptI)
|
||||
{
|
||||
// const point& p = localPoints[ptI];
|
||||
|
||||
// const vector& pNormal = pointNormals[ptI];
|
||||
|
||||
const labelList& pEdges = pointEdges[ptI];
|
||||
|
||||
DynamicList<label> tempFeatEdges(0);
|
||||
|
||||
forAll(pEdges, pE)
|
||||
{
|
||||
label edgeI = pEdges[pE];
|
||||
|
||||
const labelList& eFaces = edgeFaces[edgeI];
|
||||
|
||||
if
|
||||
(
|
||||
eFaces.size() == 2
|
||||
&& (faceNormals[eFaces[0]] & faceNormals[eFaces[1]]) < featCos
|
||||
)
|
||||
{
|
||||
tempFeatEdges.append(edgeI);
|
||||
}
|
||||
}
|
||||
|
||||
tempFeatEdges.shrink();
|
||||
|
||||
if(tempFeatEdges.size())
|
||||
{
|
||||
tempFeatPoints.append(ptI);
|
||||
|
||||
tempFeatPointFeatEdges.append(tempFeatEdges);
|
||||
}
|
||||
}
|
||||
|
||||
featPoints.transfer(tempFeatPoints.shrink());
|
||||
|
||||
tempFeatPointFeatEdges.shrink();
|
||||
|
||||
featPointFeatEdges.setSize(tempFeatPointFeatEdges.size());
|
||||
|
||||
forAll(featPointFeatEdges, fPFE)
|
||||
{
|
||||
featPointFeatEdges[fPFE].transfer(tempFeatPointFeatEdges[fPFE].shrink());
|
||||
}
|
||||
}
|
||||
|
||||
Foam::indexedOctree<Foam::treeDataTriSurface>::volumeType
|
||||
Foam::querySurface::insideOutside
|
||||
(
|
||||
const scalar searchSpan2,
|
||||
const point& pt
|
||||
) const
|
||||
{
|
||||
if (!bb_.contains(pt))
|
||||
{
|
||||
return indexedOctree<treeDataTriSurface>::OUTSIDE;
|
||||
}
|
||||
|
||||
pointIndexHit pHit = tree_.findNearest(pt, searchSpan2);
|
||||
|
||||
if (!pHit.hit())
|
||||
{
|
||||
return tree_.getVolumeType(pt);
|
||||
}
|
||||
else
|
||||
{
|
||||
return indexedOctree<treeDataTriSurface>::MIXED;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Check if point is inside surface
|
||||
bool Foam::querySurface::inside(const point& pt) const
|
||||
{
|
||||
if (!bb_.contains(pt))
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
return
|
||||
(
|
||||
tree_.getVolumeType(pt) == indexedOctree<treeDataTriSurface>::INSIDE
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// Check if point is outside surface
|
||||
bool Foam::querySurface::outside(const point& pt) const
|
||||
{
|
||||
if (!bb_.contains(pt))
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
return
|
||||
(
|
||||
tree_.getVolumeType(pt) == indexedOctree<treeDataTriSurface>::OUTSIDE
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// Check if point is inside surface by at least dist2
|
||||
bool Foam::querySurface::wellInside(const point& pt, const scalar dist2) const
|
||||
{
|
||||
if (!bb_.contains(pt))
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
pointIndexHit pHit = tree_.findNearest(pt, dist2);
|
||||
|
||||
if (pHit.hit())
|
||||
{
|
||||
return false;
|
||||
}
|
||||
else
|
||||
{
|
||||
return
|
||||
tree_.getVolumeType(pt)
|
||||
== indexedOctree<treeDataTriSurface>::INSIDE;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Check if point is outside surface by at least dist2
|
||||
bool Foam::querySurface::wellOutside(const point& pt, const scalar dist2) const
|
||||
{
|
||||
if (!bb_.contains(pt))
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
pointIndexHit pHit = tree_.findNearest(pt, dist2);
|
||||
|
||||
if (pHit.hit())
|
||||
{
|
||||
return false;
|
||||
}
|
||||
else
|
||||
{
|
||||
return
|
||||
tree_.getVolumeType(pt)
|
||||
== indexedOctree<treeDataTriSurface>::OUTSIDE;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::querySurface::writeTreeOBJ() const
|
||||
{
|
||||
OFstream str("tree.obj");
|
||||
label vertI = 0;
|
||||
|
||||
const List<indexedOctree<treeDataTriSurface>::node>& nodes = tree_.nodes();
|
||||
|
||||
forAll(nodes, nodeI)
|
||||
{
|
||||
const indexedOctree<treeDataTriSurface>::node& nod = nodes[nodeI];
|
||||
|
||||
const treeBoundBox& bb = nod.bb_;
|
||||
|
||||
const pointField points(bb.points());
|
||||
|
||||
label startVertI = vertI;
|
||||
|
||||
forAll(points, i)
|
||||
{
|
||||
meshTools::writeOBJ(str, points[i]);
|
||||
vertI++;
|
||||
}
|
||||
|
||||
const edgeList edges(treeBoundBox::edges);
|
||||
|
||||
forAll(edges, i)
|
||||
{
|
||||
const edge& e = edges[i];
|
||||
|
||||
str << "l " << e[0]+startVertI+1 << ' ' << e[1]+startVertI+1
|
||||
<< nl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,96 @@
|
||||
IOobject io
|
||||
(
|
||||
Foam::polyMesh::defaultRegion,
|
||||
runTime.constant(),
|
||||
runTime,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
);
|
||||
|
||||
polyMesh pMesh(io);
|
||||
|
||||
IOobject io2
|
||||
(
|
||||
Foam::polyMesh::defaultRegion,
|
||||
runTime.timeName(),
|
||||
runTime,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
);
|
||||
|
||||
polyMesh pMesh2(io2, pointField(0), faceList(0), cellList(0));
|
||||
|
||||
polyTopoChange meshCreator(pMesh.boundaryMesh().size());
|
||||
|
||||
meshCreator.addPoint(point(-1,-1,-1), -1, -1, true);
|
||||
meshCreator.addPoint(point(1,-1,-1), -1, -1, true);
|
||||
meshCreator.addPoint(point(1,1,-1), -1, -1, true);
|
||||
meshCreator.addPoint(point(-1,1,-1), -1, -1, true);
|
||||
meshCreator.addPoint(point(-1,-1,1), -1, -1, true);
|
||||
meshCreator.addPoint(point(1,-1,1), -1, -1, true);
|
||||
meshCreator.addPoint(point(1,1,1), -1, -1, true);
|
||||
meshCreator.addPoint(point(-1,1,1), -1, -1, true);
|
||||
|
||||
meshCreator.addCell(-1, -1, -1, -1, -1);
|
||||
|
||||
labelList faceConstruct(4);
|
||||
|
||||
faceConstruct[0] = 1;
|
||||
faceConstruct[1] = 2;
|
||||
faceConstruct[2] = 6;
|
||||
faceConstruct[3] = 5;
|
||||
|
||||
meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
|
||||
|
||||
faceConstruct[0] = 0;
|
||||
faceConstruct[1] = 4;
|
||||
faceConstruct[2] = 7;
|
||||
faceConstruct[3] = 3;
|
||||
|
||||
meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
|
||||
|
||||
faceConstruct[0] = 2;
|
||||
faceConstruct[1] = 3;
|
||||
faceConstruct[2] = 7;
|
||||
faceConstruct[3] = 6;
|
||||
|
||||
meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
|
||||
|
||||
faceConstruct[0] = 0;
|
||||
faceConstruct[1] = 1;
|
||||
faceConstruct[2] = 5;
|
||||
faceConstruct[3] = 4;
|
||||
|
||||
meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
|
||||
|
||||
faceConstruct[0] = 0;
|
||||
faceConstruct[1] = 3;
|
||||
faceConstruct[2] = 2;
|
||||
faceConstruct[3] = 1;
|
||||
|
||||
meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
|
||||
|
||||
faceConstruct[0] = 4;
|
||||
faceConstruct[1] = 5;
|
||||
faceConstruct[2] = 6;
|
||||
faceConstruct[3] = 7;
|
||||
|
||||
meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
|
||||
|
||||
Info<< meshCreator.points() << endl;
|
||||
Info<< meshCreator.faces() << endl;
|
||||
Info<< meshCreator.faceOwner() << endl;
|
||||
Info<< meshCreator.faceNeighbour() << endl;
|
||||
|
||||
// calcDualMesh(meshCreator);
|
||||
|
||||
autoPtr<mapPolyMesh> map = meshCreator.changeMesh(pMesh, false);
|
||||
|
||||
pMesh.updateMesh(map);
|
||||
|
||||
if (!pMesh.write())
|
||||
{
|
||||
FatalErrorIn("CV3D::writeMesh(const Time& runTime)")
|
||||
<< "Failed writing polyMesh."
|
||||
<< exit(FatalError);
|
||||
}
|
||||
@ -391,8 +391,6 @@ void Foam::CV3D::calcDualMesh
|
||||
|
||||
patchOwners[p].shrink();
|
||||
|
||||
// TODO SORT BOUNDARY FACES AND NEIGHBOURS?
|
||||
|
||||
patchSizes[p] = patchFaces[p].size();
|
||||
|
||||
patchStarts[p] = nInternalFaces + nBoundaryFaces;
|
||||
|
||||
@ -39,13 +39,17 @@ void Foam::CV3D::insertFeaturePoints()
|
||||
labelList featPoints(0);
|
||||
labelListList featPointFeatEdges(0);
|
||||
|
||||
qSurf_.extractFeatures3D
|
||||
qSurf_.extractFeatures
|
||||
(
|
||||
controls_.featAngle,
|
||||
featPoints,
|
||||
featPointFeatEdges
|
||||
);
|
||||
|
||||
scalar planeErrorAngle = 0.1*controls_.featAngle;
|
||||
|
||||
scalar planeErrorAngleCos = cos(mathematicalConstant::pi*planeErrorAngle/180.0);
|
||||
|
||||
forAll(featPoints, i)
|
||||
{
|
||||
label ptI = featPoints[i];
|
||||
@ -53,25 +57,17 @@ void Foam::CV3D::insertFeaturePoints()
|
||||
|
||||
const labelList& featEdges = featPointFeatEdges[i];
|
||||
|
||||
DynamicList<label> convexEdges(0);
|
||||
DynamicList<label> concaveEdges(0);
|
||||
|
||||
List<Pair<vector> > planeNormalPairs(featEdges.size());
|
||||
|
||||
// Classify the edges around the feature point.
|
||||
forAll(featEdges, fE)
|
||||
{
|
||||
label edgeI = featEdges[fE];
|
||||
const edge& featEdge = edges[edgeI];
|
||||
|
||||
// Check direction of edge, if the feature point is at the end()
|
||||
// the reverse direction.
|
||||
|
||||
scalar edgeDirection = 1.0;
|
||||
|
||||
if (ptI == featEdge.end())
|
||||
{
|
||||
edgeDirection = -1.0;
|
||||
}
|
||||
|
||||
point edgeLocalFeatPt = featPt
|
||||
+ 2.0*tols_.ppDist*edgeDirection
|
||||
* featEdge.vec(localPts)/featEdge.mag(localPts);
|
||||
|
||||
// Pick up the two faces adjacent to the feature edge
|
||||
const labelList& eFaces = qSurf_.edgeFaces()[edgeI];
|
||||
|
||||
@ -81,91 +77,108 @@ void Foam::CV3D::insertFeaturePoints()
|
||||
label faceB = eFaces[1];
|
||||
vector nB = qSurf_.faceNormals()[faceB];
|
||||
|
||||
// Intersect planes parallel to faceA and faceB offset by ppDist
|
||||
// and the plane defined by edgeLocalFeatPt and the edge vector.
|
||||
plane planeA(edgeLocalFeatPt - tols_.ppDist*nA, nA);
|
||||
plane planeB(edgeLocalFeatPt - tols_.ppDist*nB, nB);
|
||||
|
||||
plane planeF(edgeLocalFeatPt, featEdge.vec(localPts));
|
||||
|
||||
point refPt = planeF.planePlaneIntersect(planeA,planeB);
|
||||
|
||||
point faceAVert =
|
||||
localPts[triSurfaceTools::oppositeVertex(qSurf_, faceA, edgeI)];
|
||||
|
||||
// Determine convex or concave angle
|
||||
if (((faceAVert - edgeLocalFeatPt) & nB) < 0)
|
||||
if (((faceAVert - featPt) & nB) < 0)
|
||||
{
|
||||
// Convex. So refPt will be inside domain and hence a master point
|
||||
|
||||
// Insert the master point refering the the first slave
|
||||
label masterPtIndex = insertPoint(refPt, number_of_vertices() + 1);
|
||||
|
||||
// Insert the slave points by reflecting refPt in both faces.
|
||||
// with each slave refering to the master
|
||||
|
||||
point reflectedA = refPt + 2*((edgeLocalFeatPt - refPt) & nA)*nA;
|
||||
insertPoint(reflectedA, masterPtIndex);
|
||||
|
||||
point reflectedB = refPt + 2*((edgeLocalFeatPt - refPt) & nB)*nB;
|
||||
insertPoint(reflectedB, masterPtIndex);
|
||||
// Convex feature edge
|
||||
convexEdges.append(edgeI);
|
||||
}
|
||||
else
|
||||
{
|
||||
// Concave. master and reflected points inside the domain.
|
||||
// Generate reflected master to be outside.
|
||||
point reflMasterPt = refPt + 2*(edgeLocalFeatPt - refPt);
|
||||
|
||||
// Reflect refPt in both faces.
|
||||
point reflectedA =
|
||||
reflMasterPt + 2*((edgeLocalFeatPt - reflMasterPt) & nA)*nA;
|
||||
|
||||
point reflectedB =
|
||||
reflMasterPt + 2*((edgeLocalFeatPt - reflMasterPt) & nB)*nB;
|
||||
|
||||
scalar totalAngle =
|
||||
180*(mathematicalConstant::pi + acos(mag(nA & nB)))
|
||||
/mathematicalConstant::pi;
|
||||
|
||||
// Number of quadrants the angle should be split into
|
||||
int nQuads = int(totalAngle/controls_.maxQuadAngle) + 1;
|
||||
|
||||
// The number of additional master points needed to obtain the
|
||||
// required number of quadrants.
|
||||
int nAddPoints = min(max(nQuads - 2, 0), 2);
|
||||
|
||||
// index of reflMaster
|
||||
label reflectedMaster = number_of_vertices() + 2 + nAddPoints;
|
||||
|
||||
// Master A is inside.
|
||||
label reflectedAI = insertPoint(reflectedA, reflectedMaster);
|
||||
|
||||
// Master B is inside.
|
||||
insertPoint(reflectedB, reflectedMaster);
|
||||
|
||||
if (nAddPoints == 1)
|
||||
{
|
||||
// One additinal point is the reflection of the slave point,
|
||||
// i.e. the original reference point
|
||||
insertPoint(refPt, reflectedMaster);
|
||||
}
|
||||
else if (nAddPoints == 2)
|
||||
{
|
||||
point reflectedAa = refPt
|
||||
- ((edgeLocalFeatPt - reflMasterPt) & nB)*nB;
|
||||
insertPoint(reflectedAa, reflectedMaster);
|
||||
|
||||
point reflectedBb = refPt
|
||||
- ((edgeLocalFeatPt - reflMasterPt) & nA)*nA;
|
||||
insertPoint(reflectedBb, reflectedMaster);
|
||||
}
|
||||
|
||||
// Slave is outside.
|
||||
insertPoint(reflMasterPt, reflectedAI);
|
||||
// Concave feature edge
|
||||
concaveEdges.append(edgeI);
|
||||
}
|
||||
|
||||
// Identify the normals of the faces attached to feature edges to
|
||||
// identify the unique planes to be reconstructed. The triangulated
|
||||
// surface will usually not mean that two feature edges that should
|
||||
// bound a plane are attached to the same face.
|
||||
|
||||
planeNormalPairs[fE].first() = nA;
|
||||
planeNormalPairs[fE].second() = nB;
|
||||
}
|
||||
|
||||
convexEdges.shrink();
|
||||
concaveEdges.shrink();
|
||||
|
||||
Info<< nl << convexEdges
|
||||
<< nl << concaveEdges
|
||||
<< nl << planeNormalPairs
|
||||
<< endl;
|
||||
|
||||
List<vector> uniquePlaneNormals(featEdges.size());
|
||||
label uniquePlaneNormalI = 0;
|
||||
|
||||
List<Pair<bool> > planeMatchedStatus(featEdges.size(), Pair<bool>(false,false));
|
||||
|
||||
Info<< planeMatchedStatus << endl;
|
||||
|
||||
// Examine the plane normals to identify unique planes.
|
||||
forAll(planeNormalPairs, nA)
|
||||
{
|
||||
const Pair<vector>& normalPairA = planeNormalPairs[nA];
|
||||
|
||||
if (!planeMatchedStatus[nA].first())
|
||||
{
|
||||
const vector& nAf = normalPairA.first();
|
||||
|
||||
scalar minNormalDotProduct = 1 + SMALL;
|
||||
|
||||
forAll(planeNormalPairs, nB)
|
||||
{
|
||||
if (nA == nB)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
const Pair<vector>& normalPairB = planeNormalPairs[nB];
|
||||
|
||||
if (!planeMatchedStatus[nB].first())
|
||||
{
|
||||
const vector& nBf = normalPairB.first();
|
||||
|
||||
scalar normalDotProduct = nAf & nBf;
|
||||
|
||||
if (normalDotProduct < minNormalDotProduct)
|
||||
{
|
||||
minNormalDotProduct = normalDotProduct;
|
||||
}
|
||||
}
|
||||
|
||||
if (!planeMatchedStatus[nB].second())
|
||||
{
|
||||
const vector& nBs = normalPairA.second();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (!planeMatchedStatus[nA].second())
|
||||
{
|
||||
const vector& nAs = normalPairA.second();
|
||||
}
|
||||
|
||||
// if (minNormalDotProduct > planeErrorAngleCos)
|
||||
// {
|
||||
// FatalErrorIn("insertFeaturePoints")
|
||||
// << "Could not find unique planes matching feature edges "
|
||||
// << "at point located at "
|
||||
// << featPt << nl
|
||||
// << exit(FatalError);
|
||||
// }
|
||||
|
||||
// const vector& matchingPNB = planeNormalsB[matchingPB];
|
||||
|
||||
// uniquePlaneNormals[pA] =
|
||||
// (pNA + matchingPNB)/(mag(pNA + matchingPNB) + VSMALL);
|
||||
}
|
||||
|
||||
Info<< uniquePlaneNormals << endl;
|
||||
}
|
||||
|
||||
|
||||
if (controls_.writeFeatureTriangulation)
|
||||
{
|
||||
writePoints("feat_allPoints.obj", false);
|
||||
|
||||
@ -53,43 +53,7 @@ Foam::querySurface::~querySurface()
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::labelList Foam::querySurface::extractFeatures2D
|
||||
(
|
||||
const scalar featAngle
|
||||
) const
|
||||
{
|
||||
scalar featCos = cos(mathematicalConstant::pi*featAngle/180.0);
|
||||
|
||||
const labelListList& edgeFaces = this->edgeFaces();
|
||||
const pointField& localPoints = this->localPoints();
|
||||
const edgeList& edges = this->edges();
|
||||
const vectorField& faceNormals = this->faceNormals();
|
||||
|
||||
DynamicList<label> featEdges(edges.size());
|
||||
|
||||
forAll(edgeFaces, edgeI)
|
||||
{
|
||||
const edge& e = edges[edgeI];
|
||||
|
||||
if (magSqr(e.vec(localPoints) & vector(1,1,0)) < SMALL)
|
||||
{
|
||||
const labelList& eFaces = edgeFaces[edgeI];
|
||||
|
||||
if
|
||||
(
|
||||
eFaces.size() == 2
|
||||
&& (faceNormals[eFaces[0]] & faceNormals[eFaces[1]]) < featCos
|
||||
)
|
||||
{
|
||||
featEdges.append(edgeI);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return featEdges.shrink();
|
||||
}
|
||||
|
||||
void Foam::querySurface::extractFeatures3D
|
||||
void Foam::querySurface::extractFeatures
|
||||
(
|
||||
const scalar featAngle,
|
||||
labelList& featPoints,
|
||||
@ -136,7 +100,7 @@ void Foam::querySurface::extractFeatures3D
|
||||
|
||||
tempFeatEdges.shrink();
|
||||
|
||||
if(tempFeatEdges.size())
|
||||
if(tempFeatEdges.size() > 2)
|
||||
{
|
||||
tempFeatPoints.append(ptI);
|
||||
|
||||
|
||||
@ -104,13 +104,9 @@ public:
|
||||
|
||||
// Query
|
||||
|
||||
//- Extract feature edges/points(2D)
|
||||
//- Extract feature edges/points
|
||||
// using the given feature angle in deg
|
||||
labelList extractFeatures2D(const scalar featAngle) const;
|
||||
|
||||
//- Extract feature edges/points(3D)
|
||||
// using the given feature angle in deg
|
||||
void extractFeatures3D
|
||||
void extractFeatures
|
||||
(
|
||||
const scalar featAngle,
|
||||
labelList& featPoints,
|
||||
|
||||
@ -1,56 +0,0 @@
|
||||
/*--------------------------------*- C++ -*----------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.5 |
|
||||
| \\ / A nd | Web: http://www.OpenFOAM.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class volScalarField;
|
||||
object G;
|
||||
}
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
dimensions [1 0 -3 0 0 0 0];
|
||||
|
||||
internalField uniform 0;
|
||||
|
||||
boundaryField
|
||||
{
|
||||
floor
|
||||
{
|
||||
type MarshakRadiation;
|
||||
T T;
|
||||
emissivity 1;
|
||||
value uniform 0;
|
||||
}
|
||||
|
||||
fixedWalls
|
||||
{
|
||||
type MarshakRadiation;
|
||||
T T;
|
||||
emissivity 1;
|
||||
value uniform 0;
|
||||
}
|
||||
|
||||
ceiling
|
||||
{
|
||||
type MarshakRadiation;
|
||||
T T;
|
||||
emissivity 1;
|
||||
value uniform 0;
|
||||
}
|
||||
|
||||
box
|
||||
{
|
||||
type MarshakRadiation;
|
||||
T T;
|
||||
emissivity 1;
|
||||
value uniform 0;
|
||||
}
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,47 +0,0 @@
|
||||
/*--------------------------------*- C++ -*----------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.5 |
|
||||
| \\ / A nd | Web: http://www.OpenFOAM.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class volScalarField;
|
||||
object T;
|
||||
}
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
dimensions [0 0 0 1 0 0 0];
|
||||
|
||||
internalField uniform 300;
|
||||
|
||||
boundaryField
|
||||
{
|
||||
floor
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform 300.0;
|
||||
}
|
||||
|
||||
ceiling
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform 300.0;
|
||||
}
|
||||
|
||||
fixedWalls
|
||||
{
|
||||
type zeroGradient;
|
||||
}
|
||||
|
||||
box
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform 500.0;
|
||||
}
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,48 +0,0 @@
|
||||
/*--------------------------------*- C++ -*----------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.5 |
|
||||
| \\ / A nd | Web: http://www.OpenFOAM.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class volVectorField;
|
||||
object U;
|
||||
}
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
dimensions [0 1 -1 0 0 0 0];
|
||||
|
||||
internalField uniform (0 0 0);
|
||||
|
||||
boundaryField
|
||||
{
|
||||
floor
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform (0 0 0);
|
||||
}
|
||||
|
||||
ceiling
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform (0 0 0);
|
||||
}
|
||||
|
||||
fixedWalls
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform (0 0 0);
|
||||
}
|
||||
|
||||
box
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform (0 0 0);
|
||||
}
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,44 +0,0 @@
|
||||
/*--------------------------------*- C++ -*----------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.5 |
|
||||
| \\ / A nd | Web: http://www.OpenFOAM.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class volScalarField;
|
||||
object epsilon;
|
||||
}
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
dimensions [0 2 -3 0 0 0 0];
|
||||
|
||||
internalField uniform 0.01;
|
||||
|
||||
boundaryField
|
||||
{
|
||||
floor
|
||||
{
|
||||
type zeroGradient;
|
||||
}
|
||||
|
||||
ceiling
|
||||
{
|
||||
type zeroGradient;
|
||||
}
|
||||
|
||||
fixedWalls
|
||||
{
|
||||
type zeroGradient;
|
||||
}
|
||||
|
||||
box
|
||||
{
|
||||
type zeroGradient;
|
||||
}
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,44 +0,0 @@
|
||||
/*--------------------------------*- C++ -*----------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.5 |
|
||||
| \\ / A nd | Web: http://www.OpenFOAM.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class volScalarField;
|
||||
object k;
|
||||
}
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
dimensions [0 2 -2 0 0 0 0];
|
||||
|
||||
internalField uniform 0.1;
|
||||
|
||||
boundaryField
|
||||
{
|
||||
floor
|
||||
{
|
||||
type zeroGradient;
|
||||
}
|
||||
|
||||
ceiling
|
||||
{
|
||||
type zeroGradient;
|
||||
}
|
||||
|
||||
fixedWalls
|
||||
{
|
||||
type zeroGradient;
|
||||
}
|
||||
|
||||
box
|
||||
{
|
||||
type zeroGradient;
|
||||
}
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,48 +0,0 @@
|
||||
/*--------------------------------*- C++ -*----------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.5 |
|
||||
| \\ / A nd | Web: http://www.OpenFOAM.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class volScalarField;
|
||||
object p;
|
||||
}
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
dimensions [1 -1 -2 0 0 0 0];
|
||||
|
||||
internalField uniform 100000;
|
||||
|
||||
boundaryField
|
||||
{
|
||||
floor
|
||||
{
|
||||
type calculated;
|
||||
value uniform 100000;
|
||||
}
|
||||
|
||||
ceiling
|
||||
{
|
||||
type calculated;
|
||||
value uniform 100000;
|
||||
}
|
||||
|
||||
fixedWalls
|
||||
{
|
||||
type calculated;
|
||||
value uniform 100000;
|
||||
}
|
||||
|
||||
box
|
||||
{
|
||||
type calculated;
|
||||
value uniform 100000;
|
||||
}
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,48 +0,0 @@
|
||||
/*--------------------------------*- C++ -*----------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.5 |
|
||||
| \\ / A nd | Web: http://www.OpenFOAM.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class volScalarField;
|
||||
object pd;
|
||||
}
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
dimensions [1 -1 -2 0 0 0 0];
|
||||
|
||||
internalField uniform 0;
|
||||
|
||||
boundaryField
|
||||
{
|
||||
floor
|
||||
{
|
||||
type fixedFluxBuoyantPressure;
|
||||
value uniform 0;
|
||||
}
|
||||
|
||||
ceiling
|
||||
{
|
||||
type fixedFluxBuoyantPressure;
|
||||
value uniform 0;
|
||||
}
|
||||
|
||||
fixedWalls
|
||||
{
|
||||
type fixedFluxBuoyantPressure;
|
||||
value uniform 0;
|
||||
}
|
||||
|
||||
box
|
||||
{
|
||||
type fixedFluxBuoyantPressure;
|
||||
value uniform 0;
|
||||
}
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,7 +1,7 @@
|
||||
/*--------------------------------*- C++ -*----------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: 1.5 |
|
||||
| \\ / O peration | Version: dev |
|
||||
| \\ / A nd | Web: http://www.OpenFOAM.org |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
@ -10,6 +10,7 @@ FoamFile
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class polyBoundaryMesh;
|
||||
location "constant/polyMesh";
|
||||
object boundary;
|
||||
}
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
Reference in New Issue
Block a user