From fbe607811a18d71343dc42ae5fd4a7da2da4749d Mon Sep 17 00:00:00 2001
From: graham <graham.macpherson@strath.ac.uk>
Date: Wed, 21 Jan 2009 17:48:27 +0000
Subject: [PATCH] Adding relaxation calculation that adjusts to changes in
 runTime controls during a run.  Experimenting with primary alignment
 rotations to achieve boundary alignment.  Altering template depth back to 60
 after master merge conflict resolution.

---
 .../mesh/generation/CV3DMesher/CV3D.C         | 190 +++++++++++++++---
 .../mesh/generation/CV3DMesher/CV3DMesher.C   |  35 ++--
 wmake/rules/linux64Gcc/c++                    |   2 +-
 3 files changed, 177 insertions(+), 50 deletions(-)

diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3D.C b/applications/utilities/mesh/generation/CV3DMesher/CV3D.C
index 79b3f0424c..81cde7813b 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/CV3D.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/CV3D.C
@@ -228,6 +228,8 @@ void Foam::CV3D::setVertexAlignmentDirections()
             {
                 alignmentDirections.setSize(0);
             }
+
+            vit->alignmentDirections() = alignmentDirections;
         }
     }
 }
@@ -259,9 +261,9 @@ Foam::scalar Foam::CV3D::alignmentDistanceWeight(scalar dist) const
 
 Foam::scalar Foam::CV3D::faceAreaWeight(scalar faceArea) const
 {
-    scalar fl2 = 0.2;
+    scalar fl2 = 0.36;
 
-    scalar fu2 = 1.0;
+    scalar fu2 = 0.8;
 
     scalar m2 = controls_.minCellSize2;
 
@@ -533,11 +535,95 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
 
     dualVertices.setSize(dualVerti);
 
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // loop around the Delaunay edges to construct the dual faces.
     // Find the face-centre and use it to calculate the displacement vector
     // contribution to the Delaunay vertices (Dv) attached to the edge.  Add the
     // contribution to the running displacement vector of each Dv.
 
+    // 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());
+
+    //         point dualFaceCentre(dualFace.centre(dualVertices));
+
+    //         vector rAB = dVA - dVB;
+
+    //         scalar rABMag = mag(rAB);
+
+    //         scalar faceArea = dualFace.mag(dualVertices);
+
+    //         scalar directStiffness = 2.0;
+
+    //         scalar transverseStiffness = 0.0001;
+
+    //         scalar r0 = 0.9*controls_.minCellSize;
+
+    //         vector dA = -directStiffness*(1 - r0/rABMag)
+    //         *faceAreaWeight(faceArea)*rAB;
+
+    //         vector dT = transverseStiffness*faceAreaWeight(faceArea)
+    //         *(dualFaceCentre - 0.5*(dVA - dVB));
+
+    //         if (vA->internalPoint())
+    //         {
+    //             displacementAccumulator[vA->index()] += (dA + dT);
+    //         }
+    //         if (vB->internalPoint())
+    //         {
+    //             displacementAccumulator[vB->index()] += (-dA + dT);
+    //         }
+    //     }
+    // }
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Rotate faces that are sufficiently large and well enough aligned with the
+    // cell alignment direction(s)
+
     for
     (
         Triangulation::Finite_edges_iterator eit = finite_edges_begin();
@@ -563,7 +649,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);
@@ -584,37 +670,83 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
             point dVA = topoint(vA->point());
             point dVB = topoint(vB->point());
 
-            point dualFaceCentre(dualFace.centre(dualVertices));
-
-            vector rAB = dVA - dVB;
-
-            scalar rABMag = mag(rAB);
-
-            scalar faceArea = dualFace.mag(dualVertices);
-
-            scalar directStiffness = 2.0;
-
-            scalar transverseStiffness = 0.0001;
-
-            scalar r0 = 0.9*controls_.minCellSize;
-
-            vector dA = -directStiffness*(1 - r0/rABMag)
-            *faceAreaWeight(faceArea)*rAB;
-
-            vector dT = transverseStiffness*faceAreaWeight(faceArea)
-            *(dualFaceCentre - 0.5*(dVA - dVB));
-
-            if (vA->internalPoint())
+            if
+            (
+                vA->alignmentDirections().size() > 0
+             || vB->alignmentDirections().size() > 0
+            )
             {
-                displacementAccumulator[vA->index()] += (dA + dT);
-            }
-            if (vB->internalPoint())
-            {
-                displacementAccumulator[vB->index()] += (-dA + dT);
+                vector alignmentDir;
+
+                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
+
+                    alignmentDir = vA->alignmentDirections()[0]
+                    - vB->alignmentDirections()[0];
+
+                    if (mag(alignmentDir) < SMALL)
+                    {
+                        alignmentDir = vA->alignmentDirections()[0];
+                    }
+
+                    alignmentDir /= mag(alignmentDir);
+                }
+
+                vector rAB = dVA - dVB;
+
+                scalar rABMag = mag(rAB);
+
+                if ((rAB & alignmentDir) < 0)
+                {
+                    // swap the direction of the alignment so that has the same
+                    // sense as rAB
+                    alignmentDir *= -1;
+                }
+
+                scalar cosAcceptanceAngle = 0.743;
+
+                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;
+                    }
+                }
             }
         }
     }
 
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     vector totalDisp = sum(displacementAccumulator);
     scalar totalDist = sum(mag(displacementAccumulator));
 
diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C b/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C
index 60a60bb545..e8e93cc79f 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C
@@ -92,31 +92,16 @@ int main(int argc, char *argv[])
         mesh.boundaryConform();
     }
 
-    scalar relaxation =
-    mesh.meshingControls().relaxationFactorStart;
-
-    label nIterations = label
-    (
-        (runTime.endTime().value() - runTime.startTime().value())
-        /runTime.deltaT().value()
-    );
-
-    scalar relaxationDelta =
-    (
-        mesh.meshingControls().relaxationFactorStart
-      - mesh.meshingControls().relaxationFactorEnd
-    )
-    /max(nIterations, 1);
+    scalar relaxation = mesh.meshingControls().relaxationFactorStart;
 
     while (runTime.run())
     {
         runTime++;
 
         Info<< nl
-            << "Relaxation iteration " << runTime.timeIndex() << nl
-            << "~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
-
-        Info<< "relaxation = " << relaxation << endl;
+            << "Time = " << runTime.timeName()
+            << nl << "relaxation = " << relaxation
+            << endl;
 
         mesh.relaxPoints(relaxation);
 
@@ -124,7 +109,17 @@ int main(int argc, char *argv[])
         mesh.insertSurfacePointPairs();
         mesh.boundaryConform();
 
-        relaxation -= relaxationDelta;
+        relaxation -= (relaxation - mesh.meshingControls().relaxationFactorEnd)
+        /max
+        (
+            (runTime.endTime().value() - runTime.timeOutputValue())
+            /runTime.deltaT().value(),
+            1
+        );
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
     }
 
     mesh.write();
diff --git a/wmake/rules/linux64Gcc/c++ b/wmake/rules/linux64Gcc/c++
index 36b5c29b20..677b1e5bb4 100644
--- a/wmake/rules/linux64Gcc/c++
+++ b/wmake/rules/linux64Gcc/c++
@@ -6,7 +6,7 @@ CC          = g++ -m64
 
 include $(RULES)/c++$(WM_COMPILE_OPTION)
 
-ptFLAGS     = -DNoRepository -ftemplate-depth-40
+ptFLAGS     = -DNoRepository -ftemplate-depth-60
 
 c++FLAGS    = $(GFLAGS) $(c++WARN) $(c++OPT) $(c++DBUG) $(ptFLAGS) $(LIB_HEADER_DIRS) -fPIC