From 4a87e2d0b4ec2e6a13adf7370119cde407bbe8d9 Mon Sep 17 00:00:00 2001
From: Mattijs Janssens <ext-mjanssens@esi-group.com>
Date: Wed, 26 May 2021 12:22:24 +0000
Subject: [PATCH] INT: splineEdge: allowing usage in extrudeMesh. See #1983.

---
 CONTRIBUTORS.md                               |  2 +-
 .../blockEdges/splineEdge/CatmullRomSpline.C  | 96 ++++++++++++++++---
 .../blockEdges/splineEdge/CatmullRomSpline.H  | 12 ++-
 tutorials/mesh/extrudeMesh/polyline/Allrun    | 12 ++-
 .../polyline/system/extrudeMeshDict           | 20 +++-
 5 files changed, 125 insertions(+), 17 deletions(-)

diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md
index 30d191014a..c3ca4001ef 100644
--- a/CONTRIBUTORS.md
+++ b/CONTRIBUTORS.md
@@ -50,7 +50,7 @@ It is likely incomplete...
 - Norbert Weber
 - Henry Weller
 - Niklas Wikstrom
+- Guanyang Xue
 - Thorsten Zirwes
 
-
 <!----------------------------------------------------------------------------->
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C b/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C
index bfa0aaf27b..681b6c0316 100644
--- a/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C
+++ b/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C
@@ -27,6 +27,56 @@ License
 
 #include "CatmullRomSpline.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::CatmullRomSpline::derivative
+(
+    const label segment,
+    const scalar mu
+) const
+{
+    const point& p0 = points()[segment];
+    const point& p1 = points()[segment+1];
+
+    // determine the end points
+    point e0;
+    point e1;
+
+    if (segment == 0)
+    {
+        // end: simple reflection
+        e0 = 2*p0 - p1;
+    }
+    else
+    {
+        e0 = points()[segment-1];
+    }
+
+    if (segment+1 == nSegments())
+    {
+        // end: simple reflection
+        e1 = 2*p1 - p0;
+    }
+    else
+    {
+        e1 = points()[segment+2];
+    }
+    const point derivativePoint
+    (
+        0.5 *
+        (
+            (-e0 + p1)
+          + mu *
+            (
+                2 * (2*e0 - 5*p0 + 4*p1 - e1)
+              + mu * 3 * (-e0 + 3*p0 - 3*p1 + e1)
+            )
+        )
+    );
+    return mag(derivativePoint);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::CatmullRomSpline::CatmullRomSpline
@@ -114,27 +164,51 @@ Foam::point Foam::CatmullRomSpline::position
     }
 
 
-    return 0.5 *
-    (
-        ( 2*p0 )
-      + mu *
+    return
+        0.5 *
         (
-            ( -e0 + p1 )
+            (2*p0)
           + mu *
             (
-                ( 2*e0 - 5*p0 + 4*p1 - e1 )
+                (-e0 + p1)
               + mu *
-                ( -e0 + 3*p0 - 3*p1 + e1 )
+                (
+                    (2*e0 - 5*p0 + 4*p1 - e1)
+                  + mu*(-e0 + 3*p0 - 3*p1 + e1)
+                )
             )
-        )
-    );
+        );
 }
 
 
 Foam::scalar Foam::CatmullRomSpline::length() const
 {
-    NotImplemented;
-    return 1;
+    const solveScalar xi[5]=
+    {
+        -0.9061798459386639927976,
+        -0.5384693101056830910363,
+        0,
+        0.5384693101056830910363,
+        0.9061798459386639927976
+    };
+    const solveScalar wi[5]=
+    {
+        0.2369268850561890875143,
+        0.4786286704993664680413,
+        0.5688888888888888888889,
+        0.4786286704993664680413,
+        0.2369268850561890875143
+    };
+    scalar sum=0;
+    for (label segment=0;segment<nSegments();segment++)
+    {
+        for (int i=0;i<5;i++)
+        {
+            sum+=wi[i]*derivative(segment,(xi[i]+1.0)/2.0)/2.0;
+        }
+    }
+
+    return sum;
 }
 
 
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.H b/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.H
index 27093f9cbd..a102584182 100644
--- a/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.H
+++ b/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -80,6 +80,15 @@ class CatmullRomSpline
 :
     public polyLine
 {
+
+    //- Derivative of spline
+    scalar derivative
+    (
+        const label segment,
+        const scalar mu
+    ) const;
+
+
 public:
 
     // Constructors
@@ -103,7 +112,6 @@ public:
         point position(const label segment, const scalar lambda) const;
 
         //- The length of the curve
-        //  \note NotImplemented
         scalar length() const;
 };
 
diff --git a/tutorials/mesh/extrudeMesh/polyline/Allrun b/tutorials/mesh/extrudeMesh/polyline/Allrun
index ee1f26eec4..86a89615c4 100755
--- a/tutorials/mesh/extrudeMesh/polyline/Allrun
+++ b/tutorials/mesh/extrudeMesh/polyline/Allrun
@@ -3,12 +3,20 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
-runApplication extrudeMesh
+# Extrude using arcs
+USE_ARC=true runApplication -s arc extrudeMesh
 
 # For output fields from checkMesh
 mkdir -p 1
+runApplication -s arc checkMesh -writeAllFields -time 1
 
-runApplication checkMesh -writeAllFields
+
+# Extrude using splines
+USE_ARC=false runApplication -s spline extrudeMesh
+
+# For output fields from checkMesh
+mkdir -p 2
+runApplication -s spline checkMesh -writeAllFields -time 2
 
 paraFoam -touch -vtk
 
diff --git a/tutorials/mesh/extrudeMesh/polyline/system/extrudeMeshDict b/tutorials/mesh/extrudeMesh/polyline/system/extrudeMeshDict
index 03accb0e5e..69f5f0e97e 100644
--- a/tutorials/mesh/extrudeMesh/polyline/system/extrudeMeshDict
+++ b/tutorials/mesh/extrudeMesh/polyline/system/extrudeMeshDict
@@ -42,7 +42,9 @@ polylineCoeffs
         ( 1.5 -1 1.183974596 )
     );
 
-    edges 9
+#if ${USE_ARC:-true}
+    //- Using arcs
+    edges
     (
         line 0 1
         arc  1 2 ( 0 -0.087867966 0.962132034 )
@@ -54,6 +56,22 @@ polylineCoeffs
         arc  7 8 ( 0.976794919 -1 2.009807621 )
         line 8 9
     );
+#else
+    //- Using spline (different shape):
+    edges
+    (
+        spline 0 9
+        (
+            ( 0 0 0.05 )
+            ( 0 0 0.5 )
+            ( 0.5 0  0.5 )
+            ( 0.5 0  0 )
+            ( 1.5 0  0 )
+            ( 1.5 -1 0 )
+            ( 1.5 -1 1.183974596 )
+        )
+    );
+#endif
 
     toleranceCheck  1e-6;
 }