Adding locally deduced alignment directions to each vertex to supply the target alignment to the cell/face relaxation method.

This commit is contained in:
graham
2008-11-26 20:14:24 +00:00
parent efd6f447a2
commit 93413854ed
6 changed files with 251 additions and 45 deletions

View File

@ -65,6 +65,198 @@ void Foam::CV3D::reinsertPoints(const pointField& points)
}
void Foam::CV3D::setVertexAlignmentDirections()
{
for
(
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
{
if (vit->internalOrBoundaryPoint())
{
List<vector> alignmentDirections(vit->alignmentDirections());
point vert(topoint(vit->point()));
pointIndexHit pHit = qSurf_.tree().findNearest
(
vert,
controls_.nearWallAlignedDist2
);
if (pHit.hit())
{
// Primary alignment
vector np = qSurf_.faceNormals()[pHit.index()];
// Generate equally spaced 'spokes' in a circle normal to the
// direction from the vertex to the closest point on the surface
// and look for a secondary intersection.
vector d = pHit.hitPoint() - vert;
tensor R = rotationTensor(vector(0,0,1), np);
label s = 36;
scalar closestSpokeHitDistance = GREAT;
point closestSpokeHitPoint = point(GREAT,GREAT,GREAT);
label closestSpokeHitDistanceIndex = -1;
for(label i = 0; i < s; i++)
{
vector spoke
(
Foam::cos(i*mathematicalConstant::twoPi/s),
Foam::sin(i*mathematicalConstant::twoPi/s),
0
);
spoke *= controls_.nearWallAlignedDist;
spoke = R & spoke;
pointIndexHit spokeHit;
// internal spoke
spokeHit = qSurf_.tree().findLine
(
vert,
vert + spoke
);
if (spokeHit.hit())
{
scalar spokeHitDistance = mag
(
spokeHit.hitPoint() - vert
);
if (spokeHitDistance < closestSpokeHitDistance)
{
closestSpokeHitPoint = spokeHit.hitPoint();
closestSpokeHitDistanceIndex = spokeHit.index();
}
}
//external spoke
point mirrorVert = vert + 2*d;
spokeHit = qSurf_.tree().findLine
(
mirrorVert,
mirrorVert + spoke
);
if (spokeHit.hit())
{
scalar spokeHitDistance = mag
(
spokeHit.hitPoint() - mirrorVert
);
if (spokeHitDistance < closestSpokeHitDistance)
{
closestSpokeHitPoint = spokeHit.hitPoint();
closestSpokeHitDistanceIndex = spokeHit.index();
}
}
}
if (closestSpokeHitDistanceIndex > -1)
{
// Auxiliary alignment generated by spoke intersection
// normal.
vector na =
qSurf_.faceNormals()[closestSpokeHitDistanceIndex];
// Secondary alignment
vector ns = np ^ na;
if (mag(ns) < SMALL)
{
FatalErrorIn("Foam::CV3D::setVertexAlignmentDirections")
<< "Parallel normals detected in spoke search."
<< nl << exit(FatalError);
}
ns /= mag(ns);
// Tertiary alignment
vector nt = ns ^ np;
alignmentDirections.setSize(3);
alignmentDirections[0] = np;
alignmentDirections[1] = ns;
alignmentDirections[2] = nt;
Info<< "internal " << vit->internalPoint()
<< nl << alignmentDirections
<< nl << "v " << vert + alignmentDirections[0]
<< nl << "v " << vert + alignmentDirections[1]
<< nl << "v " << vert + alignmentDirections[2]
<< nl << "v " << vert
<< nl << "v " << pHit.hitPoint()
<< nl << "v " << closestSpokeHitPoint
<< nl << "f 4 1"
<< nl << "f 4 2"
<< nl << "f 4 3"
<< nl << endl;
}
else
{
// Using on the primary alignment and the closest of the
// globalAlignmentDirections
alignmentDirections.setSize(0);
}
}
else
{
alignmentDirections.setSize(0);
}
}
}
}
Foam::scalar Foam::CV3D::alignmentDistanceWeight(scalar dist) const
{
scalar w;
scalar x = dist/controls_.nearWallAlignedDist;
if (x < 0.5)
{
w = 1;
}
else if (x < 1)
{
w = 2*(1 - x);
}
else
{
w = 0;
}
return w;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::CV3D::CV3D
@ -208,36 +400,6 @@ void Foam::CV3D::insertGrid()
Random rndGen(1321);
scalar pert = controls_.randomPerturbation*cmptMin(delta);
// for (int i=0; i<ni; i++)
// {
// for (int j=0; j<nj; j++)
// {
// for (int k=0; k<nk; k++)
// {
// point p
// (
// x0 + i*delta.x(),
// y0 + j*delta.y(),
// z0 + k*delta.z()
// );
// if (controls_.randomiseInitialGrid)
// {
// p.x() += pert*(rndGen.scalar01() - 0.5);
// p.y() += pert*(rndGen.scalar01() - 0.5);
// p.z() += pert*(rndGen.scalar01() - 0.5);
// }
// if (qSurf_.wellInside(p, 0.5*controls_.minCellSize2))
// {
// insert(Point(p.x(), p.y(), p.z()))->index() = nVert++;
// }
// }
// }
// }
// Info<< nVert - startOfInternalPoints_ << " vertices inserted" << endl;
std::vector<Point> initialPoints;
for (int i=0; i<ni; i++)
@ -263,6 +425,7 @@ void Foam::CV3D::insertGrid()
if (qSurf_.wellInside(p, 0.5*controls_.minCellSize2))
{
initialPoints.push_back(Point(p.x(), p.y(), p.z()));
// insert(Point(p.x(), p.y(), p.z()))->index() = nVert++;
}
}
}
@ -341,12 +504,10 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
dualVerti++;
}
// else
// {
// cit->cellIndex() = -1;
// }
}
setVertexAlignmentDirections();
dualVertices.setSize(dualVerti);
// loop around the Delaunay edges to construct the dual faces.

View File

@ -113,6 +113,9 @@ public:
//- Square of nearWallAlignedDist
scalar nearWallAlignedDist2;
//- Alignment direction vectors
vectorField globalAlignmentDirections;
//- Chose if the cell orientation should relax during the iterations
// or remain fixed to the x-y directions
Switch relaxOrientation;
@ -331,6 +334,10 @@ private:
labelList& patchStarts
);
void setVertexAlignmentDirections();
scalar alignmentDistanceWeight(scalar dist) const;
public:

View File

@ -166,6 +166,7 @@ void Foam::CV3D::writeMesh(bool writeToTimestep)
patchStarts
);
writeDual(points, faces, "dualMesh.obj");
IOobject io
(
@ -198,11 +199,19 @@ void Foam::CV3D::writeMesh(bool writeToTimestep)
polyMesh pMesh
(
io,
points,
faces,
owner,
neighbour
xferMove(points),
xferMove(faces),
xferMove(owner),
xferMove(neighbour)
);
// polyMesh pMesh
// (
// io,
// points,
// faces,
// owner,
// neighbour
// );
List<polyPatch*> patches(patchStarts.size());
@ -227,7 +236,7 @@ void Foam::CV3D::writeMesh(bool writeToTimestep)
<< exit(FatalError);
}
writeDual(points, faces, "dualMesh.obj");
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //

View File

@ -73,7 +73,7 @@ int main(int argc, char *argv[])
<< "Read surface with " << surf.size() << " triangles from file "
<< args.args()[1] << nl << endl;
// Read and triangulation
// Read and triangulation
// ~~~~~~~~~~~~~~~~~~~~~~
CV3D mesh(runTime, surf);

View File

@ -45,6 +45,7 @@ Foam::CV3D::controls::controls(const IOdictionary& controlDict)
readScalar(controlDict.lookup("nearWallAlignedDist"))*minCellSize
),
nearWallAlignedDist2(Foam::sqr(nearWallAlignedDist)),
globalAlignmentDirections(controlDict.lookup("globalAlignmentDirections")),
relaxOrientation(controlDict.lookup("relaxOrientation")),
insertSurfaceNearestPointPairs
(
@ -62,7 +63,17 @@ Foam::CV3D::controls::controls(const IOdictionary& controlDict)
writeFinalTriangulation(controlDict.lookup("writeFinalTriangulation")),
randomiseInitialGrid(controlDict.lookup("randomiseInitialGrid")),
randomPerturbation(readScalar(controlDict.lookup("randomPerturbation")))
{}
{
if(globalAlignmentDirections.size() < 3)
{
FatalErrorIn("Foam::CV3D::controls::controls")
<< "Too many alignment directions specified - 3 at maximum."
<< exit(FatalError);
}
globalAlignmentDirections /= mag(globalAlignmentDirections);
}
// ************************************************************************* //

View File

@ -35,6 +35,8 @@ Description
#define indexedVertex_H
#include <CGAL/Triangulation_3.h>
#include "List.H"
#include "vector.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -63,6 +65,8 @@ class indexedVertex
// Lowest numbered is inside one (master).
int type_;
Foam::List<Foam::vector> alignmentDirections_;
public:
@ -89,28 +93,32 @@ public:
:
Vb(),
index_(INTERNAL_POINT),
type_(INTERNAL_POINT)
type_(INTERNAL_POINT),
alignmentDirections_(0)
{}
indexedVertex(const Point& p)
:
Vb(p),
index_(INTERNAL_POINT),
type_(INTERNAL_POINT)
type_(INTERNAL_POINT),
alignmentDirections_(0)
{}
indexedVertex(const Point& p, Cell_handle f)
:
Vb(f, p),
index_(INTERNAL_POINT),
type_(INTERNAL_POINT)
type_(INTERNAL_POINT),
alignmentDirections_(0)
{}
indexedVertex(Cell_handle f)
:
Vb(f),
index_(INTERNAL_POINT),
type_(INTERNAL_POINT)
type_(INTERNAL_POINT),
alignmentDirections_(0)
{}
@ -135,6 +143,16 @@ public:
return type_;
}
Foam::List<Foam::vector>& alignmentDirections()
{
return alignmentDirections_;
}
const Foam::List<Foam::vector>& alignmentDirections() const
{
return alignmentDirections_;
}
inline bool uninitialised() const
{
return type_ == INTERNAL_POINT && index_ == INTERNAL_POINT;