From 176a0b250727f0d9130d47ab7a476e731174c806 Mon Sep 17 00:00:00 2001
From: graham <graham.macpherson@strath.ac.uk>
Date: Thu, 31 Jul 2008 16:22:14 +0100
Subject: [PATCH] Working on feature point handling by properly constructing
planes
---
.../insertFeaturePoints.C | 178 +++++++++++
.../querySurface.C | 296 ++++++++++++++++++
.../hardCodedSimpleCubeForPolyTopoChange.H | 96 ++++++
.../mesh/generation/CV3DMesher/calcDualMesh.C | 2 -
.../CV3DMesher/insertFeaturePoints.C | 163 +++++-----
.../mesh/generation/CV3DMesher/querySurface.C | 40 +--
.../mesh/generation/CV3DMesher/querySurface.H | 8 +-
.../hotRadiationRoom/0/G | 56 ----
.../hotRadiationRoom/0/T | 47 ---
.../hotRadiationRoom/0/U | 48 ---
.../hotRadiationRoom/0/epsilon | 44 ---
.../hotRadiationRoom/0/k | 44 ---
.../hotRadiationRoom/0/p | 48 ---
.../hotRadiationRoom/0/pd | 48 ---
.../constant/polyMesh/boundary | 3 +-
15 files changed, 664 insertions(+), 457 deletions(-)
create mode 100644 applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/insertFeaturePoints.C
create mode 100644 applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/querySurface.C
create mode 100644 applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/hardCodedSimpleCubeForPolyTopoChange.H
delete mode 100644 tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/G
delete mode 100644 tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/T
delete mode 100644 tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/U
delete mode 100644 tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/epsilon
delete mode 100644 tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/k
delete mode 100644 tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p
delete mode 100644 tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/pd
diff --git a/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/insertFeaturePoints.C b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/insertFeaturePoints.C
new file mode 100644
index 00000000000..18f3db880af
--- /dev/null
+++ b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/insertFeaturePoints.C
@@ -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);
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/querySurface.C b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/querySurface.C
new file mode 100644
index 00000000000..311c6350c11
--- /dev/null
+++ b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/querySurface.C
@@ -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;
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/hardCodedSimpleCubeForPolyTopoChange.H b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/hardCodedSimpleCubeForPolyTopoChange.H
new file mode 100644
index 00000000000..8fd0e855ff7
--- /dev/null
+++ b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/hardCodedSimpleCubeForPolyTopoChange.H
@@ -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);
+ }
\ No newline at end of file
diff --git a/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C b/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
index 74c66438f02..018abb495d5 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
@@ -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;
diff --git a/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C b/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
index 18f3db880af..af97d348cfc 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
@@ -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);
+ // Concave feature edge
+ concaveEdges.append(edgeI);
+ }
- // Reflect refPt in both faces.
- point reflectedA =
- reflMasterPt + 2*((edgeLocalFeatPt - reflMasterPt) & nA)*nA;
+ // 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.
- point reflectedB =
- reflMasterPt + 2*((edgeLocalFeatPt - reflMasterPt) & nB)*nB;
+ planeNormalPairs[fE].first() = nA;
+ planeNormalPairs[fE].second() = nB;
+ }
- scalar totalAngle =
- 180*(mathematicalConstant::pi + acos(mag(nA & nB)))
- /mathematicalConstant::pi;
+ convexEdges.shrink();
+ concaveEdges.shrink();
- // Number of quadrants the angle should be split into
- int nQuads = int(totalAngle/controls_.maxQuadAngle) + 1;
+ Info<< nl << convexEdges
+ << nl << concaveEdges
+ << nl << planeNormalPairs
+ << endl;
- // The number of additional master points needed to obtain the
- // required number of quadrants.
- int nAddPoints = min(max(nQuads - 2, 0), 2);
+ List<vector> uniquePlaneNormals(featEdges.size());
+ label uniquePlaneNormalI = 0;
- // index of reflMaster
- label reflectedMaster = number_of_vertices() + 2 + nAddPoints;
+ List<Pair<bool> > planeMatchedStatus(featEdges.size(), Pair<bool>(false,false));
- // Master A is inside.
- label reflectedAI = insertPoint(reflectedA, reflectedMaster);
+ Info<< planeMatchedStatus << endl;
- // Master B is inside.
- insertPoint(reflectedB, reflectedMaster);
+ // Examine the plane normals to identify unique planes.
+ forAll(planeNormalPairs, nA)
+ {
+ const Pair<vector>& normalPairA = planeNormalPairs[nA];
- 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)
+ if (!planeMatchedStatus[nA].first())
+ {
+ const vector& nAf = normalPairA.first();
+
+ scalar minNormalDotProduct = 1 + SMALL;
+
+ forAll(planeNormalPairs, nB)
{
- point reflectedAa = refPt
- - ((edgeLocalFeatPt - reflMasterPt) & nB)*nB;
- insertPoint(reflectedAa, reflectedMaster);
+ if (nA == nB)
+ {
+ continue;
+ }
+
+ const Pair<vector>& normalPairB = planeNormalPairs[nB];
+
+ if (!planeMatchedStatus[nB].first())
+ {
+ const vector& nBf = normalPairB.first();
- point reflectedBb = refPt
- - ((edgeLocalFeatPt - reflMasterPt) & nA)*nA;
- insertPoint(reflectedBb, reflectedMaster);
+ scalar normalDotProduct = nAf & nBf;
+
+ if (normalDotProduct < minNormalDotProduct)
+ {
+ minNormalDotProduct = normalDotProduct;
+ }
+ }
+
+ if (!planeMatchedStatus[nB].second())
+ {
+ const vector& nBs = normalPairA.second();
+ }
}
+ }
- // Slave is outside.
- insertPoint(reflMasterPt, reflectedAI);
+ 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);
diff --git a/applications/utilities/mesh/generation/CV3DMesher/querySurface.C b/applications/utilities/mesh/generation/CV3DMesher/querySurface.C
index 311c6350c11..4575abc1cbc 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/querySurface.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/querySurface.C
@@ -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);
diff --git a/applications/utilities/mesh/generation/CV3DMesher/querySurface.H b/applications/utilities/mesh/generation/CV3DMesher/querySurface.H
index 8e4a945c407..44e89380bf8 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/querySurface.H
+++ b/applications/utilities/mesh/generation/CV3DMesher/querySurface.H
@@ -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,
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/G b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/G
deleted file mode 100644
index cb9783649fe..00000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/G
+++ /dev/null
@@ -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;
- }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/T b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/T
deleted file mode 100644
index 19ee7d9d08e..00000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/T
+++ /dev/null
@@ -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;
- }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/U b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/U
deleted file mode 100644
index 9a0eaf66b54..00000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/U
+++ /dev/null
@@ -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);
- }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/epsilon b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/epsilon
deleted file mode 100644
index 5f8de5792ef..00000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/epsilon
+++ /dev/null
@@ -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;
- }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/k b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/k
deleted file mode 100644
index a7d2d3f4f07..00000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/k
+++ /dev/null
@@ -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;
- }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p
deleted file mode 100644
index f200e3eeb0f..00000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p
+++ /dev/null
@@ -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;
- }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/pd b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/pd
deleted file mode 100644
index 1841d7882f5..00000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/pd
+++ /dev/null
@@ -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;
- }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/constant/polyMesh/boundary b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/constant/polyMesh/boundary
index add6399d3f9..aebfe3f9a72 100644
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/constant/polyMesh/boundary
+++ b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/constant/polyMesh/boundary
@@ -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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
--
GitLab