NURBS3DVolume.H

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2007-2021 PCOpt/NTUA
    Copyright (C) 2013-2021 FOSS GP
    Copyright (C) 2019-2020 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/>.

Class
    Foam::NURBS3DVolume

Description
    NURBS3DVolume morpher. Includes support functions for gradient computations
    Base class providing support for different coordinate systems

    Reference:
    \verbatim
        For a short introduction to a volumetric B-Splines morpher and its use
        in shape optimisation

            Papoutsis-Kiachagias, E., Magoulas, N., Mueller, J.,
            Othmer, C., & Giannakoglou, K. (2015).
            Noise reduction in car aerodynamics using a surrogate objective
            function and the continuous adjoint method with wall functions.
            Computers & Fluids, 122, 223-232.
            http://doi.org/10.1016/j.compfluid.2015.09.002
    \endverbatim

SourceFiles
    NURBS3DVolume.C
    NURBS3DVolumeI.H

\*---------------------------------------------------------------------------*/

#ifndef NURBS3DVolume_H
#define NURBS3DVolume_H

#include "NURBSbasis.H"
#include "boolVector.H"
#include "localIOdictionary.H"
#include "pointMesh.H"
#include "pointPatchField.H"
#include "pointPatchFieldsFwd.H"
#include "boundaryFieldsFwd.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
                        Class NURBS3DVolume Declaration
\*---------------------------------------------------------------------------*/

class NURBS3DVolume
:
    public localIOdictionary
{
protected:

    // Protected Data

        typedef List<boolVector> boolVectorList;

        const fvMesh& mesh_;
        dictionary dict_;
        word name_;
        //- NURBS basis functions
        NURBSbasis basisU_;
        NURBSbasis basisV_;
        NURBSbasis basisW_;

        //- Max iterations for Newton-Raphson
        label maxIter_;

        //- Tolerance for Newton-Raphson
        scalar tolerance_;

        //- How many times to bound parametric coordinates until deciding it is
        //- outside the box
        label nMaxBound_;

        //- The volumetric B-Splines control points
        vectorField cps_;

        //- Map of points-in-box to mesh points
        autoPtr<labelList> mapPtr_;

        //- Map of mesh points to points-in-box
        //  Return -1 if point not within the box
        autoPtr<labelList> reverseMapPtr_;

        //- Parametric coordinates of pointsInBox
        autoPtr<pointVectorField> parametricCoordinatesPtr_;

        //- Coordinates in the local system for which CPs are defined
        vectorField localSystemCoordinates_;

        //- Confine movement in certain directions and control points. Refers
        //- to the local system
        label confineUMovement_;

        label confineVMovement_;

        label confineWMovement_;

        label confineBoundaryControlPoints_;

        //- Which movement components to freeze in each plane
        boolVectorList confineUMinCPs_;

        boolVectorList confineUMaxCPs_;

        boolVectorList confineVMinCPs_;

        boolVectorList confineVMaxCPs_;

        boolVectorList confineWMinCPs_;

        boolVectorList confineWMaxCPs_;

        //- Which of the cps are moved in an optimisation
        boolList activeControlPoints_;

        //- Which design variables are changed in an optimisation
        boolList activeDesignVariables_;

        //- Folder to store control points
        string cpsFolder_;

        //- Read parametric coordinates from file if present in the folder
        bool readStoredData_;


    // Protected Member Functions

        //- Find points within control points box
        void findPointsInBox(const vectorField& meshPoints);

        //- Compute parametric coordinates in order to match a given set
        //- of coordinates, based on the cps of the class
        // Uses a Newton-Raphson loop.
        // Argument is the points residing in the box
        void computeParametricCoordinates(const vectorField& points);

        void computeParametricCoordinates(tmp<vectorField> tPoints);

        //- Bound components to certain limits
        bool bound
        (
            vector& vec,
            scalar minValue = 1e-7,
            scalar maxValue = 0.999999
        );

        //- Create lists with active design variables and control points
        void determineActiveDesignVariablesAndPoints();

        //- Confine movement in boundary control points if necessary
        void confineBoundaryControlPoints();

        //- Confine control point movement to maintain user-defined continuity
        void continuityRealatedConfinement();

        //- Confine movement in all control points for user-defined directions
        void confineControlPointsDirections();

        //- Confine all three movements for a prescribed control point
        void confineControlPoint(const label cpI);

        //- Confine specific movements for a prescribed control point
        void confineControlPoint(const label cpI, const boolVector&);

        //- Create folders to store cps and derivatives
        void makeFolders();

        //- Transform a point from its coordinate system to a cartesian system
        virtual vector transformPointToCartesian
        (
            const vector& localCoordinates
        ) const = 0;

        //- Transformation tensor for dxdb, from local coordinate system to
        //- cartesian
        virtual tensor transformationTensorDxDb(label globalPointIndex) = 0;

        //- Update coordinates in the local system based on the cartesian points
        virtual void updateLocalCoordinateSystem
        (
            const vectorField& cartesianPoints
        ) = 0;


public:

    //- Runtime type information
    TypeName("NURBS3DVolume");


    // Declare run-time constructor selection table

        declareRunTimeSelectionTable
        (
            autoPtr,
            NURBS3DVolume,
            dictionary,
            (
                const dictionary& dict,
                const fvMesh& mesh,
                bool computeParamCoors
            ),
            (dict, mesh, computeParamCoors)
        );

    // Constructors

        //- Construct from dictionary
        NURBS3DVolume
        (
            const dictionary& dict,
            const fvMesh& mesh,
            bool computeParamCoors = true
        );

        //- Construct as copy
        NURBS3DVolume(const NURBS3DVolume&);


    // Selectors

        //- Return a reference to the selected NURBS model
        static autoPtr<NURBS3DVolume> New
        (
            const dictionary& dict,
            const fvMesh& mesh,
            bool computeParamCoors = true
        );


    //- Destructor
    virtual ~NURBS3DVolume() = default;


    // Member Functions

        // Derivatives wrt parametric coordinates

            //- Volume point derivative wrt u at point u,v,w
            vector volumeDerivativeU
            (
                const scalar u,
                const scalar v,
                const scalar w
            ) const;

            //- Volume point derivative wrt v at point u,v,w
            vector volumeDerivativeV
            (
                const scalar u,
                const scalar v,
                const scalar w
            ) const;

            //- Volume point derivative wrt w at point u,v,w
            vector volumeDerivativeW
            (
                const scalar u,
                const scalar v,
                const scalar w
            ) const;

            //- Jacobian matrix wrt to the volume parametric coordinates
            tensor JacobianUVW(const vector& u) const;


        // Derivatives wrt to control points

            //- Volume point derivative wrt to control point cpI at point u,v,w
            //  Scalar since in the local system!
            scalar volumeDerivativeCP(const vector& u, const label cpI) const;

            //- Control point sensitivities computed using point-based surface
            //- sensitivities
            vectorField computeControlPointSensitivities
            (
                const pointVectorField& pointSens,
                const labelList& sensitivityPatchIDs
            );

            //- Control point sensitivities computed using face-based surface
            //- sensitivities
            vectorField computeControlPointSensitivities
            (
                const volVectorField& faceSens,
                const labelList& sensitivityPatchIDs
            );

            //- Control point sensitivities computed using face-based surface
            //- sensitivities
            vectorField computeControlPointSensitivities
            (
                const boundaryVectorField& faceSens,
                const labelList& sensitivityPatchIDs
            );

            //- Control point sensitivities computed using face-based surface
            //- sensitivities
            vector computeControlPointSensitivities
            (
                const vectorField& faceSens,
                const label patchI,
                const label cpI
            );

            //- Part of control point sensitivities related to the face normal
            //- variations
            tmp<tensorField> dndbBasedSensitivities
            (
                const label patchI,
                const label cpI,
                bool DimensionedNormalSens = true
            );

            //- Get patch dx/db
            tmp<tensorField> patchDxDb
            (
                const label patchI,
                const label cpI
            );

            //- Get patch dx/db
            tmp<tensorField> patchDxDbFace
            (
                const label patchI,
                const label cpI
            );


        // Cartesian coordinates

            //- Compute cartesian coordinates based on control points
            //- and parametric coordinates
            tmp<vectorField> coordinates(const vectorField& uVector) const;

            //- The same, for a specific point
            vector coordinates(const vector& uVector) const;

            //- Mesh movement based on given control point movement
            tmp<vectorField> computeNewPoints
            (
                const vectorField& controlPointsMovement
            );

            //- Boundary mesh movement based on given control point movement
            tmp<vectorField> computeNewBoundaryPoints
            (
                const vectorField& controlPointsMovement,
                const labelList& patchesToBeMoved,
                const bool moveCPs = true
            );


        // Control points management

            //- Get control point ID from its I-J-K coordinates
            label getCPID(const label i, const label j, const label k) const;

            //- Set new control points
            //  New values should be on the coordinates system original CPs
            //  were defined
            void setControlPoints(const vectorField& newCps);


            //- Bound control points movement in the boundary control points
            //- and in certain directions if needed
            void boundControlPointMovement
            (
                vectorField& controlPointsMovement
            ) const;

            //- Compute max. displacement at the boundary
            scalar computeMaxBoundaryDisplacement
            (
                const vectorField& controlPointsMovement,
                const labelList& patchesToBeMoved
            );


        // Access Functions

            // Functions triggering calculations

                //- Get mesh points that reside within the control points box
                tmp<vectorField> getPointsInBox();

                //- Get map of points in box to mesh points
                const labelList& getMap();

                //- Get map of mesh points to points in box.
                //- Return -1 if point is outside the box
                const labelList& getReverseMap();

                //- Get parametric coordinates
                const pointVectorField& getParametricCoordinates();

                //- Get dxCartesiandb for a certain control point
                tmp<pointTensorField> getDxDb(const label cpI);

                //- Get dxCartesiandb for a certain control point on cells
                tmp<volTensorField> getDxCellsDb(const label cpI);

                //- Get number of variables if CPs are moved symmetrically in U
                label nUSymmetry() const;

                //- Get number of variables if CPs are moved symmetrically in V
                label nVSymmetry() const;

                //- 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

                //- Get box name
                inline const word& name() const;

                //- Which control points are active?
                //  A control point is active if at least one component can move
                inline const boolList& getActiveCPs() const;

                //- Which design variables are active?
                //  Numbering is (X1,Y1,Z1), (X2,Y2,Z2) ...
                inline const boolList& getActiveDesignVariables() const;

                //- Get control points
                inline const vectorField& getControlPoints() const;

                inline vectorField& getControlPoints();

                //- Get confine movements
                inline bool confineUMovement() const;

                inline bool confineVMovement() const;

                inline bool confineWMovement() const;

                //- Get basis functions
                inline const NURBSbasis& basisU() const;

                inline const NURBSbasis& basisV() const;

                inline const NURBSbasis& basisW() const;

                //- Get number of control points per direction
                inline Vector<label> nCPsPerDirection() const;

                //- Get mesh
                inline const fvMesh& mesh() const;

                //- Get dictionary
                inline const dictionary& dict() const;


        // Write Functions

            //- Write control points on a cartesian coordinates system for
            //- visualization
            void writeCps
            (
                const fileName& = "cpsFile",
                const bool transform = true
            ) const;

            //- Write parametric coordinates
            void writeParamCoordinates() const;

            //- Write the control points to support restart
            virtual bool writeData(Ostream& os) const;
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#include "NURBS3DVolumeI.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //