mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
Changed mistake in declaring alignmentDirections, now a non-const reference. Added initial loop outline for cell based looping.
This commit is contained in:
graham
@ -76,7 +76,7 @@ void Foam::CV3D::setVertexAlignmentDirections()
|
||||
{
|
||||
if (vit->internalOrBoundaryPoint())
|
||||
{
|
||||
List<vector> alignmentDirections(vit->alignmentDirections());
|
||||
List<vector>& alignmentDirections(vit->alignmentDirections());
|
||||
|
||||
point vert(topoint(vit->point()));
|
||||
|
||||
@ -233,8 +233,6 @@ void Foam::CV3D::setVertexAlignmentDirections()
|
||||
{
|
||||
alignmentDirections.setSize(0);
|
||||
}
|
||||
|
||||
vit->alignmentDirections() = alignmentDirections;
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -500,7 +498,7 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
|
||||
|
||||
pointField displacementAccumulator(startOfSurfacePointPairs_, vector::zero);
|
||||
|
||||
//scalarField weightAccumulator(startOfSurfacePointPairs_,0);
|
||||
scalarField weightAccumulator(startOfSurfacePointPairs_, 0);
|
||||
|
||||
label dualVerti = 0;
|
||||
|
||||
@ -631,311 +629,19 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
|
||||
// Rotate faces that are sufficiently large and well enough aligned with the
|
||||
// cell alignment direction(s)
|
||||
|
||||
vector n2 = vector(1,1,1);
|
||||
// vector n2 = vector(1,1,1);
|
||||
|
||||
n2 /= mag(n2);
|
||||
// n2 /= mag(n2);
|
||||
|
||||
tensor R = rotationTensor(vector(1,0,0), n2);
|
||||
// tensor R = rotationTensor(vector(1,0,0), n2);
|
||||
|
||||
List<vector> globalAlignmentDirs(3);
|
||||
// List<vector> globalAlignmentDirs(3);
|
||||
|
||||
globalAlignmentDirs[0] = R & vector(1,0,0);
|
||||
globalAlignmentDirs[1] = R & vector(0,1,0);
|
||||
globalAlignmentDirs[2] = R & vector(0,0,1);
|
||||
// globalAlignmentDirs[0] = R & vector(1,0,0);
|
||||
// globalAlignmentDirs[1] = R & vector(0,1,0);
|
||||
// globalAlignmentDirs[2] = R & vector(0,0,1);
|
||||
|
||||
Info<< "globalAlignmentDirs " << globalAlignmentDirs << endl;
|
||||
|
||||
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::relaxPoints")
|
||||
<< "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 dualFace(verticesOnFace);
|
||||
|
||||
Cell_handle c = eit->first;
|
||||
Vertex_handle vA = c->vertex(eit->second);
|
||||
Vertex_handle vB = c->vertex(eit->third);
|
||||
|
||||
point dVA = topoint(vA->point());
|
||||
point dVB = topoint(vB->point());
|
||||
|
||||
Field<vector> alignmentDirs;
|
||||
|
||||
if
|
||||
(
|
||||
vA->alignmentDirections().size() > 0
|
||||
|| vB->alignmentDirections().size() > 0
|
||||
)
|
||||
{
|
||||
if
|
||||
(
|
||||
vA->alignmentDirections().size() == 3
|
||||
|| vB->alignmentDirections().size() == 3
|
||||
)
|
||||
{
|
||||
alignmentDirs.setSize(3);
|
||||
|
||||
if (vB->alignmentDirections().size() == 0)
|
||||
{
|
||||
alignmentDirs = vA->alignmentDirections();
|
||||
}
|
||||
else if (vA->alignmentDirections().size() == 0)
|
||||
{
|
||||
alignmentDirs = vB->alignmentDirections();
|
||||
}
|
||||
else if
|
||||
(
|
||||
vA->alignmentDirections().size() == 3
|
||||
&& vB->alignmentDirections().size() == 3
|
||||
)
|
||||
{
|
||||
forAll(vA->alignmentDirections(), aA)
|
||||
{
|
||||
const vector& a(vA->alignmentDirections()[aA]);
|
||||
|
||||
scalar maxDotProduct = 0.0;
|
||||
|
||||
forAll(vB->alignmentDirections(), aB)
|
||||
{
|
||||
const vector& b(vB->alignmentDirections()[aB]);
|
||||
|
||||
if (mag(a & b) > maxDotProduct)
|
||||
{
|
||||
maxDotProduct = mag(a & b);
|
||||
|
||||
alignmentDirs[aA] = a + sign(a & b)*b;
|
||||
|
||||
alignmentDirs[aA] /= mag(alignmentDirs[aA]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// One of the vertices has 3 alignments and the other
|
||||
// has 1
|
||||
|
||||
vector otherAlignment;
|
||||
|
||||
if (vA->alignmentDirections().size() == 3)
|
||||
{
|
||||
alignmentDirs = vA->alignmentDirections();
|
||||
|
||||
otherAlignment = vB->alignmentDirections()[0];
|
||||
}
|
||||
else
|
||||
{
|
||||
alignmentDirs = vB->alignmentDirections();
|
||||
|
||||
otherAlignment = vA->alignmentDirections()[0];
|
||||
}
|
||||
|
||||
label matchingDirection = -1;
|
||||
|
||||
scalar maxDotProduct = 0.0;
|
||||
|
||||
forAll(alignmentDirs, aD)
|
||||
{
|
||||
const vector& a(alignmentDirs[aD]);
|
||||
|
||||
if (mag(a & otherAlignment) > maxDotProduct)
|
||||
{
|
||||
maxDotProduct = mag(a & otherAlignment);
|
||||
|
||||
matchingDirection = aD;
|
||||
}
|
||||
}
|
||||
|
||||
vector& matchingAlignment =
|
||||
alignmentDirs[matchingDirection];
|
||||
|
||||
matchingAlignment = matchingAlignment
|
||||
+ sign(matchingAlignment & otherAlignment)
|
||||
*otherAlignment;
|
||||
|
||||
matchingAlignment /= mag(matchingAlignment);
|
||||
}
|
||||
|
||||
// vector midpoint = 0.5*(dVA + dVB);
|
||||
|
||||
// Info<< "midpoint " << midpoint
|
||||
// << nl << alignmentDirs
|
||||
// << nl << "v " << midpoint + alignmentDirs[0]
|
||||
// << nl << "v " << midpoint + alignmentDirs[1]
|
||||
// << nl << "v " << midpoint + alignmentDirs[2]
|
||||
// << nl << "v " << midpoint
|
||||
// << nl << "f 4 1"
|
||||
// << nl << "f 4 2"
|
||||
// << nl << "f 4 3"
|
||||
// << nl << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
alignmentDirs.setSize(1);
|
||||
|
||||
vector& alignmentDir = alignmentDirs[0];
|
||||
|
||||
if
|
||||
(
|
||||
vA->alignmentDirections().size() > 0
|
||||
&& vB->alignmentDirections().size() == 0
|
||||
)
|
||||
{
|
||||
alignmentDir = vA->alignmentDirections()[0];
|
||||
}
|
||||
else if
|
||||
(
|
||||
vA->alignmentDirections().size() == 0
|
||||
&& vB->alignmentDirections().size() > 0
|
||||
)
|
||||
{
|
||||
alignmentDir = vB->alignmentDirections()[0];
|
||||
}
|
||||
else
|
||||
{
|
||||
// Both vertices have an alignment
|
||||
|
||||
const vector& a(vA->alignmentDirections()[0]);
|
||||
|
||||
const vector& b(vB->alignmentDirections()[0]);
|
||||
|
||||
if (mag(a & b) < 1 - SMALL)
|
||||
{
|
||||
alignmentDirs.setSize(3);
|
||||
|
||||
alignmentDirs[0] = a + b;
|
||||
|
||||
alignmentDirs[1] = a - b;
|
||||
|
||||
alignmentDirs[2] =
|
||||
alignmentDirs[0] ^ alignmentDirs[1];
|
||||
|
||||
alignmentDirs /= mag(alignmentDirs);
|
||||
}
|
||||
else
|
||||
{
|
||||
alignmentDir = a + sign(a & b)*b;
|
||||
|
||||
alignmentDir /= mag(alignmentDir);
|
||||
}
|
||||
}
|
||||
if (alignmentDirs.size() ==1)
|
||||
{
|
||||
// Use the least similar of globalAlignmentDirs as the
|
||||
// 2nd alignment and then generate the third.
|
||||
|
||||
scalar minDotProduct = 1+SMALL;
|
||||
|
||||
alignmentDirs.setSize(3);
|
||||
|
||||
forAll(alignmentDirs, aD)
|
||||
{
|
||||
if
|
||||
(
|
||||
mag(alignmentDir & globalAlignmentDirs[aD])
|
||||
< minDotProduct
|
||||
)
|
||||
{
|
||||
minDotProduct = mag
|
||||
(
|
||||
alignmentDir & globalAlignmentDirs[aD]
|
||||
);
|
||||
|
||||
alignmentDirs[1] = globalAlignmentDirs[aD];
|
||||
}
|
||||
}
|
||||
|
||||
alignmentDirs[2] = alignmentDirs[0] ^ alignmentDirs[1];
|
||||
|
||||
alignmentDirs[2] /= mag(alignmentDirs[2]);
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
alignmentDirs = globalAlignmentDirs;
|
||||
}
|
||||
|
||||
// alignmentDirs found, now apply them
|
||||
|
||||
vector rAB = dVA - dVB;
|
||||
|
||||
scalar rABMag = mag(rAB);
|
||||
|
||||
forAll(alignmentDirs, aD)
|
||||
{
|
||||
vector& alignmentDir = alignmentDirs[aD];
|
||||
|
||||
if ((rAB & alignmentDir) < 0)
|
||||
{
|
||||
// swap the direction of the alignment so that has the
|
||||
// same sense as rAB
|
||||
alignmentDir *= -1;
|
||||
}
|
||||
|
||||
//scalar cosAcceptanceAngle = 0.743;
|
||||
scalar cosAcceptanceAngle = 0.67;
|
||||
|
||||
if (((rAB/rABMag) & alignmentDir) > cosAcceptanceAngle)
|
||||
{
|
||||
alignmentDir *= 0.5*controls_.minCellSize;
|
||||
|
||||
vector delta = alignmentDir - 0.5*rAB;
|
||||
|
||||
scalar faceArea = dualFace.mag(dualVertices);
|
||||
|
||||
delta *= faceAreaWeight(faceArea);
|
||||
|
||||
if (vA->internalPoint())
|
||||
{
|
||||
displacementAccumulator[vA->index()] += delta;
|
||||
}
|
||||
if (vB->internalPoint())
|
||||
{
|
||||
displacementAccumulator[vB->index()] += -delta;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
// Simple isotropic forcing using tension between the Delaunay vertex and
|
||||
// and the dual face centre
|
||||
// Info<< "globalAlignmentDirs " << globalAlignmentDirs << endl;
|
||||
|
||||
// for
|
||||
// (
|
||||
@ -962,7 +668,7 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
|
||||
// if (cc->cellIndex() < 0)
|
||||
// {
|
||||
// FatalErrorIn("Foam::CV3D::relaxPoints")
|
||||
// << "Dual face uses circumcenter defined by a "
|
||||
// << "Dual face uses circumcenter defined by a "
|
||||
// << " Delaunay tetrahedron with no internal "
|
||||
// << "or boundary points."
|
||||
// << exit(FatalError);
|
||||
@ -983,26 +689,321 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
|
||||
// point dVA = topoint(vA->point());
|
||||
// point dVB = topoint(vB->point());
|
||||
|
||||
// scalar faceArea = dualFace.mag(dualVertices);
|
||||
// Field<vector> alignmentDirs;
|
||||
|
||||
// point dualFaceCentre(dualFace.centre(dualVertices));
|
||||
|
||||
// if (vA->internalPoint())
|
||||
// if
|
||||
// (
|
||||
// vA->alignmentDirections().size() > 0
|
||||
// || vB->alignmentDirections().size() > 0
|
||||
// )
|
||||
// {
|
||||
// displacementAccumulator[vA->index()] +=
|
||||
// faceArea*(dualFaceCentre - dVA);
|
||||
// if
|
||||
// (
|
||||
// vA->alignmentDirections().size() == 3
|
||||
// || vB->alignmentDirections().size() == 3
|
||||
// )
|
||||
// {
|
||||
// alignmentDirs.setSize(3);
|
||||
|
||||
// weightAccumulator[vA->index()] += faceArea;
|
||||
// if (vB->alignmentDirections().size() == 0)
|
||||
// {
|
||||
// alignmentDirs = vA->alignmentDirections();
|
||||
// }
|
||||
// else if (vA->alignmentDirections().size() == 0)
|
||||
// {
|
||||
// alignmentDirs = vB->alignmentDirections();
|
||||
// }
|
||||
// else if
|
||||
// (
|
||||
// vA->alignmentDirections().size() == 3
|
||||
// && vB->alignmentDirections().size() == 3
|
||||
// )
|
||||
// {
|
||||
// forAll(vA->alignmentDirections(), aA)
|
||||
// {
|
||||
// const vector& a(vA->alignmentDirections()[aA]);
|
||||
|
||||
// scalar maxDotProduct = 0.0;
|
||||
|
||||
// forAll(vB->alignmentDirections(), aB)
|
||||
// {
|
||||
// const vector& b(vB->alignmentDirections()[aB]);
|
||||
|
||||
// if (mag(a & b) > maxDotProduct)
|
||||
// {
|
||||
// maxDotProduct = mag(a & b);
|
||||
|
||||
// alignmentDirs[aA] = a + sign(a & b)*b;
|
||||
|
||||
// alignmentDirs[aA] /= mag(alignmentDirs[aA]);
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// // One of the vertices has 3 alignments and the other
|
||||
// // has 1
|
||||
|
||||
// vector otherAlignment;
|
||||
|
||||
// if (vA->alignmentDirections().size() == 3)
|
||||
// {
|
||||
// alignmentDirs = vA->alignmentDirections();
|
||||
|
||||
// otherAlignment = vB->alignmentDirections()[0];
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// alignmentDirs = vB->alignmentDirections();
|
||||
|
||||
// otherAlignment = vA->alignmentDirections()[0];
|
||||
// }
|
||||
|
||||
// label matchingDirection = -1;
|
||||
|
||||
// scalar maxDotProduct = 0.0;
|
||||
|
||||
// forAll(alignmentDirs, aD)
|
||||
// {
|
||||
// const vector& a(alignmentDirs[aD]);
|
||||
|
||||
// if (mag(a & otherAlignment) > maxDotProduct)
|
||||
// {
|
||||
// maxDotProduct = mag(a & otherAlignment);
|
||||
|
||||
// matchingDirection = aD;
|
||||
// }
|
||||
// }
|
||||
|
||||
// vector& matchingAlignment =
|
||||
// alignmentDirs[matchingDirection];
|
||||
|
||||
// matchingAlignment = matchingAlignment
|
||||
// + sign(matchingAlignment & otherAlignment)
|
||||
// *otherAlignment;
|
||||
|
||||
// matchingAlignment /= mag(matchingAlignment);
|
||||
// }
|
||||
|
||||
// // vector midpoint = 0.5*(dVA + dVB);
|
||||
|
||||
// // Info<< "midpoint " << midpoint
|
||||
// // << nl << alignmentDirs
|
||||
// // << nl << "v " << midpoint + alignmentDirs[0]
|
||||
// // << nl << "v " << midpoint + alignmentDirs[1]
|
||||
// // << nl << "v " << midpoint + alignmentDirs[2]
|
||||
// // << nl << "v " << midpoint
|
||||
// // << nl << "f 4 1"
|
||||
// // << nl << "f 4 2"
|
||||
// // << nl << "f 4 3"
|
||||
// // << nl << endl;
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// alignmentDirs.setSize(1);
|
||||
|
||||
// vector& alignmentDir = alignmentDirs[0];
|
||||
|
||||
// if
|
||||
// (
|
||||
// vA->alignmentDirections().size() > 0
|
||||
// && vB->alignmentDirections().size() == 0
|
||||
// )
|
||||
// {
|
||||
// alignmentDir = vA->alignmentDirections()[0];
|
||||
// }
|
||||
// else if
|
||||
// (
|
||||
// vA->alignmentDirections().size() == 0
|
||||
// && vB->alignmentDirections().size() > 0
|
||||
// )
|
||||
// {
|
||||
// alignmentDir = vB->alignmentDirections()[0];
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// // Both vertices have an alignment
|
||||
|
||||
// const vector& a(vA->alignmentDirections()[0]);
|
||||
|
||||
// const vector& b(vB->alignmentDirections()[0]);
|
||||
|
||||
// if (mag(a & b) < 1 - SMALL)
|
||||
// {
|
||||
// alignmentDirs.setSize(3);
|
||||
|
||||
// alignmentDirs[0] = a + b;
|
||||
|
||||
// alignmentDirs[1] = a - b;
|
||||
|
||||
// alignmentDirs[2] =
|
||||
// alignmentDirs[0] ^ alignmentDirs[1];
|
||||
|
||||
// alignmentDirs /= mag(alignmentDirs);
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// alignmentDir = a + sign(a & b)*b;
|
||||
|
||||
// alignmentDir /= mag(alignmentDir);
|
||||
// }
|
||||
// }
|
||||
// if (alignmentDirs.size() ==1)
|
||||
// {
|
||||
// // Use the least similar of globalAlignmentDirs as the
|
||||
// // 2nd alignment and then generate the third.
|
||||
|
||||
// scalar minDotProduct = 1+SMALL;
|
||||
|
||||
// alignmentDirs.setSize(3);
|
||||
|
||||
// forAll(alignmentDirs, aD)
|
||||
// {
|
||||
// if
|
||||
// (
|
||||
// mag(alignmentDir & globalAlignmentDirs[aD])
|
||||
// < minDotProduct
|
||||
// )
|
||||
// {
|
||||
// minDotProduct = mag
|
||||
// (
|
||||
// alignmentDir & globalAlignmentDirs[aD]
|
||||
// );
|
||||
|
||||
// alignmentDirs[1] = globalAlignmentDirs[aD];
|
||||
// }
|
||||
// }
|
||||
|
||||
// alignmentDirs[2] = alignmentDirs[0] ^ alignmentDirs[1];
|
||||
|
||||
// alignmentDirs[2] /= mag(alignmentDirs[2]);
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// if (vB->internalPoint())
|
||||
// else
|
||||
// {
|
||||
// displacementAccumulator[vB->index()] +=
|
||||
// faceArea*(dualFaceCentre - dVB);
|
||||
// alignmentDirs = globalAlignmentDirs;
|
||||
// }
|
||||
|
||||
// weightAccumulator[vB->index()] += faceArea;
|
||||
// // alignmentDirs found, now apply them
|
||||
|
||||
// vector rAB = dVA - dVB;
|
||||
|
||||
// scalar rABMag = mag(rAB);
|
||||
|
||||
// forAll(alignmentDirs, aD)
|
||||
// {
|
||||
// vector& alignmentDir = alignmentDirs[aD];
|
||||
|
||||
// if ((rAB & alignmentDir) < 0)
|
||||
// {
|
||||
// // swap the direction of the alignment so that has the
|
||||
// // same sense as rAB
|
||||
// alignmentDir *= -1;
|
||||
// }
|
||||
|
||||
// //scalar cosAcceptanceAngle = 0.743;
|
||||
// scalar cosAcceptanceAngle = 0.67;
|
||||
|
||||
// if (((rAB/rABMag) & alignmentDir) > cosAcceptanceAngle)
|
||||
// {
|
||||
// alignmentDir *= 0.5*controls_.minCellSize;
|
||||
|
||||
// vector delta = alignmentDir - 0.5*rAB;
|
||||
|
||||
// scalar faceArea = dualFace.mag(dualVertices);
|
||||
|
||||
// delta *= faceAreaWeight(faceArea);
|
||||
|
||||
// if (vA->internalPoint())
|
||||
// {
|
||||
// displacementAccumulator[vA->index()] += delta;
|
||||
// }
|
||||
// if (vB->internalPoint())
|
||||
// {
|
||||
// displacementAccumulator[vB->index()] += -delta;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
// Simple isotropic forcing using tension between the Delaunay vertex and
|
||||
// and the dual face centre
|
||||
|
||||
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::relaxPoints")
|
||||
<< "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 dualFace(verticesOnFace);
|
||||
|
||||
Cell_handle c = eit->first;
|
||||
Vertex_handle vA = c->vertex(eit->second);
|
||||
Vertex_handle vB = c->vertex(eit->third);
|
||||
|
||||
point dVA = topoint(vA->point());
|
||||
point dVB = topoint(vB->point());
|
||||
|
||||
scalar faceArea = dualFace.mag(dualVertices);
|
||||
|
||||
// point dualFaceCentre(dualFace.centre(dualVertices));
|
||||
point dEMid = 0.5*(dVA + dVB);
|
||||
|
||||
if (vA->internalPoint())
|
||||
{
|
||||
displacementAccumulator[vA->index()] +=
|
||||
faceArea*(dEMid - dVA);
|
||||
//faceArea*(dualFaceCentre - dVA);
|
||||
|
||||
weightAccumulator[vA->index()] += faceArea;
|
||||
}
|
||||
if (vB->internalPoint())
|
||||
{
|
||||
displacementAccumulator[vB->index()] +=
|
||||
faceArea*(dEMid - dVB);
|
||||
//faceArea*(dualFaceCentre - dVB);
|
||||
|
||||
weightAccumulator[vB->index()] += faceArea;
|
||||
}
|
||||
}
|
||||
}
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
@ -1010,15 +1011,46 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
|
||||
|
||||
// Loop over all Delaunay vertices (Dual cells)
|
||||
|
||||
// Pick up all incident vertices to the Dv- check the number to see if this
|
||||
// is a reasonable result
|
||||
// for
|
||||
// (
|
||||
// Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
|
||||
// vit != finite_vertices_end();
|
||||
// vit++
|
||||
// )
|
||||
// {
|
||||
// if (vit->internalOrBoundaryPoint())
|
||||
// {
|
||||
// std::list<Vertex_handle> incidentVertices;
|
||||
|
||||
// Form the edge from the current Dv and iterate around the incident cells,
|
||||
// finding their circumcentres, then form (and store?) the dual face.
|
||||
// incident_vertices(vit, std::back_inserter(incidentVertices));
|
||||
|
||||
// for
|
||||
// (
|
||||
// std::list<Vertex_handle>::iterator ivit =
|
||||
// incidentVertices.begin();
|
||||
// ivit != incidentVertices.end();
|
||||
// ++ivit
|
||||
// )
|
||||
// {
|
||||
// // Pick up all incident vertices to the Dv - check the number to
|
||||
// // see if this is a reasonable result
|
||||
|
||||
// // Form the edge from the current Dv and iterate around the
|
||||
// // incident cells, finding their circumcentres, then form (and
|
||||
// // store?) the dual face.
|
||||
|
||||
// // Wait!
|
||||
|
||||
// // There is no sensible way in CGAL-3.3.1 to turn the Delaunay
|
||||
// // vertex in question and the incident vertex into an edge.
|
||||
// // There is, however, an incident_edges function in CGAL-3.4
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
//displacementAccumulator /= weightAccumulator;
|
||||
displacementAccumulator /= weightAccumulator;
|
||||
|
||||
vector totalDisp = sum(displacementAccumulator);
|
||||
scalar totalDist = sum(mag(displacementAccumulator));
|
||||
|
||||
Reference in New Issue
Block a user