NURBS3DVolume.C

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-2019 PCOpt/NTUA
    Copyright (C) 2013-2019 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/>.
 
\*---------------------------------------------------------------------------*/
 
#include "addToRunTimeSelectionTable.H"
#include "NURBS3DVolume.H"
 
#include "OFstream.H"
#include "Time.H"
#include "deltaBoundary.H"
#include "coupledFvPatch.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
    defineTypeNameAndDebug(NURBS3DVolume, 0);
    defineRunTimeSelectionTable(NURBS3DVolume, dictionary);
}
 
 
// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
Foam::label Foam::NURBS3DVolume::getCPID
(
    const label i,
    const label j,
    const label k
) const
{
    const label nCPsU = basisU_.nCPs();
    const label nCPsV = basisV_.nCPs();
 
    return k*nCPsU*nCPsV + j*nCPsU + i;
}
 
 
void Foam::NURBS3DVolume::findPointsInBox(const vectorField& meshPoints)
{
    // It is considered an error to recompute points in the control boxes
    if (mapPtr_.valid() || reverseMapPtr_.valid())
    {
        FatalErrorInFunction
           << "Attempting to recompute points residing within control boxes"
           << exit(FatalError);
    }
 
    mapPtr_.reset(new labelList(meshPoints.size(), -1));
    reverseMapPtr_.reset(new labelList(meshPoints.size(), -1));
    labelList& map = mapPtr_();
    labelList& reverseMap = reverseMapPtr_();
 
    // Identify points inside morphing boxes
    scalar lowerX = min(cps_.component(0));
    scalar upperX = max(cps_.component(0));
    scalar lowerY = min(cps_.component(1));
    scalar upperY = max(cps_.component(1));
    scalar lowerZ = min(cps_.component(2));
    scalar upperZ = max(cps_.component(2));
 
    Info<< "Control Points bounds \n"
        << "\tX1 : (" << lowerX << " " << upperX <<  ")\n"
        << "\tX2 : (" << lowerY << " " << upperY <<  ")\n"
        << "\tX3 : (" << lowerZ << " " << upperZ <<  ")\n" << endl;
 
    label count(0);
    forAll(meshPoints, pI)
    {
        const vector& pointI = meshPoints[pI];
        if
        (
            pointI.x() >= lowerX && pointI.x() <= upperX
         && pointI.y() >= lowerY && pointI.y() <= upperY
         && pointI.z() >= lowerZ && pointI.z() <= upperZ
        )
        {
            map[count] = pI;
            reverseMap[pI] = count;
            ++count;
        }
    }
 
    // Resize lists
    map.setSize(count);
 
    reduce(count, sumOp<label>());
    Info<< "Initially found " << count << " points inside control boxes"
        << endl;
}
 
 
void Foam::NURBS3DVolume::computeParametricCoordinates
(
    const vectorField& points
)
{
    scalar timeBef = mesh_.time().elapsedCpuTime();
 
    if (parametricCoordinatesPtr_.valid())
    {
        FatalErrorInFunction
           << "Attempting to recompute parametric coordinates"
           << exit(FatalError);
    }
 
    labelList& map = mapPtr_();
    labelList& reverseMap = reverseMapPtr_();
 
    parametricCoordinatesPtr_.reset
    (
        new pointVectorField
        (
            IOobject
            (
                "parametricCoordinates" + name_,
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            pointMesh::New(mesh_),
            dimensionedVector(dimless, Zero)
        )
    );
    vectorField& paramCoors = parametricCoordinatesPtr_().primitiveFieldRef();
 
    // If already present, read from file
    if
    (
        parametricCoordinatesPtr_().typeHeaderOk<pointVectorField>(true)
     && readStoredData_
    )
    {
        // These are the parametric coordinates that have been computed
        // correctly (hopefully!) in a previous run.
        // findPointsInBox has possibly found more points and lists need
        // resizing
        Info<< "Reading parametric coordinates from file" << endl;
        IOobject& header = parametricCoordinatesPtr_().ref();
        parametricCoordinatesPtr_() =
            pointVectorField
            (
                header,
                pointMesh::New(mesh_)
            );
 
        // Initialize intermediate fields with sizes from findPointInBox
        labelList actualMap(map.size());
 
        // Read and store
        label curIndex(0);
        forAll(map, pI)
        {
            const label globalPointIndex = map[pI];
            if (paramCoors[globalPointIndex] != vector::zero)
            {
                actualMap[curIndex] = map[pI];
                reverseMap[globalPointIndex] = curIndex;
                ++curIndex;
            }
            else
            {
                reverseMap[globalPointIndex] = -1;
            }
        }
 
        // Resize intermediate
        actualMap.setSize(curIndex);
 
        reduce(curIndex, sumOp<label>());
        Info<< "Read non-zero parametric coordinates for " << curIndex
            << " points" << endl;
 
        // Update lists with the appropriate entries
        map = actualMap;
    }
    // Else, compute parametric coordinates iteratively
    else
    {
        // Initialize parametric coordinates based on min/max of control points
        scalar minX1 = min(cps_.component(0));
        scalar maxX1 = max(cps_.component(0));
        scalar minX2 = min(cps_.component(1));
        scalar maxX2 = max(cps_.component(1));
        scalar minX3 = min(cps_.component(2));
        scalar maxX3 = max(cps_.component(2));
 
        scalar oneOverDenomX(1./(maxX1 - minX1));
        scalar oneOverDenomY(1./(maxX2 - minX2));
        scalar oneOverDenomZ(1./(maxX3 - minX3));
 
        forAll(map, pI)
        {
            const label globalPI = map[pI];
            paramCoors[globalPI].x() = (points[pI].x() - minX1)*oneOverDenomX;
            paramCoors[globalPI].y() = (points[pI].y() - minX2)*oneOverDenomY;
            paramCoors[globalPI].z() = (points[pI].z() - minX3)*oneOverDenomZ;
        }
 
        // Indices of points that failed to converge
        // (i.e. are bounded for nMaxBound iters)
        boolList dropOffPoints(map.size(), false);
        label nDropedPoints(0);
 
        // Initial cartesian coordinates
        tmp<vectorField> tsplinesBasedCoors(coordinates(paramCoors));
        vectorField& splinesBasedCoors = tsplinesBasedCoors.ref();
 
        // Newton-Raphson loop to compute parametric coordinates
        // based on cartesian coordinates and the known control points
        Info<< "Mapping of mesh points to parametric space for box " << name_
            << " ..." << endl;
        // Do loop on a point-basis and check residual of each point equation
        label maxIterNeeded(0);
        forAll(points, pI)
        {
            label iter(0);
            label nBoundIters(0);
            vector res(GREAT, GREAT, GREAT);
            do
            {
                const label globalPI = map[pI];
                vector& uVec = paramCoors[globalPI];
                vector& coorPointI = splinesBasedCoors[pI];
                uVec += ((inv(JacobianUVW(uVec))) & (points[pI] - coorPointI));
                // Bounding might be needed for the first iterations
                // If multiple bounds happen, point is outside of the control
                // boxes and should be discarded
                if (bound(uVec))
                {
                    ++nBoundIters;
                }
                if (nBoundIters > nMaxBound_)
                {
                    dropOffPoints[pI] = true;
                    ++nDropedPoints;
                    break;
                }
                // Update current cartesian coordinates based on parametric ones
                coorPointI = coordinates(uVec);
                // Compute residual
                res = cmptMag(points[pI] - coorPointI);
            }
            while
            (
                (iter++ < maxIter_)
             && (
                       res.component(0) > tolerance_
                    || res.component(1) > tolerance_
                    || res.component(2) > tolerance_
                )
            );
            if (iter > maxIter_)
            {
                WarningInFunction
                    << "Mapping to parametric space for point " << pI
                    << " failed." << endl
                    << "Residual after " << maxIter_ + 1 << " iterations : "
                    << res << endl
                    << "parametric coordinates " <<  paramCoors[map[pI]]
                    <<  endl
                    << "Local system coordinates " <<  points[pI] <<  endl
                    << "Threshold residual per direction : " << tolerance_
                    << endl;
            }
            maxIterNeeded = max(maxIterNeeded, iter);
        }
        reduce(maxIterNeeded, maxOp<label>());
 
        label nParameterizedPoints = map.size() - nDropedPoints;
 
        // Resize mapping lists and parametric coordinates after dropping some
        // points
        labelList mapOld(map);
 
        map.setSize(nParameterizedPoints);
 
        label curIndex(0);
        forAll(dropOffPoints, pI)
        {
            if (!dropOffPoints[pI])
            {
                map[curIndex] = mapOld[pI];
                reverseMap[mapOld[pI]] = curIndex;
                ++curIndex;
            }
            else
            {
                paramCoors[mapOld[pI]] = vector::zero;
                reverseMap[mapOld[pI]] = -1;
            }
        }
 
        reduce(nDropedPoints, sumOp<label>());
        reduce(nParameterizedPoints, sumOp<label>());
        Info<< "Found " << nDropedPoints
            << " to discard from morphing boxes" << endl;
        Info<< "Keeping " << nParameterizedPoints
            << " parameterized points in boxes"  << endl;
 
        splinesBasedCoors = coordinates(paramCoors)();
        scalar maxDiff(-GREAT);
        forAll(splinesBasedCoors, pI)
        {
            scalar diff =
                mag(splinesBasedCoors[pI] - localSystemCoordinates_[map[pI]]);
            if (diff > maxDiff)
            {
                maxDiff = diff;
            }
        }
        reduce(maxDiff, maxOp<scalar>());
        scalar timeAft = mesh_.time().elapsedCpuTime();
        Info<< "\tMapping completed in " << timeAft - timeBef << " seconds"
            << endl;
        Info<< "\tMax iterations per point needed to compute parametric "
            << "coordinates : "
            << maxIterNeeded << endl;
        Info<< "\tMax difference between original mesh points and "
            << "parameterized ones "
            << maxDiff << endl;
    }
}
 
 
void Foam::NURBS3DVolume::computeParametricCoordinates
(
    tmp<vectorField> tPoints
)
{
    const vectorField& points = tPoints();
    computeParametricCoordinates(points);
}
 
 
bool Foam::NURBS3DVolume::bound
(
    vector& vec,
    scalar minValue,
    scalar maxValue
)
{
    bool boundPoint(false);
    // Lower value bounding
    if (vec.x() < scalar(0))
    {
        vec.x() = minValue;
        boundPoint = true;
    }
    if (vec.y() < scalar(0))
    {
        vec.y() = minValue;
        boundPoint = true;
    }
    if (vec.z() < scalar(0))
    {
        vec.z() = minValue;
        boundPoint = true;
    }
    // Upper value bounding
    if (vec.x() > 1)
    {
        vec.x() = maxValue;
        boundPoint = true;
    }
    if (vec.y() > 1)
    {
         vec.y() = maxValue;
        boundPoint = true;
    }
    if (vec.z() > 1)
    {
        vec.z() = maxValue;
        boundPoint = true;
    }
    return boundPoint;
}
 
 
void Foam::NURBS3DVolume::makeFolders()
{
    if (Pstream::master())
    {
        mkDir(mesh_.time().globalPath()/"optimisation"/cpsFolder_);
    }
}
 
 
void Foam::NURBS3DVolume::determineActiveDesignVariablesAndPoints()
{
    label nCPs = cps_.size();
    activeControlPoints_ = boolList(nCPs, true);
    activeDesignVariables_ = boolList(3*nCPs, true);
 
    // Check whether all boundary control points should be confined
    confineBoundaryControlPoints();
 
    // Apply confinement to maintain continuity
    continuityRealatedConfinement();
 
    // Confine user-specified directions
    confineControlPointsDirections();
 
    // Determine active control points. A control point is considered active
    // if at least one of its components is free to move
    forAll(activeControlPoints_, cpI)
    {
        if
        (
            !activeDesignVariables_[3*cpI]
         && !activeDesignVariables_[3*cpI + 1]
         && !activeDesignVariables_[3*cpI + 2]
        )
        {
            activeControlPoints_[cpI] = false;
        }
    }
}
 
 
void Foam::NURBS3DVolume::confineBoundaryControlPoints()
{
    const label nCPsU = basisU_.nCPs();
    const label nCPsV = basisV_.nCPs();
    const label nCPsW = basisW_.nCPs();
 
    // Zero movement of the boundary control points. Active by default
    if (confineBoundaryControlPoints_)
    {
        // Side patches
        for (label iCPw = 0; iCPw < nCPsW; iCPw += nCPsW - 1)
        {
            for (label iCPv = 0; iCPv < nCPsV; iCPv++)
            {
                for (label iCPu = 0; iCPu < nCPsU; iCPu++)
                {
                    confineControlPoint(getCPID(iCPu, iCPv, iCPw));
                }
            }
        }
        // Front-back patches
        for (label iCPw = 0; iCPw < nCPsW; iCPw++)
        {
            for (label iCPv = 0; iCPv < nCPsV; iCPv++)
            {
                for (label iCPu = 0; iCPu < nCPsU; iCPu += nCPsU - 1)
                {
                    confineControlPoint(getCPID(iCPu, iCPv, iCPw));
                }
            }
        }
        // Top-bottom patches
        for (label iCPw = 0; iCPw < nCPsW; iCPw++)
        {
            for (label iCPv = 0; iCPv < nCPsV; iCPv += nCPsV - 1)
            {
                for (label iCPu = 0; iCPu < nCPsU; iCPu++)
                {
                    confineControlPoint(getCPID(iCPu, iCPv, iCPw));
                }
            }
        }
    }
}
 
 
void Foam::NURBS3DVolume::continuityRealatedConfinement()
{
    const label nCPsU = basisU_.nCPs();
    const label nCPsV = basisV_.nCPs();
    const label nCPsW = basisW_.nCPs();
 
    // Zero movement to a number of x-constant slices of cps in order to
    // preserve continuity at the boundary of the parameterized space
    forAll(confineUMinCPs_, iCPu)
    {
        const FixedList<bool, 3>& confineSlice = confineUMinCPs_[iCPu];
        // Control points at the start of the parameterized space
        for (label iCPw = 0; iCPw < nCPsW; iCPw++)
        {
            for (label iCPv = 0; iCPv < nCPsV; iCPv++)
            {
                confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
            }
        }
    }
 
    forAll(confineUMaxCPs_, sliceI)
    {
        const FixedList<bool, 3>& confineSlice = confineUMaxCPs_[sliceI];
        label iCPu = nCPsU - 1 - sliceI;
        // Control points at the end of the parameterized space
        for (label iCPw = 0; iCPw < nCPsW; iCPw++)
        {
            for (label iCPv = 0; iCPv < nCPsV; iCPv++)
            {
                confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
            }
        }
    }
 
    // Zero movement to a number of y-constant slices of cps in order to
    // preserve continuity at the boundary of the parameterized space
    forAll(confineVMinCPs_, iCPv)
    {
        const FixedList<bool, 3>& confineSlice = confineVMinCPs_[iCPv];
        // Control points at the start of the parameterized space
        for (label iCPw = 0; iCPw < nCPsW; iCPw++)
        {
            for (label iCPu = 0; iCPu < nCPsU; iCPu++)
            {
                confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
            }
        }
    }
 
    forAll(confineVMaxCPs_, sliceI)
    {
        const FixedList<bool, 3>& confineSlice = confineVMaxCPs_[sliceI];
        label iCPv = nCPsV - 1 - sliceI;
        // Control points at the end of the parameterized space
        for (label iCPw = 0; iCPw < nCPsW; iCPw++)
        {
            for (label iCPu = 0; iCPu < nCPsU; iCPu++)
            {
                confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
            }
        }
    }
 
    // Zero movement to a number of w-constant slices of cps in order to
    // preserve continuity at the boundary of the parameterized space
    forAll(confineWMinCPs_, iCPw)
    {
        const FixedList<bool, 3>& confineSlice = confineWMinCPs_[iCPw];
        // Control points at the start of the parameterized space
        for (label iCPv = 0; iCPv < nCPsV; iCPv++)
        {
            for (label iCPu = 0; iCPu < nCPsU; iCPu++)
            {
                confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
            }
        }
    }
 
    forAll(confineWMaxCPs_, sliceI)
    {
        const FixedList<bool, 3>& confineSlice = confineWMaxCPs_[sliceI];
        label iCPw = nCPsW - 1 - sliceI;
        // Control points at the end of the parameterized space
        for (label iCPv = 0; iCPv < nCPsV; iCPv++)
        {
            for (label iCPu = 0; iCPu < nCPsU; iCPu++)
            {
                confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
            }
        }
    }
}
 
 
void Foam::NURBS3DVolume::confineControlPointsDirections()
{
    for (label cpI = 0; cpI < cps_.size(); ++cpI)
    {
        if (confineUMovement_) activeDesignVariables_[3*cpI] = false;
        if (confineVMovement_) activeDesignVariables_[3*cpI + 1] = false;
        if (confineWMovement_) activeDesignVariables_[3*cpI + 2] = false;
    }
}
 
 
void Foam::NURBS3DVolume::confineControlPoint(const label cpI)
{
    if (cpI < 0 || cpI > cps_.size() -1)
    {
        FatalErrorInFunction
           << "Attempted to confine control point movement for a control point "
           << " ID which is out of bounds"
           << exit(FatalError);
    }
    else
    {
        activeDesignVariables_[3*cpI] = false;
        activeDesignVariables_[3*cpI + 1] = false;
        activeDesignVariables_[3*cpI + 2] = false;
    }
}
 
 
void Foam::NURBS3DVolume::confineControlPoint
(
    const label cpI,
    const FixedList<bool, 3>& confineDirections
)
{
    if (cpI < 0 || cpI > cps_.size() -1)
    {
        FatalErrorInFunction
           << "Attempted to confine control point movement for a control point "
           << " ID which is out of bounds"
           << exit(FatalError);
    }
    else
    {
        if (confineDirections[0]) activeDesignVariables_[3*cpI] = false;
        if (confineDirections[1]) activeDesignVariables_[3*cpI + 1] = false;
        if (confineDirections[2]) activeDesignVariables_[3*cpI + 2] = false;
    }
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::NURBS3DVolume::NURBS3DVolume
(
    const dictionary& dict,
    const fvMesh& mesh,
    bool computeParamCoors
)
:
    mesh_(mesh),
    name_(dict.dictName()),
    basisU_(dict.get<label>("nCPsU"), dict.get<label>("degreeU")),
    basisV_(dict.get<label>("nCPsV"), dict.get<label>("degreeV")),
    basisW_(dict.get<label>("nCPsW"), dict.get<label>("degreeW")),
    maxIter_(dict.getOrDefault<label>("maxIterations", 10)),
    tolerance_(dict.getOrDefault<scalar>("tolerance", 1.e-10)),
    nMaxBound_(dict.getOrDefault<scalar>("nMaxBoundIterations", 4)),
    cps_(0),
    mapPtr_(nullptr),
    reverseMapPtr_(nullptr),
    parametricCoordinatesPtr_(nullptr),
    localSystemCoordinates_(mesh_.nPoints(), Zero),
    confineUMovement_
    (
        dict.getOrDefaultCompat<bool>
        (
            "confineUMovement", {{"confineX1movement", 1912}}, false
        )
    ),
    confineVMovement_
    (
        dict.getOrDefaultCompat<bool>
        (
            "confineVMovement", {{"confineX2movement", 1912}}, false
        )
    ),
    confineWMovement_
    (
        dict.getOrDefaultCompat<bool>
        (
            "confineWMovement", {{"confineX3movement", 1912}}, false
        )
    ),
    confineBoundaryControlPoints_
    (
        dict.getOrDefault<bool>("confineBoundaryControlPoints", true)
    ),
    confineUMinCPs_
    (
        dict.getOrDefaultCompat<boolListList3>
        (
            "confineUMinCPs", {{"boundUMinCPs", 1912}}, boolListList3(0)
        )
    ),
    confineUMaxCPs_
    (
        dict.getOrDefaultCompat<boolListList3>
        (
            "confineUMaxCPs", {{"boundUMaxCPs", 1912}}, boolListList3(0)
        )
    ),
    confineVMinCPs_
    (
        dict.getOrDefaultCompat<boolListList3>
        (
            "confineVMinCPs", {{"boundVMinCPs", 1912}}, boolListList3(0)
        )
    ),
    confineVMaxCPs_
    (
        dict.getOrDefaultCompat<boolListList3>
        (
            "confineVMaxCPs", {{"boundVMaxCPs", 1912}}, boolListList3(0)
        )
    ),
    confineWMinCPs_
    (
        dict.getOrDefaultCompat<boolListList3>
        (
            "confineWMinCPs", {{"boundWMinCPs", 1912}}, boolListList3(0)
        )
    ),
    confineWMaxCPs_
    (
        dict.getOrDefaultCompat<boolListList3>
        (
            "confineWMaxCPs", {{"boundWMaxCPs", 1912}}, boolListList3(0)
        )
    ),
    activeControlPoints_(0), //zero here, execute sanity checks first
    activeDesignVariables_(0), //zero here, execute sanity checks first
    cpsFolder_("controlPoints"),
    readStoredData_(dict.getOrDefault<bool>("readStoredData", true))
{
    // Create folders
    makeFolders();
 
    // Sanity checks
    if
    (
        (confineUMinCPs_.size() + confineUMaxCPs_.size() >= basisU_.nCPs())
     || (confineVMinCPs_.size() + confineVMaxCPs_.size() >= basisV_.nCPs())
     || (confineWMinCPs_.size() + confineWMaxCPs_.size() >= basisW_.nCPs())
    )
    {
        FatalErrorInFunction
           << "Number of control point slices to be kept frozen at "
           << "the boundaries is invalid \n"
           << "Number of control points in u " << basisU_.nCPs() << "\n"
           << "Number of control points in v " << basisV_.nCPs() << "\n"
           << "Number of control points in w " << basisW_.nCPs() << "\n"
           << exit(FatalError);
    }
 
    word controlPointsDefinition(dict.get<word>("controlPointsDefinition"));
 
    // Read specific control points if given
    if (controlPointsDefinition == "fromFile")
    {
        Info<< "Reading control points from file " << endl;
        IOdictionary cpsDict
        (
            IOobject
            (
                name_ + "cpsBsplines" + mesh_.time().timeName(),
                mesh_.time().caseConstant(),
                cpsFolder_,
                mesh_,
                IOobject::MUST_READ,
                IOobject::NO_WRITE,
                false
            )
        );
 
        cpsDict.readEntry("controlPoints", cps_);
        if (cps_.size() != (basisU_.nCPs()*basisV_.nCPs()*basisW_.nCPs()))
        {
            FatalErrorInFunction
               << "Number of control points does not agree with "
               << "nCPsU*nCPv*nCPsW"
               << exit(FatalError);
        }
    }
    // Else, construct control points based on their min, max coordinates and
    // their number
    else if (controlPointsDefinition == "axisAligned")
    {
        const label nCPsU = basisU_.nCPs();
        const label nCPsV = basisV_.nCPs();
        const label nCPsW = basisW_.nCPs();
 
        cps_.setSize(nCPsU * nCPsV * nCPsW);
 
        vector lowerBounds(dict.get<vector>("lowerCpBounds"));
        vector upperBounds(dict.get<vector>("upperCpBounds"));
 
        // Equidistribute cps in th u,v,w directions
        for (label iCPw = 0; iCPw < nCPsW; ++iCPw)
        {
            for (label iCPv = 0; iCPv < nCPsV; ++iCPv)
            {
                for (label iCPu = 0; iCPu < nCPsU; ++iCPu)
                {
                    cps_[getCPID(iCPu, iCPv, iCPw)] =
                        (
                            lowerBounds.x()
                           +scalar(iCPu)/scalar(nCPsU - 1)
                           *(upperBounds.x()-lowerBounds.x())
                        )*vector(1, 0, 0)
                      + (
                            lowerBounds.y()
                           +scalar(iCPv)/scalar(nCPsV - 1)
                           *(upperBounds.y()-lowerBounds.y())
                        )*vector(0, 1, 0)
                      + (
                            lowerBounds.z()
                           +scalar(iCPw)/scalar(nCPsW - 1)
                           *(upperBounds.z()-lowerBounds.z())
                        )*vector(0, 0, 1);
                }
            }
        }
    }
    else
    {
        FatalErrorInFunction
            << "Unknown controlPointsDefinition type "
            << controlPointsDefinition
            << endl << endl
            << "Valid types are : axisAligned, fromFile" << endl
            << exit(FatalError);
    }
    determineActiveDesignVariablesAndPoints();
    writeCpsInDict();
}
 
 
// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
 
Foam::autoPtr<Foam::NURBS3DVolume> Foam::NURBS3DVolume::New
(
    const dictionary& dict,
    const fvMesh& mesh,
    bool computeParamCoors
)
{
    const word modelType(dict.get<word>("type"));
 
    Info<< "NURBS3DVolume type : " << modelType << endl;
 
    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
 
    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
            dict,
            "type",
            modelType,
            *dictionaryConstructorTablePtr_
        ) << exit(FatalIOError);
    }
 
    return autoPtr<NURBS3DVolume>(cstrIter()(dict, mesh, computeParamCoors));
}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
Foam::vector Foam::NURBS3DVolume::volumeDerivativeU
(
    const scalar u,
    const scalar v,
    const scalar w
) const
{
    const label degreeU = basisU_.degree();
    const label degreeV = basisV_.degree();
    const label degreeW = basisW_.degree();
 
    const label nCPsU = basisU_.nCPs();
    const label nCPsV = basisV_.nCPs();
    const label nCPsW = basisW_.nCPs();
 
    vector derivative(Zero);
 
    for (label iCPw = 0; iCPw < nCPsW; ++iCPw)
    {
        for (label iCPv = 0; iCPv < nCPsV; ++iCPv)
        {
            for (label iCPu = 0; iCPu < nCPsU; ++iCPu)
            {
                derivative +=
                    cps_[getCPID(iCPu, iCPv, iCPw)]
                   *basisU_.basisDerivativeU(iCPu, degreeU, u)
                   *basisV_.basisValue(iCPv, degreeV, v)
                   *basisW_.basisValue(iCPw, degreeW, w);
            }
        }
    }
 
    return derivative;
}
 
 
Foam::vector Foam::NURBS3DVolume::volumeDerivativeV
(
    const scalar u,
    const scalar v,
    const scalar w
) const
{
    const label degreeU = basisU_.degree();
    const label degreeV = basisV_.degree();
    const label degreeW = basisW_.degree();
 
    const label nCPsU = basisU_.nCPs();
    const label nCPsV = basisV_.nCPs();
    const label nCPsW = basisW_.nCPs();
 
    vector derivative(Zero);
 
    for (label iCPw = 0; iCPw < nCPsW; ++iCPw)
    {
        for (label iCPv = 0; iCPv < nCPsV; ++iCPv)
        {
            for (label iCPu = 0; iCPu < nCPsU; ++iCPu)
            {
                derivative +=
                    cps_[getCPID(iCPu, iCPv, iCPw)]
                   *basisU_.basisValue(iCPu, degreeU, u)
                   *basisV_.basisDerivativeU(iCPv, degreeV, v)
                   *basisW_.basisValue(iCPw, degreeW, w);
            }
        }
    }
 
    return derivative;
}
 
 
Foam::vector Foam::NURBS3DVolume::volumeDerivativeW
(
    const scalar u,
    const scalar v,
    const scalar w
) const
{
    const label degreeU = basisU_.degree();
    const label degreeV = basisV_.degree();
    const label degreeW = basisW_.degree();
 
    const label nCPsU = basisU_.nCPs();
    const label nCPsV = basisV_.nCPs();
    const label nCPsW = basisW_.nCPs();
 
    vector derivative(Zero);
 
    for (label iCPw = 0; iCPw < nCPsW; iCPw++)
    {
        for (label iCPv = 0; iCPv < nCPsV; iCPv++)
        {
            for (label iCPu = 0; iCPu < nCPsU; iCPu++)
            {
                derivative +=
                    cps_[getCPID(iCPu, iCPv, iCPw)]
                   *basisU_.basisValue(iCPu, degreeU, u)
                   *basisV_.basisValue(iCPv, degreeV, v)
                   *basisW_.basisDerivativeU(iCPw, degreeW, w);
            }
        }
    }
 
    return derivative;
}
 
 
Foam::tensor Foam::NURBS3DVolume::JacobianUVW
(
    const vector& uVector
) const
{
    const scalar u = uVector.x();
    const scalar v = uVector.y();
    const scalar w = uVector.z();
 
    vector uDeriv = volumeDerivativeU(u, v, w);
    vector vDeriv = volumeDerivativeV(u, v, w);
    vector wDeriv = volumeDerivativeW(u, v, w);
 
    tensor Jacobian(Zero);
 
    Jacobian[0] = uDeriv.component(0);
    Jacobian[1] = vDeriv.component(0);
    Jacobian[2] = wDeriv.component(0);
    Jacobian[3] = uDeriv.component(1);
    Jacobian[4] = vDeriv.component(1);
    Jacobian[5] = wDeriv.component(1);
    Jacobian[6] = uDeriv.component(2);
    Jacobian[7] = vDeriv.component(2);
    Jacobian[8] = wDeriv.component(2);
 
    return Jacobian;
}
 
 
Foam::scalar Foam::NURBS3DVolume::volumeDerivativeCP
(
    const vector& uVector,
    const label cpI
) const
{
    const scalar u = uVector.x();
    const scalar v = uVector.y();
    const scalar w = uVector.z();
 
    const label nCPsU = basisU_.nCPs();
    const label nCPsV = basisV_.nCPs();
 
    const label degreeU = basisU_.degree();
    const label degreeV = basisV_.degree();
    const label degreeW = basisW_.degree();
 
    label iCPw = cpI/label(nCPsU*nCPsV);
    label iCPv = (cpI - iCPw*nCPsU*nCPsV)/nCPsU;
    label iCPu = (cpI - iCPw*nCPsU*nCPsV - iCPv*nCPsU);
 
    // Normally, this should be a tensor, however the parameterization is
    // isotropic. Hence the tensor degenerates to a diagonal tensor with all
    // diagonal elements being equal. This returns the (unique) diag element
    scalar derivative =
        basisU_.basisValue(iCPu, degreeU, u)
       *basisV_.basisValue(iCPv, degreeV, v)
       *basisW_.basisValue(iCPw, degreeW, w);
 
    return derivative;
}
 
 
Foam::vectorField Foam::NURBS3DVolume::computeControlPointSensitivities
(
    const pointVectorField& pointSens,
    const labelList& sensitivityPatchIDs
)
{
    vectorField controlPointDerivs(cps_.size(), Zero);
 
    // Get parametric coordinates
    const vectorField& parametricCoordinates = getParametricCoordinates();
 
    forAll(controlPointDerivs, cpI)
    {
        forAll(sensitivityPatchIDs, pI)
        {
            const label patchI = sensitivityPatchIDs[pI];
            const polyPatch& patch = mesh_.boundaryMesh()[patchI];
            const labelList& meshPoints = patch.meshPoints();
 
            forAll(meshPoints, mpI)
            {
                const label globalIndex = meshPoints[mpI];
                const label whichPointInBox = reverseMapPtr_()[globalIndex];
 
                // If point resides within control points box,
                // add contribution to cp derivative
                if (whichPointInBox != -1)
                {
                    controlPointDerivs[cpI] +=
                    (
                        pointSens[globalIndex]
                      & transformationTensorDxDb(globalIndex)
                    )
                   *volumeDerivativeCP
                    (
                        parametricCoordinates[globalIndex],
                        cpI
                    );
                }
            }
        }
    }
 
    // Sum contributions from all processors
    Pstream::listCombineGather(controlPointDerivs, plusEqOp<vector>());
    Pstream::listCombineScatter(controlPointDerivs);
 
    return controlPointDerivs;
}
 
 
Foam::vectorField Foam::NURBS3DVolume::computeControlPointSensitivities
(
    const volVectorField& faceSens,
    const labelList& sensitivityPatchIDs
)
{
    return
        computeControlPointSensitivities
        (
            faceSens.boundaryField(),
            sensitivityPatchIDs
        );
}
 
 
Foam::vectorField Foam::NURBS3DVolume::computeControlPointSensitivities
(
    const boundaryVectorField& faceSens,
    const labelList& sensitivityPatchIDs
)
{
    // Return field
    vectorField controlPointDerivs(cps_.size(), Zero);
 
    // Get parametric coordinates
    const vectorField& parametricCoordinates = getParametricCoordinates();
 
    // Auxiliary quantities
    deltaBoundary deltaBoundary(mesh_);
    const labelList& reverseMap = reverseMapPtr_();
 
    forAll(controlPointDerivs, cpI)
    {
        forAll(sensitivityPatchIDs, pI)
        {
            const label patchI = sensitivityPatchIDs[pI];
            const polyPatch& patch = mesh_.boundaryMesh()[patchI];
            const label patchStart = patch.start();
            const fvPatchVectorField& patchSens = faceSens[patchI];
 
            // loop over patch faces
            forAll(patch, fI)
            {
                const face& fGlobal = mesh_.faces()[fI + patchStart];
                const pointField facePoints = fGlobal.points(mesh_.points());
                // loop over face points
                tensorField facePointDerivs(facePoints.size(), Zero);
                forAll(fGlobal, pI)
                {
                    const label globalIndex = fGlobal[pI]; //global point index
                    const label whichPointInBox = reverseMap[globalIndex];
                    // if point resides within control points box,
                    // add contribution to d( facePoints )/db
                    if (whichPointInBox != -1)
                    {
                        // TENSOR-BASED
                        //~~~~~~~~~~~~~
                        facePointDerivs[pI] =
                            transformationTensorDxDb(globalIndex)
                          * volumeDerivativeCP
                            (
                                parametricCoordinates[globalIndex],
                                cpI
                            );
 
                    }
                }
 
                tensor fCtrs_d =
                    deltaBoundary.makeFaceCentresAndAreas_d
                    (
                        facePoints,
                        facePointDerivs
                    )[0];
                controlPointDerivs[cpI] += patchSens[fI] & fCtrs_d;
            }
        }
    }
    // Sum contributions from all processors
    Pstream::listCombineGather(controlPointDerivs, plusEqOp<vector>());
    Pstream::listCombineScatter(controlPointDerivs);
 
    return controlPointDerivs;
}
 
 
Foam::vector Foam::NURBS3DVolume::computeControlPointSensitivities
(
    const vectorField& faceSens,
    const label patchI,
    const label cpI
)
{
    // Return vector
    vector cpSens(Zero);
    // Get parametric coordinates
    const vectorField& parametricCoordinates = getParametricCoordinates();
 
    // Auxiliary quantities
    deltaBoundary deltaBoundary(mesh_);
    const labelList& reverseMap = reverseMapPtr_();
 
    const polyPatch& patch = mesh_.boundaryMesh()[patchI];
    const label patchStart = patch.start();
    // Loop over patch faces
    forAll(patch, fI)
    {
        const face& fGlobal = mesh_.faces()[fI + patchStart];
        const pointField facePoints = fGlobal.points(mesh_.points());
        // Loop over face points
        tensorField facePointDerivs(facePoints.size(), Zero);
        forAll(fGlobal, pI)
        {
            const label globalIndex = fGlobal[pI];  //global point index
            const label whichPointInBox = reverseMap[globalIndex];
            // If point resides within control points box,
            // add contribution to d( facePoints )/db
            if (whichPointInBox != -1)
            {
                // TENSOR-BASED
                //~~~~~~~~~~~~~
                facePointDerivs[pI] =
                    transformationTensorDxDb(globalIndex)
                   *volumeDerivativeCP
                    (
                        parametricCoordinates[globalIndex],
                        cpI
                    );
            }
        }
 
        tensor fCtrs_d =
            deltaBoundary.makeFaceCentresAndAreas_d
            (
                facePoints,
                facePointDerivs
            )[0];
        cpSens += faceSens[fI] & fCtrs_d;
    }
    // Sum contributions from all processors
    reduce(cpSens, sumOp<vector>());
 
    return cpSens;
}
 
 
Foam::tmp<Foam::tensorField> Foam::NURBS3DVolume::dndbBasedSensitivities
(
    const label patchI,
    const label cpI,
    bool DimensionedNormalSens
)
{
    const fvPatch& patch = mesh_.boundary()[patchI];
    const polyPatch& ppatch = patch.patch();
    // Return field
    tmp<tensorField> tdndbSens(new tensorField(patch.size(), Zero));
    tensorField& dndbSens = tdndbSens.ref();
    // Auxiliary quantities
    deltaBoundary deltaBoundary(mesh_);
    const label patchStart = ppatch.start();
    const labelList& reverseMap = reverseMapPtr_();
 
    // Get parametric coordinates
    const vectorField& parametricCoordinates = getParametricCoordinates();
 
    // Loop over patch faces
    forAll(patch, fI)
    {
        const face& fGlobal = mesh_.faces()[fI + patchStart];
        const pointField facePoints = fGlobal.points(mesh_.points());
        // Loop over face points
        tensorField facePointDerivs(facePoints.size(), Zero);
        forAll(fGlobal, pI)
        {
            const label globalIndex = fGlobal[pI];  //global point index
            const label whichPointInBox = reverseMap[globalIndex];
            // If point resides within control points box,
            // add contribution to d( facePoints )/db
            if (whichPointInBox != -1)
            {
                // TENSOR-BASED
                //~~~~~~~~~~~~~
                facePointDerivs[pI] =
                    transformationTensorDxDb(globalIndex)
                   *volumeDerivativeCP
                    (
                        parametricCoordinates[globalIndex],
                        cpI
                    );
            }
        }
 
        // Determine whether to return variance of dimensioned or unit normal
        tensorField dNdbSens =
            deltaBoundary.makeFaceCentresAndAreas_d
            (
                facePoints,
                facePointDerivs
            );
 
        if (DimensionedNormalSens)
        {
            dndbSens[fI] = dNdbSens[1];
        }
        else
        {
            dndbSens[fI] = dNdbSens[2];
        }
    }
 
    return tdndbSens;
}
 
 
Foam::tmp<Foam::tensorField> Foam::NURBS3DVolume::patchDxDb
(
    const label patchI,
    const label cpI
)
{
    // Get parametric coordinates
    const vectorField& parametricCoordinates = getParametricCoordinates();
 
    // Patch data
    const polyPatch& patch = mesh_.boundaryMesh()[patchI];
    const labelList& meshPoints = patch.meshPoints();
 
    // Return field
    auto tdxdb = tmp<tensorField>::New(patch.nPoints(), Zero);
    auto& dxdb = tdxdb.ref();
 
    forAll(meshPoints, pI)
    {
        const label globalIndex = meshPoints[pI];  //global point index
        const label whichPointInBox = reverseMapPtr_()[globalIndex];
 
        // If point resides within control points box, find dxdb
        if (whichPointInBox != -1)
        {
            dxdb[pI] =
                transformationTensorDxDb(globalIndex)
               *volumeDerivativeCP
                (
                    parametricCoordinates[globalIndex],
                    cpI
                );
        }
    }
 
    return tdxdb;
}
 
 
Foam::tmp<Foam::tensorField> Foam::NURBS3DVolume::patchDxDbFace
(
    const label patchI,
    const label cpI
)
{
    // get parametric coordinates
    const vectorField& parametricCoordinates = getParametricCoordinates();
 
    // Patch data
    const polyPatch& patch = mesh_.boundaryMesh()[patchI];
    const label patchStart = patch.start();
 
    // Return field
    auto tdxdb = tmp<tensorField>::New(patch.size(), Zero);
    auto& dxdb = tdxdb.ref();
 
    // Mesh differentiation engine
    deltaBoundary deltaBound(mesh_);
 
    forAll(patch, fI)
    {
        const face& fGlobal = mesh_.faces()[fI + patchStart];
        const pointField facePoints = fGlobal.points(mesh_.points());
        // Loop over face points
        tensorField facePointDerivs(facePoints.size(), Zero);
        forAll(fGlobal, pI)
        {
            const label globalIndex = fGlobal[pI];  //global point index
            const label whichPointInBox = reverseMapPtr_()[globalIndex];
            // If point resides within control points box,
            // add contribution to d( facePoints )/db
            if (whichPointInBox != -1)
            {
                // TENSOR-BASED
                //~~~~~~~~~~~~~
                facePointDerivs[pI] =
                    transformationTensorDxDb(globalIndex)
                   *volumeDerivativeCP
                    (
                        parametricCoordinates[globalIndex],
                        cpI
                    );
 
            }
        }
        dxdb[fI] =
            deltaBound.makeFaceCentresAndAreas_d
            (
                facePoints,
                facePointDerivs
            )[0];
    }
 
    return tdxdb;
}
 
 
Foam::tmp<Foam::vectorField> Foam::NURBS3DVolume::coordinates
(
    const vectorField& uVector
) const
{
    const label degreeU = basisU_.degree();
    const label degreeV = basisV_.degree();
    const label degreeW = basisW_.degree();
 
    const label nCPsU = basisU_.nCPs();
    const label nCPsV = basisV_.nCPs();
    const label nCPsW = basisW_.nCPs();
 
    const label nPoints = mapPtr_().size();
    auto tpoints = tmp<vectorField>::New(nPoints, Zero);
    auto& points = tpoints.ref();
 
    forAll(points, pI)
    {
        const label globalPI = mapPtr_()[pI];
        const scalar u = uVector[globalPI].x();
        const scalar v = uVector[globalPI].y();
        const scalar w = uVector[globalPI].z();
        for (label iCPw = 0; iCPw < nCPsW; iCPw++)
        {
            for (label iCPv = 0; iCPv < nCPsV; iCPv++)
            {
                for (label iCPu = 0; iCPu < nCPsU; iCPu++)
                {
                    points[pI] +=
                       cps_[getCPID(iCPu, iCPv, iCPw)]
                      *basisU_.basisValue(iCPu, degreeU, u)
                      *basisV_.basisValue(iCPv, degreeV, v)
                      *basisW_.basisValue(iCPw, degreeW, w);
                }
            }
        }
    }
 
    return tpoints;
}
 
 
Foam::vector Foam::NURBS3DVolume::coordinates
(
    const vector& uVector
) const
{
    const label degreeU = basisU_.degree();
    const label degreeV = basisV_.degree();
    const label degreeW = basisW_.degree();
 
    const label nCPsU = basisU_.nCPs();
    const label nCPsV = basisV_.nCPs();
    const label nCPsW = basisW_.nCPs();
 
    const scalar u = uVector.x();
    const scalar v = uVector.y();
    const scalar w = uVector.z();
 
    vector point(Zero);
    for (label iCPw = 0; iCPw < nCPsW; iCPw++)
    {
        for (label iCPv = 0; iCPv < nCPsV; iCPv++)
        {
            for (label iCPu = 0; iCPu < nCPsU; iCPu++)
            {
                point +=
                    cps_[getCPID(iCPu, iCPv, iCPw)]
                   *basisU_.basisValue(iCPu, degreeU, u)
                   *basisV_.basisValue(iCPv, degreeV, v)
                   *basisW_.basisValue(iCPw, degreeW, w);
            }
        }
    }
 
    return point;
}
 
 
Foam::tmp<Foam::vectorField> Foam::NURBS3DVolume::computeNewPoints
(
    const vectorField& controlPointsMovement
)
{
    // Get parametric coordinates and map
    const vectorField& paramCoors = getParametricCoordinates();
    const labelList& map = mapPtr_();
 
    // 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);
    const vectorField& parameterizedPoints = tparameterizedPoints();
 
    // Return field. Initialized with current mesh points
    tmp<vectorField> tnewPoints(new vectorField(mesh_.points()));
    vectorField& newPoints = tnewPoints.ref();
 
    // Update position of parameterized points
    forAll(parameterizedPoints, pI)
    {
        newPoints[map[pI]] = transformPointToCartesian(parameterizedPoints[pI]);
    }
 
    // Update coordinates in the local system based on the cartesian points
    updateLocalCoordinateSystem(newPoints);
    DebugInfo
        << "Max mesh movement equal to "
        << gMax(mag(newPoints - mesh_.points())) << endl;
 
    return tnewPoints;
}
 
 
Foam::tmp<Foam::vectorField> Foam::NURBS3DVolume::computeNewBoundaryPoints
(
    const vectorField& controlPointsMovement,
    const labelList& patchesToBeMoved
)
{
    // Get parametric coordinates
    const vectorField& paramCoors = getParametricCoordinates();
 
    // Update control points position
    cps_ += controlPointsMovement;
 
    writeCps("cpsBsplines"+mesh_.time().timeName());
    writeCpsInDict();
 
    // Return field. Initialized with current mesh points
    tmp<vectorField> tnewPoints(new vectorField(mesh_.points()));
    vectorField& newPoints = tnewPoints.ref();
 
    // Update position of parameterized boundary points
    for (const label patchI : patchesToBeMoved)
    {
        const polyPatch& patch = mesh_.boundaryMesh()[patchI];
        const labelList& meshPoints = patch.meshPoints();
 
        for (const label globalIndex : meshPoints)
        {
            const label whichPointInBox = reverseMapPtr_()[globalIndex];
            // If point resides within control points box,
            // compute new cartesian coordinates
            if (whichPointInBox != -1)
            {
                newPoints[globalIndex] =
                    transformPointToCartesian
                    (
                        coordinates
                        (
                            paramCoors[globalIndex]
                        )
                    );
            }
        }
    }
 
    // Update coordinates in the local system based on the cartesian points
    updateLocalCoordinateSystem(newPoints);
    DebugInfo
        << "Max mesh movement equal to "
        << gMax(mag(newPoints - mesh_.points())) << endl;
 
    return tnewPoints;
}
 
 
void Foam::NURBS3DVolume::setControlPoints(const vectorField& newCps)
{
    if (cps_.size() != newCps.size())
    {
        FatalErrorInFunction
           << "Attempting to replace control points with a set of "
           << "different size"
           << exit(FatalError);
    }
    cps_ = newCps;
}
 
 
void Foam::NURBS3DVolume::boundControlPointMovement
(
    vectorField& controlPointsMovement
)
{
    forAll(controlPointsMovement, cpI)
    {
        if (!activeDesignVariables_[3*cpI])
        {
            controlPointsMovement[cpI].x() = Zero;
        }
        if (!activeDesignVariables_[3*cpI + 1])
        {
            controlPointsMovement[cpI].y() = Zero;
        }
        if (!activeDesignVariables_[3*cpI + 2])
        {
            controlPointsMovement[cpI].z() = Zero;
        }
    }
}
 
 
Foam::scalar Foam::NURBS3DVolume::computeMaxBoundaryDisplacement
(
    const vectorField& controlPointsMovement,
    const labelList& patchesToBeMoved
)
{
    // Backup old cps
    vectorField oldCPs = cps_;
    // Get parametric coordinates
    const vectorField& paramCoors = getParametricCoordinates();
    // Update control points position
    cps_ += controlPointsMovement;
    // Update position of parameterized boundary points
    scalar maxDisplacement(Zero);
    for (const label patchI : patchesToBeMoved)
    {
        const polyPatch& patch = mesh_.boundaryMesh()[patchI];
        const labelList& meshPoints = patch.meshPoints();
 
        for (const label globalIndex : meshPoints)
        {
            const label whichPointInBox = reverseMapPtr_()[globalIndex];
            // If point resides within control points box,
            // compute new cartesian coordinates
            if (whichPointInBox != -1)
            {
                vector newPoint =
                    transformPointToCartesian
                    (
                       coordinates
                       (
                           paramCoors[globalIndex]
                       )
                    );
                maxDisplacement =
                    max
                    (
                        maxDisplacement,
                        mag(newPoint - mesh_.points()[globalIndex])
                    );
            }
        }
    }
    reduce(maxDisplacement, maxOp<scalar>());
    cps_ = oldCPs;
 
    return maxDisplacement;
}
 
 
Foam::tmp<Foam::vectorField> Foam::NURBS3DVolume::getPointsInBox()
{
    if (mapPtr_.empty())
    {
        findPointsInBox(localSystemCoordinates_);
    }
    tmp<vectorField> pointsInBox
    (
        new vectorField(localSystemCoordinates_, mapPtr_())
    );
 
    return pointsInBox;
}
 
 
const Foam::labelList& Foam::NURBS3DVolume::getMap()
{
    if (mapPtr_.empty())
    {
        findPointsInBox(localSystemCoordinates_);
    }
 
    return mapPtr_();
}
 
 
const Foam::labelList& Foam::NURBS3DVolume::getReverseMap()
{
    if (reverseMapPtr_.empty())
    {
        findPointsInBox(localSystemCoordinates_);
    }
 
    return reverseMapPtr_();
}
 
 
const Foam::pointVectorField& Foam::NURBS3DVolume::getParametricCoordinates()
{
    // If not computed yet, compute parametric coordinates
    if (parametricCoordinatesPtr_.empty())
    {
        // Find mesh points inside control points box
        // if they have been identified yet
        if (mapPtr_.empty())
        {
            findPointsInBox(localSystemCoordinates_);
        }
        computeParametricCoordinates(getPointsInBox()());
    }
 
    return parametricCoordinatesPtr_();
}
 
 
Foam::tmp<Foam::pointTensorField> Foam::NURBS3DVolume::getDxDb(const label cpI)
{
    // Get parametric coordinates
    const vectorField& parametricCoordinates = getParametricCoordinates();
 
    // Set return field to zero
    tmp<pointTensorField> tDxDb
    (
        new pointTensorField
        (
            IOobject
            (
                "DxDb",
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            pointMesh::New(mesh_),
            dimensionedTensor(dimless, Zero)
        )
    );
 
    pointTensorField& DxDb = tDxDb.ref();
 
    // All points outside the control box remain unmoved.
    // Loop over only points within the control box
    const labelList& map = mapPtr_();
    for (const label globalIndex : map)
    {
        DxDb[globalIndex] =
            transformationTensorDxDb(globalIndex)
           *volumeDerivativeCP
            (
                parametricCoordinates[globalIndex],
                cpI
            );
    }
 
    return tDxDb;
}
 
 
Foam::tmp<Foam::volTensorField> Foam::NURBS3DVolume::getDxCellsDb
(
    const label cpI
)
{
    // Get parametric coordinates
    const vectorField& parametricCoordinates = getParametricCoordinates();
 
    // Set return field to zero
    tmp<volTensorField> tDxDb
    (
        new volTensorField
        (
            IOobject
            (
                "DxDb",
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh_,
            dimensionedTensor(dimless, Zero)
        )
    );
 
    volTensorField& DxDb = tDxDb.ref();
    deltaBoundary deltaBound(mesh_);
    const labelListList& pointCells = mesh_.pointCells();
 
    // All points outside the control box remain unmoved.
    // Loop over only points within the control box
    const labelList& map = mapPtr_();
    for (const label globalIndex : map)
    {
        tensor pointDxDb =
            transformationTensorDxDb(globalIndex)
           *volumeDerivativeCP
            (
                parametricCoordinates[globalIndex],
                cpI
            );
        const labelList& pointCellsI = pointCells[globalIndex];
        tmp<tensorField> tC_d = deltaBound.cellCenters_d(globalIndex);
        const tensorField& C_d = tC_d();
        forAll(pointCellsI, cI)
        {
            const label cellI = pointCellsI[cI];
            DxDb[cellI] += C_d[cI] & pointDxDb;
        }
    }
 
    // Assign boundary values since the grad of this field is often needed
    forAll(mesh_.boundary(), pI)
    {
        const fvPatch& patch = mesh_.boundary()[pI];
        if (!isA<coupledFvPatch>(patch))
        {
            DxDb.boundaryFieldRef()[pI] = patchDxDbFace(pI, cpI);
        }
    }
 
    // Correct coupled boundaries
    DxDb.correctBoundaryConditions();
 
    return tDxDb;
}
 
 
Foam::label Foam::NURBS3DVolume::nUSymmetry() const
{
    label nU(basisU_.nCPs());
    if (nU % 2 == 0)
    {
        nU /=2;
    }
    else
    {
        nU = (nU - 1)/2 + 1;
    }
    return nU;
}
 
 
Foam::label Foam::NURBS3DVolume::nVSymmetry() const
{
    label nV(basisV_.nCPs());
    if (nV % 2 == 0)
    {
        nV /=2;
    }
    else
    {
        nV = (nV - 1)/2 + 1;
    }
    return nV;
}
 
 
Foam::label Foam::NURBS3DVolume::nWSymmetry() const
{
    label nW(basisW_.nCPs());
    if (nW % 2 == 0)
    {
        nW /=2;
    }
    else
    {
        nW = (nW - 1)/2 + 1;
    }
    return nW;
}
 
 
void Foam::NURBS3DVolume::writeCps(const fileName& baseName) const
{
    const label nCPsU = basisU_.nCPs();
    const label nCPsV = basisV_.nCPs();
 
    vectorField cpsInCartesian(cps_.size(), Zero);
    forAll(cpsInCartesian, cpI)
    {
        cpsInCartesian[cpI] = transformPointToCartesian(cps_[cpI]);
    }
 
    Info<< "Writing control point positions to file" << endl;
 
    if (Pstream::master())
    {
        OFstream cpsFile("optimisation"/cpsFolder_/name_ + baseName + ".csv");
        // Write header
        cpsFile
            << "\"Points : 0\", \"Points : 1\", \"Points : 2\","
            << "\"i\", \"j\", \"k\","
            << "\"active : 0\", \"active : 1\", \"active : 2\"" << endl;
 
        forAll(cpsInCartesian, cpI)
        {
            const label iCPw = cpI/label(nCPsU*nCPsV);
            const label iCPv = (cpI - iCPw*nCPsU*nCPsV)/nCPsU;
            const label iCPu = (cpI - iCPw*nCPsU*nCPsV - iCPv*nCPsU);
 
            cpsFile
                << cpsInCartesian[cpI].x() << ", "
                << cpsInCartesian[cpI].y() << ", "
                << cpsInCartesian[cpI].z() << ", "
                << iCPu << ", "
                << iCPv << ", "
                << iCPw << ", "
                << activeDesignVariables_[3*cpI] << ", "
                << activeDesignVariables_[3*cpI + 1] << ", "
                << activeDesignVariables_[3*cpI + 2] << endl;
        }
    }
}
 
 
void Foam::NURBS3DVolume::writeCpsInDict() const
{
    if (Pstream::master())
    {
        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();
}
 
 
// ************************************************************************* //