Merge branch 'feature-adjoint-shapeOptimisation' into 'develop'

ENH: New adjont shape optimisation functionality

See merge request !307

applications/solvers/incompressible/adjointOptimisationFoam/adjointOptimisationFoam.C

@@ -55,10 +55,6 @@ int main(int argc, char *argv[])
 
     for (om++; !om.end(); om++)
     {
-        Info<< "* * * * * * * * * * * * * * * * * * *" << endl;
-        Info<< "Time = " << runTime.timeName()         << endl;
-        Info<< "* * * * * * * * * * * * * * * * * * *" << endl;
-
         if (om.update())
         {
             // Update design variables and solve all primal equations

applications/utilities/postProcessing/optimisation/cumulativeDisplacement/Make/files

+cumulativeDisplacement.C
+
+EXE = $(FOAM_APPBIN)/cumulativeDisplacement

applications/utilities/postProcessing/optimisation/cumulativeDisplacement/Make/options

+EXE_INC = \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -ldynamicFvMesh \
+    -ldynamicMesh \
+    -lmeshTools \
+    -lfiniteVolume

applications/utilities/postProcessing/optimisation/cumulativeDisplacement/createFields.H

+    pointField points0
+    (
+        pointIOField
+        (
+            IOobject
+            (
+                "points",
+                mesh.time().constant(),
+                polyMesh::meshSubDir,
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        )
+    );

applications/utilities/postProcessing/optimisation/cumulativeDisplacement/cumulativeDisplacement.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    cumulativeDisplacement
+
+Description
+    Computes and writes the difference between the mesh points in each time
+    instance and the ones in the constant folder. Additionally, the projection
+    of this difference to the normal point vectors of the initial mesh is also
+    written
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "emptyFvPatch.H"
+#include "coupledFvPatch.H"
+#include "pointMesh.H"
+#include "pointPatchField.H"
+#include "pointPatchFieldsFwd.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    timeSelector::addOptions();
+
+    #include "addRegionOption.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    instantList timeDirs = timeSelector::select0(runTime, args);
+    #include "createNamedMesh.H"
+    #include "createFields.H"
+
+    // polyPatch::pointNormals will give the wrong result for points
+    // belonging to multiple patches or patch-processorPatch intersections.
+    // Keeping a mesh-wide field to allow easy reduction using syncTools.
+    // A bit expensive? Better way?
+    vectorField pointNormals(mesh.nPoints(), vector::zero);
+    for (const fvPatch& patch : mesh.boundary())
+    {
+        //  Each local patch point belongs to these local patch faces.
+        //  Local numbering
+        const labelListList& patchPointFaces = patch.patch().pointFaces();
+        if (!isA<coupledFvPatch>(patch) && !isA<emptyFvPatch>(patch))
+        {
+            const labelList&  meshPoints = patch.patch().meshPoints();
+            const vectorField nf(patch.nf());
+            forAll(meshPoints, ppI)
+            {
+                const labelList& pointFaces = patchPointFaces[ppI];
+                forAll(pointFaces, pfI)
+                {
+                    const label& localFaceIndex = pointFaces[pfI];
+                    pointNormals[meshPoints[ppI]] += nf[localFaceIndex];
+                }
+            }
+        }
+    }
+    // Sum from all processors
+    syncTools::syncPointList
+    (
+        mesh, pointNormals, plusEqOp<vector>(), vector::zero
+    );
+    pointNormals /= mag(pointNormals) + VSMALL;
+
+    forAll(timeDirs, timeI)
+    {
+        runTime.setTime(timeDirs[timeI], timeI);
+        Info<< "Time = " << runTime.timeName() << endl;
+        mesh.readUpdate();
+
+        const pointMesh& pMesh(pointMesh::New(mesh));
+        // Point displacement projected to the
+        // unit point normal of the initial geometry
+        pointScalarField normalDisplacement
+        (
+            IOobject
+            (
+                "normalDisplacement",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            pMesh,
+            dimensionedScalar(dimless, Zero)
+        );
+        // Point displacement as a vector
+        pointVectorField displacement
+        (
+            IOobject
+            (
+                "displacement",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            pMesh,
+            dimensionedVector(dimless, Zero)
+        );
+
+        Info<< "Calculating cumulative mesh movement for time "
+            << runTime.timeName() << endl;
+        // Normal displacement
+        const pointField& meshPoints = mesh.points();
+        forAll(mesh.boundary(), pI)
+        {
+            const polyPatch& patch = mesh.boundaryMesh()[pI];
+            const labelList& localPoints = patch.meshPoints();
+            forAll(localPoints, ppI)
+            {
+                label pointI = localPoints[ppI];
+                normalDisplacement[pointI]  =
+                    (meshPoints[pointI] - points0[pointI]) 
+                  & pointNormals[pointI];
+            }
+        }
+        normalDisplacement.write();
+        // Vectorial displacement
+        displacement.primitiveFieldRef() = meshPoints - points0;
+        displacement.write();
+    }
+
+    // Print execution time
+    mesh.time().printExecutionTime(Info);
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/optimisation/writeActiveDesignVariables/Make/files

+writeActiveDesignVariables.C
+
+EXE = $(FOAM_APPBIN)/writeActiveDesignVariables

applications/utilities/preProcessing/optimisation/writeActiveDesignVariables/Make/options

+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/optimisation/adjointOptimisation/adjoint/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -ldynamicMesh \
+    -lfvMotionSolvers \
+    -ladjointOptimisation

applications/utilities/preProcessing/optimisation/writeActiveDesignVariables/writeActiveDesignVariables.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    writeActiveDesignVariables
+
+Description
+    Writes the active design variables based on the selected parameterisation 
+    scheme, as read from the meshMovement part of optimisationDict.
+    Keeps a back-up of the original optimisationDict in 
+    system/optimisationDict.org, as comments in the dict will be lost.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "optMeshMovement.H"
+#include "updateMethod.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    IOdictionary optDict
+    (
+        IOobject
+        (
+            "optimisationDict",
+            mesh.time().caseSystem(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+    // Back-up old optimisationDict, to maintain potential comments in it
+    if (Pstream::master())
+    {
+        cp(optDict.objectPath(), optDict.objectPath() + ".org");
+    }
+
+    // Construct mesh movement object and grab active design variables
+    // Will exit with error if not implemented for this type
+    const dictionary& movementDict = 
+        optDict.subDict("optimisation").subDict("meshMovement");
+    // Empty patch list will do 
+    labelList patchIDs(0);
+    autoPtr<optMeshMovement> movementPtr
+        (optMeshMovement::New(mesh, movementDict, patchIDs));
+    const labelList activeDesignVariables = 
+        movementPtr().getActiveDesignVariables();
+
+    // Construct update method to grab the type
+    dictionary& updateMethodDict = 
+        optDict.subDict("optimisation").subDict("updateMethod");
+    autoPtr<updateMethod> updMethod(updateMethod::New(mesh, updateMethodDict));
+
+    // Add to appropriate dictionary (creates it if it does not exist)
+    dictionary& coeffsDict = updateMethodDict.subDictOrAdd(updMethod().type());
+    coeffsDict.add<labelList>
+    (
+        "activeDesignVariables", 
+        activeDesignVariables, 
+        true
+    );
+    
+    // Write modified dictionary
+    optDict.regIOobject::writeObject
+    (
+        IOstream::ASCII,
+        IOstream::currentVersion,
+        IOstream::UNCOMPRESSED,
+        true
+    );
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/optimisation/writeMorpherCPs/Make/files

+writeMorpherCPs.C
+
+EXE = $(FOAM_APPBIN)/writeMorpherCPs

applications/utilities/preProcessing/optimisation/writeMorpherCPs/Make/options

+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/optimisation/adjointOptimisation/adjoint/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -ldynamicMesh \
+    -lfvMotionSolvers \
+    -ladjointOptimisation

applications/utilities/preProcessing/optimisation/writeMorpherCPs/writeMorpherCPs.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    writeMorpherCPs
+
+Description
+    Writes the NURBS3DVolume control points corresponding to dynamicMeshDict
+    to csv files. Parametric coordinates are not computed.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "NURBS3DVolume.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    IOdictionary dict
+    (
+        IOobject
+        (
+            "dynamicMeshDict",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+    const dictionary& coeffDict =
+        dict.subDict("volumetricBSplinesMotionSolverCoeffs");
+
+    wordList controlBoxes(coeffDict.get<wordList>("controlBoxes"));
+
+    forAll(controlBoxes, iNURB)
+    {
+        // Creating an object writes the control points in the 
+        // constructor
+        NURBS3DVolume::New
+        (
+            coeffDict.subDict(controlBoxes[iNURB]),
+            mesh,
+            false // do not compute parametric coordinates
+        );
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/optimisation/adjointOptimisation/adjoint/ATCModel/zeroATCcells/pointCells/pointCells.C

@@ -52,6 +52,9 @@ pointCells::pointCells
 :
     zeroATCcells(mesh, dict)
 {
+    boolList isZeroed(mesh_.nCells(), false);
+    labelList zeroedIDs(mesh_.nCells(), -1);
+    label i(0);
     forAll(mesh_.boundary(), patchI)
     {
         const fvPatch& patch = mesh_.boundary()[patchI];
@@ -65,7 +68,14 @@ pointCells::pointCells
                 for (const label pointI : meshPoints)
                 {
                     const labelList& pointCells = mesh_.pointCells()[pointI];
-                    zeroATCcells_.append(pointCells);
+                    for (const label cellI : pointCells)
+                    {
+                        if (!isZeroed[cellI])
+                        {
+                            zeroedIDs[i++] = cellI;
+                            isZeroed[cellI] = true;
+                        }
+                    }
                 }
             }
         }
@@ -76,14 +86,22 @@ pointCells::pointCells
         if (zoneID != -1)
         {
             const labelList& zoneCells = mesh_.cellZones()[zoneID];
-            zeroATCcells_.append(zoneCells);
+            for (const label cellI : zoneCells)
+            {
+                if (!isZeroed[cellI])
+                {
+                    zeroedIDs[i++] = cellI;
+                    isZeroed[cellI] = true;
+                }
+            }
         }
     }
-    /*
+    zeroedIDs.setSize(i);
+    zeroATCcells_ = zeroedIDs;
+
     label size = zeroATCcells_.size();
     reduce(size, sumOp<label>());
     Info<< "Zeroing ATC on "<< size << " cells" << nl << endl;
-    */
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/optimisation/adjointOptimisation/adjoint/Make/files

@@ -16,6 +16,7 @@ solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVars.C
 solvers/solverControl/solverControl/solverControl.C
 solvers/solverControl/SIMPLEControl/SIMPLEControl/SIMPLEControl.C
 solvers/solverControl/SIMPLEControl/singleRun/SIMPLEControlSingleRun.C
+solvers/solverControl/SIMPLEControl/optimisation/SIMPLEControlOpt.C
 
 /* SOLVERS */
 solvers/solver/solver.C
@@ -52,6 +53,8 @@ objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C
 objectives/incompressible/objectiveForce/objectiveForce.C
 objectives/incompressible/objectiveMoment/objectiveMoment.C
 objectives/incompressible/objectivePtLosses/objectivePtLosses.C
+objectives/incompressible/objectiveForceTarget/objectiveForceTarget.C
+objectives/incompressible/objectivePartialVolume/objectivePartialVolume.C
 
 /* OBJECTIVE MANAGER*/
 objectiveManager/objectiveManager/objectiveManager.C
@@ -87,8 +90,39 @@ adjointBoundaryConditions/adjointOutletVelocityFlux/adjointOutletVelocityFluxFvP
 /* DELTA BOUNDARY */
 deltaBoundary/deltaBoundary.C
 
+/* NURBS */
+NURBS=parameterization/NURBS
+$(NURBS)/NURBSbasis/NURBSbasis.C
+$(NURBS)/NURBS3DCurve/NURBS3DCurve.C
+$(NURBS)/NURBS3DSurface/NURBS3DSurface.C
+$(NURBS)/NURBS3DVolume/NURBS3DVolume/NURBS3DVolume.C
+$(NURBS)/NURBS3DVolume/cartesian/NURBS3DVolumeCartesian.C
+$(NURBS)/NURBS3DVolume/cylindrical/NURBS3DVolumeCylindrical.C
+$(NURBS)/NURBS3DVolume/volBSplinesBase/volBSplinesBase.C
+
+/* BEZIER */
+parameterization/Bezier/Bezier.C
+
+/* MOTION SOLVERS */
+dynamicMesh/motionSolver/volumetricBSplinesMotionSolver/volumetricBSplinesMotionSolver.C
+dynamicMesh/motionSolver/elasticityMotionSolver/elasticityMotionSolver.C
+dynamicMesh/motionSolver/laplacianMotionSolver/laplacianMotionSolver.C
+
+/* DISPLACEMENT METHOD */
+displacementMethod/displacementMethod/displacementMethod.C
+displacementMethod/displacementMethodvolumetricBSplinesMotionSolver/displacementMethodvolumetricBSplinesMotionSolver.C
+displacementMethod/displacementMethoddisplacementLaplacian/displacementMethoddisplacementLaplacian.C
+displacementMethod/displacementMethodvelocityLaplacian/displacementMethodvelocityLaplacian.C
+displacementMethod/displacementMethodelasticityMotionSolver/displacementMethodelasticityMotionSolver.C
+displacementMethod/displacementMethodlaplacianMotionSolver/displacementMethodlaplacianMotionSolver.C
+
+/* INTERPOLATION SCHEMES */
+interpolation/pointVolInterpolation/pointVolInterpolation.C
+
+
 /* ADJOINT SENSITIVITY */
 optimisation/adjointSensitivity/sensitivity/sensitivity.C
+optimisation/adjointSensitivity/shapeSensitivitiesBase/shapeSensitivitiesBase.C
 incoSens=optimisation/adjointSensitivity/incompressible
 $(incoSens)/adjointSensitivity/adjointSensitivityIncompressible.C
 $(incoSens)/adjointEikonalSolver/adjointEikonalSolverIncompressible.C
@@ -96,21 +130,49 @@ $(incoSens)/adjointMeshMovementSolver/adjointMeshMovementSolverIncompressible.C
 $(incoSens)/sensitivitySurface/sensitivitySurfaceIncompressible.C
 $(incoSens)/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C
 $(incoSens)/sensitivityMultiple/sensitivityMultipleIncompressible.C
+$(incoSens)/SIBase/SIBaseIncompressible.C
+$(incoSens)/FIBase/FIBaseIncompressible.C
+$(incoSens)/sensitivityVolBSplines/sensitivityVolBSplinesIncompressible.C
+$(incoSens)/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.C
+$(incoSens)/sensitivityBezier/sensitivityBezierIncompressible.C
+$(incoSens)/sensitivityBezierFI/sensitivityBezierFIIncompressible.C
 
 /* LINE SEARCH */
 optimisation/lineSearch/lineSearch/lineSearch.C
+optimisation/lineSearch/ArmijoConditions/ArmijoConditions.C
 optimisation/lineSearch/stepUpdate/stepUpdate/stepUpdate.C
+optimisation/lineSearch/stepUpdate/bisection/bisection.C
+optimisation/lineSearch/stepUpdate/quadratic/quadratic.C
 
 /* UPDATE METHOD */
-optimisation/updateMethod/updateMethod/updateMethod.C
-optimisation/updateMethod/constrainedOptimisationMethod/constrainedOptimisationMethod.C
+updateMethod=optimisation/updateMethod/
+$(updateMethod)/updateMethod/updateMethod.C
+$(updateMethod)/constrainedOptimisationMethod/constrainedOptimisationMethod.C
+$(updateMethod)/steepestDescent/steepestDescent.C
+$(updateMethod)/BFGS/BFGS.C
+$(updateMethod)/DBFGS/DBFGS.C
+$(updateMethod)/LBFGS/LBFGS.C
+$(updateMethod)/SR1/SR1.C
+$(updateMethod)/conjugateGradient/conjugateGradient.C
+$(updateMethod)/constraintProjection/constraintProjection.C
+$(updateMethod)/SQP/SQP.C
+
+/* OPT MESH MOVEMENT */
+optMeshMovement=optimisation/optMeshMovement