diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3D.C b/applications/utilities/mesh/generation/CV3DMesher/CV3D.C
index 4e30845a92a920fa8ab7e382e4b0fc39072cbc67..cfaa91ae1be8930af62299519ed28c0aff1f2d55 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/CV3D.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/CV3D.C
@@ -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,305 +629,305 @@ 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;
+ // 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;
+ // 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;
+ // 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);
- }
+ // 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.append(cc->cellIndex());
+ // }
+ // } while (++cc != ccStart);
- verticesOnFace.shrink();
+ // verticesOnFace.shrink();
- face dualFace(verticesOnFace);
+ // face dualFace(verticesOnFace);
- Cell_handle c = eit->first;
- Vertex_handle vA = c->vertex(eit->second);
- Vertex_handle vB = c->vertex(eit->third);
+ // 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());
+ // point dVA = topoint(vA->point());
+ // point dVB = topoint(vB->point());
- Field<vector> alignmentDirs;
+ // 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
+ // (
+ // 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]);
+ // 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;
+ // scalar maxDotProduct = 0.0;
- forAll(vB->alignmentDirections(), aB)
- {
- const vector& b(vB->alignmentDirections()[aB]);
+ // forAll(vB->alignmentDirections(), aB)
+ // {
+ // const vector& b(vB->alignmentDirections()[aB]);
- if (mag(a & b) > maxDotProduct)
- {
- maxDotProduct = mag(a & b);
+ // if (mag(a & b) > maxDotProduct)
+ // {
+ // maxDotProduct = mag(a & b);
- alignmentDirs[aA] = a + sign(a & b)*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
+ // alignmentDirs[aA] /= mag(alignmentDirs[aA]);
+ // }
+ // }
+ // }
+ // }
+ // else
+ // {
+ // // One of the vertices has 3 alignments and the other
+ // // has 1
- vector otherAlignment;
+ // vector otherAlignment;
- if (vA->alignmentDirections().size() == 3)
- {
- alignmentDirs = vA->alignmentDirections();
+ // if (vA->alignmentDirections().size() == 3)
+ // {
+ // alignmentDirs = vA->alignmentDirections();
- otherAlignment = vB->alignmentDirections()[0];
- }
- else
- {
- alignmentDirs = vB->alignmentDirections();
+ // otherAlignment = vB->alignmentDirections()[0];
+ // }
+ // else
+ // {
+ // alignmentDirs = vB->alignmentDirections();
- otherAlignment = vA->alignmentDirections()[0];
- }
+ // otherAlignment = vA->alignmentDirections()[0];
+ // }
- label matchingDirection = -1;
+ // label matchingDirection = -1;
- scalar maxDotProduct = 0.0;
+ // scalar maxDotProduct = 0.0;
- forAll(alignmentDirs, aD)
- {
- const vector& a(alignmentDirs[aD]);
+ // forAll(alignmentDirs, aD)
+ // {
+ // const vector& a(alignmentDirs[aD]);
- if (mag(a & otherAlignment) > maxDotProduct)
- {
- maxDotProduct = mag(a & otherAlignment);
+ // if (mag(a & otherAlignment) > maxDotProduct)
+ // {
+ // maxDotProduct = mag(a & otherAlignment);
- matchingDirection = aD;
- }
- }
+ // matchingDirection = aD;
+ // }
+ // }
- vector& matchingAlignment =
- alignmentDirs[matchingDirection];
+ // vector& matchingAlignment =
+ // alignmentDirs[matchingDirection];
- matchingAlignment = matchingAlignment
- + sign(matchingAlignment & otherAlignment)
- *otherAlignment;
+ // matchingAlignment = matchingAlignment
+ // + sign(matchingAlignment & otherAlignment)
+ // *otherAlignment;
- matchingAlignment /= mag(matchingAlignment);
- }
+ // matchingAlignment /= mag(matchingAlignment);
+ // }
- // vector midpoint = 0.5*(dVA + dVB);
+ // // 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);
- // 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];
- 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
- 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& a(vA->alignmentDirections()[0]);
+ // const vector& b(vB->alignmentDirections()[0]);
- const vector& b(vB->alignmentDirections()[0]);
+ // if (mag(a & b) < 1 - SMALL)
+ // {
+ // alignmentDirs.setSize(3);
- if (mag(a & b) < 1 - SMALL)
- {
- alignmentDirs.setSize(3);
+ // alignmentDirs[0] = a + b;
- alignmentDirs[0] = a + b;
+ // alignmentDirs[1] = a - b;
- alignmentDirs[1] = a - b;
+ // alignmentDirs[2] =
+ // alignmentDirs[0] ^ alignmentDirs[1];
- alignmentDirs[2] =
- alignmentDirs[0] ^ alignmentDirs[1];
+ // alignmentDirs /= mag(alignmentDirs);
+ // }
+ // else
+ // {
+ // alignmentDir = a + sign(a & b)*b;
- 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.
- 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;
- scalar minDotProduct = 1+SMALL;
+ // alignmentDirs.setSize(3);
- alignmentDirs.setSize(3);
+ // forAll(alignmentDirs, aD)
+ // {
+ // if
+ // (
+ // mag(alignmentDir & globalAlignmentDirs[aD])
+ // < minDotProduct
+ // )
+ // {
+ // minDotProduct = mag
+ // (
+ // alignmentDir & globalAlignmentDirs[aD]
+ // );
- forAll(alignmentDirs, aD)
- {
- if
- (
- mag(alignmentDir & globalAlignmentDirs[aD])
- < minDotProduct
- )
- {
- minDotProduct = mag
- (
- alignmentDir & globalAlignmentDirs[aD]
- );
-
- alignmentDirs[1] = globalAlignmentDirs[aD];
- }
- }
+ // alignmentDirs[1] = globalAlignmentDirs[aD];
+ // }
+ // }
- alignmentDirs[2] = alignmentDirs[0] ^ alignmentDirs[1];
+ // alignmentDirs[2] = alignmentDirs[0] ^ alignmentDirs[1];
- alignmentDirs[2] /= mag(alignmentDirs[2]);
- }
- }
- }
- else
- {
- alignmentDirs = globalAlignmentDirs;
- }
+ // alignmentDirs[2] /= mag(alignmentDirs[2]);
+ // }
+ // }
+ // }
+ // else
+ // {
+ // alignmentDirs = globalAlignmentDirs;
+ // }
- // alignmentDirs found, now apply them
+ // // alignmentDirs found, now apply them
- vector rAB = dVA - dVB;
+ // vector rAB = dVA - dVB;
- scalar rABMag = mag(rAB);
+ // scalar rABMag = mag(rAB);
- forAll(alignmentDirs, aD)
- {
- vector& alignmentDir = alignmentDirs[aD];
+ // 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;
- }
+ // 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;
+ // //scalar cosAcceptanceAngle = 0.743;
+ // scalar cosAcceptanceAngle = 0.67;
- if (((rAB/rABMag) & alignmentDir) > cosAcceptanceAngle)
- {
- alignmentDir *= 0.5*controls_.minCellSize;
+ // if (((rAB/rABMag) & alignmentDir) > cosAcceptanceAngle)
+ // {
+ // alignmentDir *= 0.5*controls_.minCellSize;
- vector delta = alignmentDir - 0.5*rAB;
+ // vector delta = alignmentDir - 0.5*rAB;
- scalar faceArea = dualFace.mag(dualVertices);
+ // scalar faceArea = dualFace.mag(dualVertices);
- delta *= faceAreaWeight(faceArea);
+ // delta *= faceAreaWeight(faceArea);
- if (vA->internalPoint())
- {
- displacementAccumulator[vA->index()] += delta;
- }
- if (vB->internalPoint())
- {
- displacementAccumulator[vB->index()] += -delta;
- }
- }
- }
- }
- }
+ // if (vA->internalPoint())
+ // {
+ // displacementAccumulator[vA->index()] += delta;
+ // }
+ // if (vB->internalPoint())
+ // {
+ // displacementAccumulator[vB->index()] += -delta;
+ // }
+ // }
+ // }
+ // }
+ // }
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -937,72 +935,75 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
// 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;
+ 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;
+ 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);
- // }
+ 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.append(cc->cellIndex());
+ }
+ } while (++cc != ccStart);
- // verticesOnFace.shrink();
+ verticesOnFace.shrink();
- // face dualFace(verticesOnFace);
+ face dualFace(verticesOnFace);
- // Cell_handle c = eit->first;
- // Vertex_handle vA = c->vertex(eit->second);
- // Vertex_handle vB = c->vertex(eit->third);
+ 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());
+ point dVA = topoint(vA->point());
+ point dVB = topoint(vB->point());
- // scalar faceArea = dualFace.mag(dualVertices);
+ scalar faceArea = dualFace.mag(dualVertices);
- // point dualFaceCentre(dualFace.centre(dualVertices));
+ // point dualFaceCentre(dualFace.centre(dualVertices));
+ point dEMid = 0.5*(dVA + dVB);
- // if (vA->internalPoint())
- // {
- // displacementAccumulator[vA->index()] +=
- // faceArea*(dualFaceCentre - dVA);
+ if (vA->internalPoint())
+ {
+ displacementAccumulator[vA->index()] +=
+ faceArea*(dEMid - dVA);
+ //faceArea*(dualFaceCentre - dVA);
- // weightAccumulator[vA->index()] += faceArea;
- // }
- // if (vB->internalPoint())
- // {
- // displacementAccumulator[vB->index()] +=
- // faceArea*(dualFaceCentre - dVB);
+ weightAccumulator[vA->index()] += faceArea;
+ }
+ if (vB->internalPoint())
+ {
+ displacementAccumulator[vB->index()] +=
+ faceArea*(dEMid - dVB);
+ //faceArea*(dualFaceCentre - dVB);
- // weightAccumulator[vB->index()] += faceArea;
- // }
- // }
- // }
+ 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));