ENH: enable writing volumetric B-Splines control points in binary

as a step towards machine-accuracy continuation of the optimisation
loop.

Additionally, control points are now written under the time/uniform
folder, to be in-line with rest of the code structure for continuation.
As a side-effect, the controlPointsDefinition in
constant/dynamicMeshDict does not need to be changed to 'fromFile'
anymore in order to perform the continuation. The 'fromFile' option is
still valid if the user wants to supply the control points manually but,
as with all other controlPointsDefinitions, it will be disregarded if the
proper file exists under the time/uniform/volumetricBSplines folder.

src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/NURBS3DVolume/NURBS3DVolume.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -631,6 +631,19 @@ Foam::NURBS3DVolume::NURBS3DVolume
     bool computeParamCoors
 )
 :
+    localIOdictionary
+    (
+        IOobject
+        (
+            dict.dictName() + "cpsBsplines",
+            mesh.time().timeName(),
+            fileName("uniform")/fileName("volumetricBSplines"),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        word::null
+    ),
     mesh_(mesh),
     dict_(dict),
     name_(dict.dictName()),
@@ -640,7 +653,17 @@ Foam::NURBS3DVolume::NURBS3DVolume
     maxIter_(dict.getOrDefault<label>("maxIterations", 10)),
     tolerance_(dict.getOrDefault<scalar>("tolerance", 1.e-10)),
     nMaxBound_(dict.getOrDefault<scalar>("nMaxBoundIterations", 4)),
-    cps_(0),
+    cps_
+    (
+        found("controlPoints") ?
+        vectorField
+        (
+            "controlPoints",
+            *this,
+            basisU_.nCPs()*basisV_.nCPs()*basisW_.nCPs()
+        ) :
+        vectorField(0)
+    ),
     mapPtr_(nullptr),
     reverseMapPtr_(nullptr),
     parametricCoordinatesPtr_(nullptr),
@@ -737,10 +760,12 @@ Foam::NURBS3DVolume::NURBS3DVolume
            << exit(FatalError);
     }
 
-    // Define control points
-    controlPointsDefinition::New(*this);
+    // Construct control points, if not already read from file
+    if (cps_.empty())
+    {
+        controlPointsDefinition::New(*this);
+    }
     determineActiveDesignVariablesAndPoints();
-    writeCpsInDict();
 }
 
 
@@ -1382,7 +1407,6 @@ Foam::tmp<Foam::vectorField> Foam::NURBS3DVolume::computeNewPoints
     // Update control points position
     cps_ += controlPointsMovement;
     writeCps("cpsBsplines"+mesh_.time().timeName());
-    writeCpsInDict();
 
     // Compute new mesh points based on updated control points
     tmp<vectorField> tparameterizedPoints = coordinates(paramCoors);
@@ -1411,7 +1435,8 @@ Foam::tmp<Foam::vectorField> Foam::NURBS3DVolume::computeNewPoints
 Foam::tmp<Foam::vectorField> Foam::NURBS3DVolume::computeNewBoundaryPoints
 (
     const vectorField& controlPointsMovement,
-    const labelList& patchesToBeMoved
+    const labelList& patchesToBeMoved,
+    const bool updateCPs
 )
 {
     // Get parametric coordinates
@@ -1420,8 +1445,10 @@ Foam::tmp<Foam::vectorField> Foam::NURBS3DVolume::computeNewBoundaryPoints
     // Update control points position
     cps_ += controlPointsMovement;
 
-    writeCps("cpsBsplines"+mesh_.time().timeName());
-    writeCpsInDict();
+    if (updateCPs)
+    {
+        writeCps("cpsBsplines"+mesh_.time().timeName());
+    }
 
     // Return field. Initialized with current mesh points
     tmp<vectorField> tnewPoints(new vectorField(mesh_.points()));
@@ -1452,8 +1479,17 @@ Foam::tmp<Foam::vectorField> Foam::NURBS3DVolume::computeNewBoundaryPoints
         }
     }
 
-    // Update coordinates in the local system based on the cartesian points
-    updateLocalCoordinateSystem(newPoints);
+    if (updateCPs)
+    {
+        // Update coordinates in the local system based on the cartesian points
+        updateLocalCoordinateSystem(newPoints);
+    }
+    else
+    {
+        // Move control points to their initial position
+        cps_ -= controlPointsMovement;
+    }
+
     DebugInfo
         << "Max mesh movement equal to "
         << gMax(mag(newPoints - mesh_.points())) << endl;
@@ -1492,7 +1528,7 @@ void Foam::NURBS3DVolume::setControlPoints(const vectorField& newCps)
 void Foam::NURBS3DVolume::boundControlPointMovement
 (
     vectorField& controlPointsMovement
-)
+) const
 {
     forAll(controlPointsMovement, cpI)
     {
@@ -1732,45 +1768,27 @@ Foam::tmp<Foam::volTensorField> Foam::NURBS3DVolume::getDxCellsDb
 Foam::label Foam::NURBS3DVolume::nUSymmetry() const
 {
     label nU(basisU_.nCPs());
-    if (nU % 2 == 0)
-    {
-        nU /=2;
-    }
-    else
-    {
-        nU = (nU - 1)/2 + 1;
-    }
-    return nU;
+    return label(nU % 2 == 0 ? 0.5*nU : (nU - 1)/2 + 1);
 }
 
 
 Foam::label Foam::NURBS3DVolume::nVSymmetry() const
 {
     label nV(basisV_.nCPs());
-    if (nV % 2 == 0)
-    {
-        nV /=2;
-    }
-    else
-    {
-        nV = (nV - 1)/2 + 1;
-    }
-    return nV;
+    return label(nV % 2 == 0 ? 0.5*nV : 0.5*(nV - 1) + 1);
 }
 
 
 Foam::label Foam::NURBS3DVolume::nWSymmetry() const
 {
     label nW(basisW_.nCPs());
-    if (nW % 2 == 0)
-    {
-        nW /=2;
-    }
-    else
-    {
-        nW = (nW - 1)/2 + 1;
-    }
-    return nW;
+    return label(nW % 2 == 0 ? 0.5*nW : 0.5*(nW - 1) + 1);
+}
+
+
+Foam::Vector<Foam::label> Foam::NURBS3DVolume::nSymmetry() const
+{
+    return Vector<label>(nUSymmetry(), nVSymmetry(), nWSymmetry());
 }
 
 
@@ -1824,37 +1842,17 @@ void Foam::NURBS3DVolume::writeCps
 }
 
 
-void Foam::NURBS3DVolume::writeCpsInDict() const
-{
-    IOdictionary cpsDict
-    (
-        IOobject
-        (
-            name_ + "cpsBsplines" + mesh_.time().timeName(),
-            mesh_.time().caseConstant(),
-            cpsFolder_,
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        )
-    );
-
-    cpsDict.add("controlPoints", cps_);
-
-    // Always write in ASCII, but allow compression
-    cpsDict.regIOobject::writeObject
-    (
-        IOstreamOption(IOstream::ASCII, mesh_.time().writeCompression()),
-        true
-    );
-}
-
-
 void Foam::NURBS3DVolume::write() const
 {
     parametricCoordinatesPtr_().write();
 }
 
 
+bool Foam::NURBS3DVolume::writeData(Ostream& os) const
+{
+    cps_.writeEntry("controlPoints", os);
+    return true;
+}
+
+
 // ************************************************************************* //

src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/NURBS3DVolume/NURBS3DVolume.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -56,6 +56,7 @@ SourceFiles
 
 #include "NURBSbasis.H"
 #include "boolVector.H"
+#include "localIOdictionary.H"
 #include "pointMesh.H"
 #include "pointPatchField.H"
 #include "pointPatchFieldsFwd.H"
@@ -71,6 +72,8 @@ namespace Foam
 \*---------------------------------------------------------------------------*/
 
 class NURBS3DVolume
+:
+    public localIOdictionary
 {
 protected:
 
@@ -371,7 +374,8 @@ public:
             tmp<vectorField> computeNewBoundaryPoints
             (
                 const vectorField& controlPointsMovement,
-                const labelList& patchesToBeMoved
+                const labelList& patchesToBeMoved,
+                const bool moveCPs = true
             );
 
 
@@ -391,7 +395,7 @@ public:
             void boundControlPointMovement
             (
                 vectorField& controlPointsMovement
-            );
+            ) const;
 
             //- Compute max. displacement at the boundary
             scalar computeMaxBoundaryDisplacement
@@ -433,6 +437,10 @@ public:
                 //- Get number of variables if CPs are moved symmetrically in W
                 label nWSymmetry() const;
 
+                //- Get number of variables per direction,
+                //- if CPs are moved symmetrically
+                Vector<label> nSymmetry() const;
+
 
             // Inline access functions
 
@@ -466,6 +474,9 @@ public:
 
                 inline const NURBSbasis& basisW() const;
 
+                //- Get number of control points per direction
+                inline Vector<label> nCPsPerDirection() const;
+
                 //- Get mesh
                 inline const fvMesh& mesh() const;
 
@@ -483,12 +494,11 @@ public:
                 const bool transform = true
             ) const;
 
-            //- Write control points on the local coordinate system.
-            //  For continuation
-            void writeCpsInDict() const;
-
             //- Write parametric coordinates
             void write() const;
+
+            //- Write the control points to support restart
+            virtual bool writeData(Ostream& os) const;
 };
 
 

src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DVolume/volBSplinesBase/volBSplinesBase.C

@@ -315,8 +315,7 @@ void Foam::volBSplinesBase::writeControlPoints() const
 {
     for (const NURBS3DVolume& box : volume_)
     {
-        box.writeCps("cpsBsplines"+mesh_.time().timeName());
-        box.writeCpsInDict();
+        box.writeCps("cpsBsplines" + mesh_.time().timeName());
     }
 }