mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
Range insertion witn indexing for initial grid. Minor Info line tweaks and superfluous file removal.
This commit is contained in:
graham
@ -58,7 +58,7 @@ void Foam::CV3D::reinsertPoints(const pointField& points)
|
||||
{
|
||||
const point& p = points[i];
|
||||
|
||||
insert(toPoint(p))->index() = nVert++;
|
||||
insert(toPoint(p))->index() = nVert++;
|
||||
}
|
||||
|
||||
Info<< nVert << " vertices reinserted" << endl;
|
||||
@ -158,7 +158,6 @@ void Foam::CV3D::insertPoints
|
||||
assert(is_valid());
|
||||
|
||||
writeTriangles("initial_triangles.obj", true);
|
||||
// writeFaces("initial_faces.obj", true);
|
||||
}
|
||||
}
|
||||
|
||||
@ -183,11 +182,13 @@ void Foam::CV3D::insertPoints(const fileName& pointFileName)
|
||||
|
||||
void Foam::CV3D::insertGrid()
|
||||
{
|
||||
Info<< "insertInitialGrid: ";
|
||||
Info<< nl << "Inserting initial grid." << endl;
|
||||
|
||||
startOfInternalPoints_ = number_of_vertices();
|
||||
label nVert = startOfInternalPoints_;
|
||||
|
||||
Info<< nl << nVert << " existing vertices." << endl;
|
||||
|
||||
scalar x0 = qSurf_.bb().min().x();
|
||||
scalar xR = qSurf_.bb().max().x() - x0;
|
||||
int ni = int(xR/controls_.minCellSize) + 1;
|
||||
@ -207,6 +208,38 @@ 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++)
|
||||
{
|
||||
for (int j=0; j<nj; j++)
|
||||
@ -229,62 +262,32 @@ void Foam::CV3D::insertGrid()
|
||||
|
||||
if (qSurf_.wellInside(p, 0.5*controls_.minCellSize2))
|
||||
{
|
||||
insert(Point(p.x(), p.y(), p.z()))->index() = nVert++;
|
||||
initialPoints.push_back(Point(p.x(), p.y(), p.z()));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Info<< nVert << " vertices inserted" << nl << endl;
|
||||
Info<< nl << initialPoints.size() << " vertices to insert." << endl;
|
||||
|
||||
// std::vector<Point> initialPoints;
|
||||
// using the range insert (it is faster than inserting points one by one)
|
||||
insert(initialPoints.begin(), initialPoints.end());
|
||||
|
||||
// for (int i=0; i<ni; i++)
|
||||
// {
|
||||
// for (int j=0; j<nj; j++)
|
||||
// {
|
||||
// for (int k=0; k<nk; k++)
|
||||
// {
|
||||
// point p1
|
||||
// (
|
||||
// x0 + i*delta.x(),
|
||||
// y0 + j*delta.y(),
|
||||
// z0 + k*delta.z()
|
||||
// );
|
||||
Info<< nl << number_of_vertices() - startOfInternalPoints_
|
||||
<< " vertices inserted." << endl;
|
||||
|
||||
// point p2 = p1 + 0.5*delta;
|
||||
|
||||
// if (controls_.randomiseInitialGrid)
|
||||
// {
|
||||
// p1.x() += pert*(rndGen.scalar01() - 0.5);
|
||||
// p1.y() += pert*(rndGen.scalar01() - 0.5);
|
||||
// p1.z() += pert*(rndGen.scalar01() - 0.5);
|
||||
// }
|
||||
|
||||
// if (qSurf_.wellInside(p1, 0.5*controls_.minCellSize2))
|
||||
// {
|
||||
// initialPoints.push_back(Point(p1.x(), p1.y(), p1.z()));
|
||||
// }
|
||||
|
||||
// if (controls_.randomiseInitialGrid)
|
||||
// {
|
||||
// p2.x() += pert*(rndGen.scalar01() - 0.5);
|
||||
// p2.y() += pert*(rndGen.scalar01() - 0.5);
|
||||
// p2.z() += pert*(rndGen.scalar01() - 0.5);
|
||||
// }
|
||||
|
||||
// if (qSurf_.wellInside(p2, 0.5*controls_.minCellSize2))
|
||||
// {
|
||||
// initialPoints.push_back(Point(p2.x(), p2.y(), p2.z()));
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
|
||||
// Info<< initialPoints.size() << " vertices inserted" << nl << endl;
|
||||
|
||||
// // using the range insert (it is faster than inserting points one by one)
|
||||
// insert(initialPoints.begin(), initialPoints.end());
|
||||
for
|
||||
(
|
||||
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
|
||||
vit != finite_vertices_end();
|
||||
++vit
|
||||
)
|
||||
{
|
||||
if (vit->uninitialised())
|
||||
{
|
||||
vit->index() = nVert++;
|
||||
}
|
||||
}
|
||||
|
||||
if (controls_.writeInitialTriangulation)
|
||||
{
|
||||
|
||||
@ -136,8 +136,8 @@ void Foam::CV3D::writeTriangles(const fileName& fName, bool internalOnly) const
|
||||
|
||||
if(writeFacet)
|
||||
{
|
||||
str << "f "
|
||||
<< vertexMap[c->vertex(facetIndices[0])->index()] + 1
|
||||
str << "f "
|
||||
<< vertexMap[c->vertex(facetIndices[0])->index()] + 1
|
||||
<< ' ' << vertexMap[c->vertex(facetIndices[1])->index()] + 1
|
||||
<< ' ' << vertexMap[c->vertex(facetIndices[2])->index()] + 1
|
||||
<< nl;
|
||||
|
||||
@ -1,8 +1,5 @@
|
||||
|
||||
#include CGAL_FILES
|
||||
|
||||
plane.C
|
||||
|
||||
querySurface.C
|
||||
CV3D.C
|
||||
controls.C
|
||||
|
||||
@ -1,416 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / 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"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
void Foam::CV3D::calcDualMesh
|
||||
(
|
||||
pointField& points,
|
||||
faceList& faces,
|
||||
labelList& owner,
|
||||
labelList& neighbour,
|
||||
wordList& patchNames,
|
||||
labelList& patchSizes,
|
||||
labelList& patchStarts
|
||||
)
|
||||
{
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ dual points ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
points.setSize(number_of_cells());
|
||||
|
||||
label dualVerti = 0;
|
||||
|
||||
for
|
||||
(
|
||||
Triangulation::Finite_cells_iterator cit = finite_cells_begin();
|
||||
cit != finite_cells_end();
|
||||
++cit
|
||||
)
|
||||
{
|
||||
if
|
||||
(
|
||||
cit->vertex(0)->internalOrBoundaryPoint()
|
||||
|| cit->vertex(1)->internalOrBoundaryPoint()
|
||||
|| cit->vertex(2)->internalOrBoundaryPoint()
|
||||
|| cit->vertex(3)->internalOrBoundaryPoint()
|
||||
)
|
||||
{
|
||||
cit->cellIndex() = dualVerti;
|
||||
points[dualVerti] = topoint(dual(cit));
|
||||
dualVerti++;
|
||||
}
|
||||
else
|
||||
{
|
||||
cit->cellIndex() = -1;
|
||||
}
|
||||
}
|
||||
|
||||
points.setSize(dualVerti);
|
||||
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~ dual cell indexing ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
// resets type and index information for Delaunay vertices
|
||||
// assigns an index to the vertices which will be the dual cell index used
|
||||
// for owner neighbour assignment
|
||||
|
||||
label dualCelli = 0;
|
||||
|
||||
for
|
||||
(
|
||||
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
|
||||
vit != finite_vertices_end();
|
||||
++vit
|
||||
)
|
||||
{
|
||||
if (vit->internalOrBoundaryPoint())
|
||||
{
|
||||
vit->type() = Vb::INTERNAL_POINT;
|
||||
vit->index() = dualCelli;
|
||||
dualCelli++;
|
||||
}
|
||||
else
|
||||
{
|
||||
vit->type() = Vb::FAR_POINT;
|
||||
vit->index() = -1;
|
||||
}
|
||||
}
|
||||
|
||||
// ~~~~~~~~~~~~ dual face and owner neighbour construction ~~~~~~~~~~~~~~~~~
|
||||
|
||||
label nPatches = 1;
|
||||
|
||||
patchNames.setSize(nPatches);
|
||||
|
||||
patchNames[0] = "CV3D_default_patch";
|
||||
|
||||
patchSizes.setSize(nPatches);
|
||||
|
||||
patchStarts.setSize(nPatches);
|
||||
|
||||
List<DynamicList<face> > patchFaces(nPatches, DynamicList<face>(0));
|
||||
|
||||
List<DynamicList<label> > patchOwners(nPatches, DynamicList<label>(0));
|
||||
|
||||
faces.setSize(number_of_edges());
|
||||
|
||||
owner.setSize(number_of_edges());
|
||||
|
||||
neighbour.setSize(number_of_edges());
|
||||
|
||||
label dualFacei = 0;
|
||||
|
||||
for
|
||||
(
|
||||
Triangulation::Finite_edges_iterator eit = finite_edges_begin();
|
||||
eit != finite_edges_end();
|
||||
++eit
|
||||
)
|
||||
{
|
||||
if
|
||||
(
|
||||
eit->first->vertex(eit->second)->internalOrBoundaryPoint()
|
||||
|| eit->first->vertex(eit->third)->internalOrBoundaryPoint()
|
||||
)
|
||||
{
|
||||
Cell_circulator ccStart = incident_cells(*eit);
|
||||
Cell_circulator cc = ccStart;
|
||||
|
||||
DynamicList<label> verticesOnFace;
|
||||
|
||||
do
|
||||
{
|
||||
if (!is_infinite(cc))
|
||||
{
|
||||
if (cc->cellIndex() < 0)
|
||||
{
|
||||
FatalErrorIn("Foam::CV3D::calcDualMesh")
|
||||
<< "Dual face uses circumcenter defined by a "
|
||||
<< " Delaunay tetrahedron with no internal "
|
||||
<< "or boundary points."
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
verticesOnFace.append(cc->cellIndex());
|
||||
}
|
||||
} while (++cc != ccStart);
|
||||
|
||||
verticesOnFace.shrink();
|
||||
|
||||
face newDualFace(verticesOnFace);
|
||||
|
||||
Cell_handle c = eit->first;
|
||||
Vertex_handle vA = c->vertex(eit->second);
|
||||
Vertex_handle vB = c->vertex(eit->third);
|
||||
|
||||
label dcA = vA->index();
|
||||
|
||||
if (!vA->internalOrBoundaryPoint())
|
||||
{
|
||||
dcA = -1;
|
||||
}
|
||||
|
||||
label dcB = vB->index();
|
||||
|
||||
if (!vB->internalOrBoundaryPoint())
|
||||
{
|
||||
dcB = -1;
|
||||
}
|
||||
|
||||
label dcOwn = -1;
|
||||
label dcNei = -1;
|
||||
|
||||
if (dcA == -1 && dcB == -1)
|
||||
{
|
||||
FatalErrorIn("calcDualMesh")
|
||||
<< "Attempting to create a face joining "
|
||||
<< "two external dual cells "
|
||||
<< exit(FatalError);
|
||||
}
|
||||
else if (dcA == -1 || dcB == -1)
|
||||
{
|
||||
// boundary face, find which is the owner
|
||||
|
||||
if (dcA == -1)
|
||||
{
|
||||
dcOwn = dcB;
|
||||
|
||||
// reverse face order to correctly orientate normal
|
||||
reverse(newDualFace);
|
||||
}
|
||||
else
|
||||
{
|
||||
dcOwn = dcA;
|
||||
}
|
||||
|
||||
// find which patch this face is on; Hardcoded for now.
|
||||
label patchIndex = 0;
|
||||
|
||||
patchFaces[patchIndex].append(newDualFace);
|
||||
patchOwners[patchIndex].append(dcOwn);
|
||||
}
|
||||
else
|
||||
{
|
||||
// internal face, find the lower cell to be the owner
|
||||
|
||||
if (dcB > dcA)
|
||||
{
|
||||
dcOwn = dcA;
|
||||
dcNei = dcB;
|
||||
}
|
||||
else
|
||||
{
|
||||
dcOwn = dcB;
|
||||
dcNei = dcA;
|
||||
|
||||
// reverse face order to correctly orientate normal
|
||||
// unsure if CGAL always circulates consistently,
|
||||
// needs to be more rigorous
|
||||
reverse(newDualFace);
|
||||
}
|
||||
|
||||
faces[dualFacei] = newDualFace;
|
||||
|
||||
owner[dualFacei] = dcOwn;
|
||||
|
||||
neighbour[dualFacei] = dcNei;
|
||||
|
||||
dualFacei++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
label nInternalFaces = dualFacei;
|
||||
|
||||
faces.setSize(nInternalFaces);
|
||||
|
||||
owner.setSize(nInternalFaces);
|
||||
|
||||
neighbour.setSize(nInternalFaces);
|
||||
|
||||
// ~~~~~~~~ sort owner, reordinging neighbour and faces to match ~~~~~~~~~~~
|
||||
// two stage sort for upper triangular order: sort by owner first, then for
|
||||
// each block of owners sort by neighbour
|
||||
|
||||
labelList sortingIndices;
|
||||
|
||||
// Stage 1
|
||||
|
||||
{
|
||||
SortableList<label> sortedOwner(owner);
|
||||
|
||||
sortingIndices = sortedOwner.indices();
|
||||
}
|
||||
|
||||
{
|
||||
labelList copyOwner(owner.size());
|
||||
|
||||
forAll(sortingIndices, sI)
|
||||
{
|
||||
copyOwner[sI] = owner[sortingIndices[sI]];
|
||||
}
|
||||
|
||||
owner = copyOwner;
|
||||
}
|
||||
|
||||
{
|
||||
labelList copyNeighbour(neighbour.size());
|
||||
|
||||
forAll(sortingIndices, sI)
|
||||
{
|
||||
copyNeighbour[sI] = neighbour[sortingIndices[sI]];
|
||||
}
|
||||
|
||||
neighbour = copyNeighbour;
|
||||
}
|
||||
|
||||
{
|
||||
faceList copyFaces(faces.size());
|
||||
|
||||
forAll(sortingIndices, sI)
|
||||
{
|
||||
copyFaces[sI] = faces[sortingIndices[sI]];
|
||||
}
|
||||
|
||||
faces = copyFaces;
|
||||
}
|
||||
|
||||
// Stage 2
|
||||
|
||||
sortingIndices = -1;
|
||||
|
||||
DynamicList<label> ownerCellJumps;
|
||||
|
||||
// Force first owner entry to be a jump
|
||||
ownerCellJumps.append(0);
|
||||
|
||||
for (label o = 1; o < owner.size(); o++)
|
||||
{
|
||||
if (owner[o] > owner[o-1])
|
||||
{
|
||||
ownerCellJumps.append(o);
|
||||
}
|
||||
}
|
||||
|
||||
ownerCellJumps.shrink();
|
||||
|
||||
forAll(ownerCellJumps, oCJ)
|
||||
{
|
||||
label start = ownerCellJumps[oCJ];
|
||||
|
||||
label length;
|
||||
|
||||
if (oCJ == ownerCellJumps.size() - 1)
|
||||
{
|
||||
length = owner.size() - start;
|
||||
}
|
||||
else
|
||||
{
|
||||
length = ownerCellJumps[oCJ + 1] - start;
|
||||
}
|
||||
|
||||
SubList<label> neighbourBlock(neighbour, length, start);
|
||||
|
||||
SortableList<label> sortedNeighbourBlock(neighbourBlock);
|
||||
|
||||
forAll(sortedNeighbourBlock, sNB)
|
||||
{
|
||||
sortingIndices[start + sNB] =
|
||||
sortedNeighbourBlock.indices()[sNB] + start;
|
||||
}
|
||||
}
|
||||
|
||||
// Perform sort
|
||||
|
||||
{
|
||||
labelList copyOwner(owner.size());
|
||||
|
||||
forAll(sortingIndices, sI)
|
||||
{
|
||||
copyOwner[sI] = owner[sortingIndices[sI]];
|
||||
}
|
||||
|
||||
owner = copyOwner;
|
||||
}
|
||||
|
||||
{
|
||||
labelList copyNeighbour(neighbour.size());
|
||||
|
||||
forAll(sortingIndices, sI)
|
||||
{
|
||||
copyNeighbour[sI] = neighbour[sortingIndices[sI]];
|
||||
}
|
||||
|
||||
neighbour = copyNeighbour;
|
||||
}
|
||||
|
||||
{
|
||||
faceList copyFaces(faces.size());
|
||||
|
||||
forAll(sortingIndices, sI)
|
||||
{
|
||||
copyFaces[sI] = faces[sortingIndices[sI]];
|
||||
}
|
||||
|
||||
faces = copyFaces;
|
||||
}
|
||||
|
||||
// ~~~~~~~~ add patch information ~~~~~~~~~~~
|
||||
|
||||
label nBoundaryFaces = 0;
|
||||
|
||||
forAll(patchFaces, p)
|
||||
{
|
||||
patchFaces[p].shrink();
|
||||
|
||||
patchOwners[p].shrink();
|
||||
|
||||
patchSizes[p] = patchFaces[p].size();
|
||||
|
||||
patchStarts[p] = nInternalFaces + nBoundaryFaces;
|
||||
|
||||
nBoundaryFaces += patchSizes[p];
|
||||
}
|
||||
|
||||
faces.setSize(nInternalFaces + nBoundaryFaces);
|
||||
|
||||
owner.setSize(nInternalFaces + nBoundaryFaces);
|
||||
|
||||
forAll(patchFaces, p)
|
||||
{
|
||||
forAll(patchFaces[p], f)
|
||||
{
|
||||
faces[dualFacei] = patchFaces[p][f];
|
||||
|
||||
owner[dualFacei] = patchOwners[p][f];
|
||||
|
||||
dualFacei++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -135,6 +135,10 @@ public:
|
||||
return type_;
|
||||
}
|
||||
|
||||
inline bool uninitialised() const
|
||||
{
|
||||
return type_ == INTERNAL_POINT && index_ == INTERNAL_POINT;
|
||||
}
|
||||
|
||||
//- Is point a far-point
|
||||
inline bool farPoint() const
|
||||
|
||||
@ -614,7 +614,7 @@ void Foam::CV3D::reinsertFeaturePoints()
|
||||
if (nVert == number_of_vertices())
|
||||
{
|
||||
FatalErrorIn("Foam::CV3D::reinsertFeaturePoints")
|
||||
<< "Failed to reinsert feature point " << topoint(fPt)
|
||||
<< "Failed to reinsert feature point " << topoint(fPt)
|
||||
<< endl;
|
||||
}
|
||||
|
||||
|
||||
@ -814,7 +814,7 @@ void Foam::CV3D::insertEdgePointGroups
|
||||
|
||||
void Foam::CV3D::insertSurfaceNearestPointPairs()
|
||||
{
|
||||
Info<< "insertSurfaceNearestPointPairs: " << nl << endl;
|
||||
Info<< nl << "insertSurfaceNearestPointPairs: " << endl;
|
||||
|
||||
label nSurfacePointsEst = number_of_vertices();
|
||||
|
||||
|
||||
@ -1,481 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 1991-2008 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 "plane.H"
|
||||
#include "tensor.H"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
// Calculate base point and unit normal vector from plane equation
|
||||
void Foam::plane::calcPntAndVec(const scalarList& C)
|
||||
{
|
||||
if (mag(C[0]) > VSMALL)
|
||||
{
|
||||
basePoint_ = vector((-C[3]/C[0]), 0, 0);
|
||||
}
|
||||
else
|
||||
{
|
||||
if (mag(C[1]) > VSMALL)
|
||||
{
|
||||
basePoint_ = vector(0, (-C[3]/C[1]), 0);
|
||||
}
|
||||
else
|
||||
{
|
||||
if (mag(C[2]) > VSMALL)
|
||||
{
|
||||
basePoint_ = vector(0, 0, (-C[3]/C[2]));
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn("void plane::calcPntAndVec(const scalarList&)")
|
||||
<< "At least one plane coefficient must have a value"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
unitVector_ = vector(C[0], C[1], C[2]);
|
||||
scalar magUnitVector(mag(unitVector_));
|
||||
|
||||
if (magUnitVector < VSMALL)
|
||||
{
|
||||
FatalErrorIn("void plane::calcPntAndVec(const scalarList&)")
|
||||
<< "Plane normal defined with zero length"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
unitVector_ /= magUnitVector;
|
||||
}
|
||||
|
||||
|
||||
void Foam::plane::calcPntAndVec
|
||||
(
|
||||
const point& point1,
|
||||
const point& point2,
|
||||
const point& point3
|
||||
)
|
||||
{
|
||||
basePoint_ = (point1 + point2 + point3)/3;
|
||||
vector line12 = point1 - point2;
|
||||
vector line23 = point2 - point3;
|
||||
|
||||
if
|
||||
(
|
||||
mag(line12) < VSMALL
|
||||
|| mag(line23) < VSMALL
|
||||
|| mag(point3-point1) < VSMALL
|
||||
)
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"void plane::calcPntAndVec\n"
|
||||
"(\n"
|
||||
" const point&,\n"
|
||||
" const point&,\n"
|
||||
" const point&\n"
|
||||
")\n"
|
||||
) << "Bad points." << abort(FatalError);
|
||||
}
|
||||
|
||||
unitVector_ = line12 ^ line23;
|
||||
scalar magUnitVector(mag(unitVector_));
|
||||
|
||||
if (magUnitVector < VSMALL)
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"void plane::calcPntAndVec\n"
|
||||
"(\n"
|
||||
" const point&,\n"
|
||||
" const point&,\n"
|
||||
" const point&\n"
|
||||
")\n"
|
||||
) << "Plane normal defined with zero length"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
unitVector_ /= magUnitVector;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
// Construct from normal vector through the origin
|
||||
Foam::plane::plane(const vector& normalVector)
|
||||
:
|
||||
unitVector_(normalVector),
|
||||
basePoint_(vector::zero)
|
||||
{}
|
||||
|
||||
|
||||
// Construct from point and normal vector
|
||||
Foam::plane::plane(const point& basePoint, const vector& normalVector)
|
||||
:
|
||||
unitVector_(normalVector),
|
||||
basePoint_(basePoint)
|
||||
{
|
||||
scalar magUnitVector(mag(unitVector_));
|
||||
|
||||
if (magUnitVector > VSMALL)
|
||||
{
|
||||
unitVector_ /= magUnitVector;
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn("plane::plane(const point&, const vector&)")
|
||||
<< "plane normal has got zero length"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Construct from plane equation
|
||||
Foam::plane::plane(const scalarList& C)
|
||||
{
|
||||
calcPntAndVec(C);
|
||||
}
|
||||
|
||||
|
||||
// Construct from three points
|
||||
Foam::plane::plane
|
||||
(
|
||||
const point& a,
|
||||
const point& b,
|
||||
const point& c
|
||||
)
|
||||
{
|
||||
calcPntAndVec(a, b, c);
|
||||
}
|
||||
|
||||
|
||||
// Construct from dictionary
|
||||
Foam::plane::plane(const dictionary& dict)
|
||||
:
|
||||
unitVector_(vector::zero),
|
||||
basePoint_(point::zero)
|
||||
{
|
||||
word planeType(dict.lookup("planeType"));
|
||||
|
||||
if (planeType == "planeEquation")
|
||||
{
|
||||
const dictionary& subDict = dict.subDict("planeEquationDict");
|
||||
scalarList C(4);
|
||||
|
||||
C[0] = readScalar(subDict.lookup("a"));
|
||||
C[1] = readScalar(subDict.lookup("b"));
|
||||
C[2] = readScalar(subDict.lookup("c"));
|
||||
C[3] = readScalar(subDict.lookup("d"));
|
||||
|
||||
calcPntAndVec(C);
|
||||
|
||||
}
|
||||
else if (planeType == "embeddedPoints")
|
||||
{
|
||||
const dictionary& subDict = dict.subDict("embeddedPoints");
|
||||
|
||||
point point1(subDict.lookup("point1"));
|
||||
point point2(subDict.lookup("point2"));
|
||||
point point3(subDict.lookup("point3"));
|
||||
|
||||
calcPntAndVec(point1, point2, point3);
|
||||
}
|
||||
else if (planeType == "pointAndNormal")
|
||||
{
|
||||
const dictionary& subDict = dict.subDict("pointAndNormalDict");
|
||||
|
||||
basePoint_ = subDict.lookup("basePoint");
|
||||
unitVector_ = subDict.lookup("normalVector");
|
||||
unitVector_ /= mag(unitVector_);
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalIOErrorIn
|
||||
(
|
||||
"plane::plane(const dictionary&)",
|
||||
dict
|
||||
)
|
||||
<< "Invalid plane type: " << planeType
|
||||
<< abort(FatalIOError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Construct from Istream. Assumes point and normal vector.
|
||||
Foam::plane::plane(Istream& is)
|
||||
:
|
||||
unitVector_(is),
|
||||
basePoint_(is)
|
||||
{
|
||||
scalar magUnitVector(mag(unitVector_));
|
||||
|
||||
if (magUnitVector > VSMALL)
|
||||
{
|
||||
unitVector_ /= magUnitVector;
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn("plane::plane(Istream& is)")
|
||||
<< "plane normal has got zero length"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
// Return plane normal vector
|
||||
const Foam::vector& Foam::plane::normal() const
|
||||
{
|
||||
return unitVector_;
|
||||
}
|
||||
|
||||
|
||||
// Return plane base point
|
||||
const Foam::point& Foam::plane::refPoint() const
|
||||
{
|
||||
return basePoint_;
|
||||
}
|
||||
|
||||
|
||||
// Return coefficcients for plane equation: ax + by + cz + d = 0
|
||||
Foam::scalarList Foam::plane::planeCoeffs() const
|
||||
{
|
||||
scalarList C(4);
|
||||
|
||||
scalar magX = mag(unitVector_.x());
|
||||
scalar magY = mag(unitVector_.y());
|
||||
scalar magZ = mag(unitVector_.z());
|
||||
|
||||
if (magX > magY)
|
||||
{
|
||||
if (magX > magZ)
|
||||
{
|
||||
C[0] = 1;
|
||||
C[1] = unitVector_.y()/unitVector_.x();
|
||||
C[2] = unitVector_.z()/unitVector_.x();
|
||||
}
|
||||
else
|
||||
{
|
||||
C[0] = 0;
|
||||
C[1] = 0;
|
||||
C[2] = 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (magY > magZ)
|
||||
{
|
||||
C[0] = 0;
|
||||
C[1] = 1;
|
||||
C[2] = unitVector_.z()/unitVector_.y();
|
||||
}
|
||||
else
|
||||
{
|
||||
C[0] = 0;
|
||||
C[1] = 0;
|
||||
C[2] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
C[3] = - C[0] * basePoint_.x()
|
||||
- C[1] * basePoint_.y()
|
||||
- C[2] * basePoint_.z();
|
||||
|
||||
return C;
|
||||
}
|
||||
|
||||
|
||||
// Return nearest point in the plane for the given point
|
||||
Foam::point Foam::plane::nearestPoint(const point& p) const
|
||||
{
|
||||
return p - unitVector_*((p - basePoint_) & unitVector_);
|
||||
}
|
||||
|
||||
|
||||
// Return distance from the given point to the plane
|
||||
Foam::scalar Foam::plane::distance(const point& p) const
|
||||
{
|
||||
return mag((p - basePoint_) & unitVector_);
|
||||
}
|
||||
|
||||
|
||||
// Cutting point for plane and line defined by origin and direction
|
||||
Foam::scalar Foam::plane::normalIntersect
|
||||
(
|
||||
const point& pnt0,
|
||||
const vector& dir
|
||||
) const
|
||||
{
|
||||
scalar denom = stabilise((dir & unitVector_), VSMALL);
|
||||
|
||||
return ((basePoint_ - pnt0) & unitVector_)/denom;
|
||||
}
|
||||
|
||||
|
||||
// Cutting line of two planes
|
||||
Foam::plane::ray Foam::plane::planeIntersect(const plane& plane2) const
|
||||
{
|
||||
// Mathworld plane-plane intersection. Assume there is a point on the
|
||||
// intersection line with z=0 and solve the two plane equations
|
||||
// for that (now 2x2 equation in x and y)
|
||||
// Better: use either z=0 or x=0 or y=0.
|
||||
|
||||
const vector& n1 = normal();
|
||||
const vector& n2 = plane2.normal();
|
||||
|
||||
const point& p1 = refPoint();
|
||||
const point& p2 = plane2.refPoint();
|
||||
|
||||
scalar n1p1 = n1&p1;
|
||||
scalar n2p2 = n2&p2;
|
||||
|
||||
vector dir = n1 ^ n2;
|
||||
|
||||
// Determine zeroed out direction (can be x,y or z) by looking at which
|
||||
// has the largest component in dir.
|
||||
scalar magX = mag(dir.x());
|
||||
scalar magY = mag(dir.y());
|
||||
scalar magZ = mag(dir.z());
|
||||
|
||||
direction iZero, i1, i2;
|
||||
|
||||
if (magX > magY)
|
||||
{
|
||||
if (magX > magZ)
|
||||
{
|
||||
iZero = 0;
|
||||
i1 = 1;
|
||||
i2 = 2;
|
||||
}
|
||||
else
|
||||
{
|
||||
iZero = 2;
|
||||
i1 = 0;
|
||||
i2 = 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (magY > magZ)
|
||||
{
|
||||
iZero = 1;
|
||||
i1 = 2;
|
||||
i2 = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
iZero = 2;
|
||||
i1 = 0;
|
||||
i2 = 1;
|
||||
}
|
||||
}
|
||||
|
||||
vector pt;
|
||||
|
||||
pt[iZero] = 0;
|
||||
pt[i1] = (n2[i2]*n1p1 - n1[i2]*n2p2) / (n1[i1]*n2[i2] - n2[i1]*n1[i2]);
|
||||
pt[i2] = (n2[i1]*n1p1 - n1[i1]*n2p2) / (n1[i2]*n2[i1] - n1[i1]*n2[i2]);
|
||||
|
||||
return ray(pt, dir);
|
||||
}
|
||||
|
||||
|
||||
// Cutting point of three planes
|
||||
Foam::point Foam::plane::planePlaneIntersect
|
||||
(
|
||||
const plane& plane2,
|
||||
const plane& plane3
|
||||
) const
|
||||
{
|
||||
List<scalarList> pcs(3);
|
||||
pcs[0]= planeCoeffs();
|
||||
pcs[1]= plane2.planeCoeffs();
|
||||
pcs[2]= plane3.planeCoeffs();
|
||||
|
||||
tensor a
|
||||
(
|
||||
pcs[0][0],pcs[0][1],pcs[0][2],
|
||||
pcs[1][0],pcs[1][1],pcs[1][2],
|
||||
pcs[2][0],pcs[2][1],pcs[2][2]
|
||||
);
|
||||
|
||||
vector b(pcs[0][3],pcs[1][3],pcs[2][3]);
|
||||
|
||||
scalar d = det(a);
|
||||
|
||||
if (mag(d) > SMALL)
|
||||
{
|
||||
return (inv(a, d) & (-b));
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"plane::planePlaneIntersect"
|
||||
"(const plane&, const plane&) const"
|
||||
) << "Planes don't have an intersection point." << endl
|
||||
<< "plane1:" << *this << endl
|
||||
<< "plane2:" << plane2 << endl
|
||||
<< "plane3:" << plane3 << endl
|
||||
<< "determinant: " << d
|
||||
<< exit(FatalError);
|
||||
|
||||
return vector::zero;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
|
||||
|
||||
bool Foam::operator==(const plane& a, const plane& b)
|
||||
{
|
||||
if (a.basePoint_ == b.basePoint_ && a.unitVector_ == b.unitVector_)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
else
|
||||
{
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
bool Foam::operator!=(const plane& a, const plane& b)
|
||||
{
|
||||
return !(a == b);
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::Ostream& Foam::operator<<(Ostream& os, const plane& a)
|
||||
{
|
||||
os << a.unitVector_ << token::SPACE << a.basePoint_;
|
||||
|
||||
return os;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
Reference in New Issue
Block a user