ThreeDGrid.cpp

Go to the documentation of this file.

/*
 *  Copyright (c) 2012, Chris Rogers
 *  All rights reserved.
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *  1. Redistributions of source code must retain the above copyright notice,
 *     this list of conditions and the following disclaimer.
 *  2. Redistributions in binary form must reproduce the above copyright notice,
 *     this list of conditions and the following disclaimer in the documentation
 *     and/or other materials provided with the distribution.
 *  3. Neither the name of STFC nor the names of its contributors may be used to
 *     endorse or promote products derived from this software without specific
 *     prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 */
 
#include "Fields/Interpolation/ThreeDGrid.h"
#include "Utilities/LogicalError.h"
 
namespace interpolation {
 
    ThreeDGrid::ThreeDGrid()
        : x_m(2, 0),
          y_m(2, 0),
          z_m(2, 0),
          xSize_m(0),
          ySize_m(0),
          zSize_m(0),
          maps_m(0),
          constantSpacing_m(false) {
        setConstantSpacing();
        x_m[1] = y_m[1] = z_m[1] = 1.;
    }
 
    ThreeDGrid::ThreeDGrid(
            int xSize, const double* x, int ySize, const double* y, int zSize, const double* z)
        : x_m(x, x + xSize),
          y_m(y, y + ySize),
          z_m(z, z + zSize),
          xSize_m(xSize),
          ySize_m(ySize),
          zSize_m(zSize),
          maps_m(),
          constantSpacing_m(false) {
        setConstantSpacing();
        if (xSize_m < 2 || ySize_m < 2 || zSize_m < 2)
            throw(LogicalError(
                    "ThreeDGrid::ThreeDGrid(...)", "3D Grid must be at least 2x2x2 grid"));
    }
 
    ThreeDGrid::ThreeDGrid(std::vector<double> x, std::vector<double> y, std::vector<double> z)
        : x_m(x),
          y_m(y),
          z_m(z),
          xSize_m(x.size()),
          ySize_m(y.size()),
          zSize_m(z.size()),
          maps_m(),
          constantSpacing_m(false) {
        setConstantSpacing();
        if (xSize_m < 2 || ySize_m < 2 || zSize_m < 2)
            throw(LogicalError(
                    "ThreeDGrid::ThreeDGrid(...)", "3D Grid must be at least 2x2x2 grid"));
    }
 
    ThreeDGrid::ThreeDGrid(
            double dX, double dY, double dZ, double minX, double minY, double minZ,
            int numberOfXCoords, int numberOfYCoords, int numberOfZCoords)
        : x_m(numberOfXCoords),
          y_m(numberOfYCoords),
          z_m(numberOfZCoords),
          xSize_m(numberOfXCoords),
          ySize_m(numberOfYCoords),
          zSize_m(numberOfZCoords),
          maps_m(),
          constantSpacing_m(true) {
        for (int i = 0; i < numberOfXCoords; i++)
            x_m[i] = minX + i * dX;
        for (int j = 0; j < numberOfYCoords; j++)
            y_m[j] = minY + j * dY;
        for (int k = 0; k < numberOfZCoords; k++)
            z_m[k] = minZ + k * dZ;
 
        setConstantSpacing();
    }
 
    ThreeDGrid::~ThreeDGrid() {}
 
    // state starts at 1,1,1
    Mesh::Iterator ThreeDGrid::begin() const {
        return Mesh::Iterator(std::vector<int>(3, 1), this);
    }
 
    Mesh::Iterator ThreeDGrid::end() const {
        std::vector<int> end(3, 1);
        end[0] = xSize_m + 1;
        return Mesh::Iterator(end, this);
    }
 
    void ThreeDGrid::getPosition(const Mesh::Iterator& it, double* position) const {
        position[0] = x(it.state_m[0]);
        position[1] = y(it.state_m[1]);
        position[2] = z(it.state_m[2]);
    }
 
    Mesh* ThreeDGrid::dual() const {
        std::vector<double> new_x(x_m.size() - 1);
        std::vector<double> new_y(y_m.size() - 1);
        std::vector<double> new_z(z_m.size() - 1);
        for (size_t i = 0; i < x_m.size() - 1; ++i) {
            new_x[i] = (x_m[i] + x_m[i + 1]) / 2.;
        }
        for (size_t i = 0; i < y_m.size() - 1; ++i) {
            new_y[i] = (y_m[i] + y_m[i + 1]) / 2.;
        }
        for (size_t i = 0; i < z_m.size() - 1; ++i) {
            new_z[i] = (z_m[i] + z_m[i + 1]) / 2.;
        }
        return new ThreeDGrid(new_x, new_y, new_z);
    }
 
    Mesh::Iterator& ThreeDGrid::addEquals(Mesh::Iterator& lhs, int difference) const {
        if (difference < 0) return subEquals(lhs, -difference);
 
        lhs.state_m[0] += difference / (ySize_m * zSize_m);
        difference = difference % (ySize_m * zSize_m);
        lhs.state_m[1] += difference / (zSize_m);
        lhs.state_m[2] += difference % (zSize_m);
 
        if (lhs.state_m[1] > ySize_m) {
            lhs.state_m[0]++;
            lhs.state_m[1] -= ySize_m;
        }
        if (lhs.state_m[2] > zSize_m) {
            lhs.state_m[1]++;
            lhs.state_m[2] -= zSize_m;
        }
 
        return lhs;
    }
 
    Mesh::Iterator& ThreeDGrid::subEquals(Mesh::Iterator& lhs, int difference) const {
        if (difference < 0) return addEquals(lhs, -difference);
 
        lhs.state_m[0] -= difference / (ySize_m * zSize_m);
        difference = difference % (ySize_m * zSize_m);
        lhs.state_m[1] -= difference / (zSize_m);
        lhs.state_m[2] -= difference % (zSize_m);
 
        while (lhs.state_m[2] < 1) {
            lhs.state_m[1]--;
            lhs.state_m[2] += zSize_m;
        }
        while (lhs.state_m[1] < 1) {
            lhs.state_m[0]--;
            lhs.state_m[1] += ySize_m;
        }
 
        return lhs;
    }
 
    void ThreeDGrid::vectorLowerBound(std::vector<double> vec, double x, int& index) {
        if (x < vec[0]) {
            index = -1;
            return;
        }
        if (x >= vec.back()) {
            index = vec.size() - 1;
            return;
        }
        size_t xLower = 0;
        size_t xUpper = vec.size() - 1;
        while (xUpper - xLower > 1) {
            index = (xUpper + xLower) / 2;
            if (x >= vec[index]) {
                xLower = index;
            } else {
                xUpper = index;
            }
        }
        index = xLower;
    }
 
    Mesh::Iterator& ThreeDGrid::addEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const {
        return addEquals(lhs, rhs.toInteger());
    }
 
    Mesh::Iterator& ThreeDGrid::subEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const {
        return subEquals(lhs, rhs.toInteger());
    }
 
    Mesh::Iterator& ThreeDGrid::addOne(Mesh::Iterator& lhs) const {
        if (lhs.state_m[1] == ySize_m && lhs.state_m[2] == zSize_m) {
            ++lhs.state_m[0];
            lhs.state_m[1] = lhs.state_m[2] = 1;
        } else if (lhs.state_m[2] == zSize_m) {
            ++lhs.state_m[1];
            lhs.state_m[2] = 1;
        } else {
            ++lhs.state_m[2];
        }
        return lhs;
    }
 
    Mesh::Iterator& ThreeDGrid::subOne(Mesh::Iterator& lhs) const {
        if (lhs.state_m[1] == 1 && lhs.state_m[2] == 1) {
            --lhs.state_m[0];
            lhs.state_m[1] = ySize_m;
            lhs.state_m[2] = zSize_m;
        } else if (lhs.state_m[2] == 1) {
            --lhs.state_m[1];
            lhs.state_m[2] = zSize_m;
        } else {
            --lhs.state_m[2];
        }
        return lhs;
    }
 
    // Check relative position
    bool ThreeDGrid::isGreater(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const {
        if (lhs.state_m[0] > rhs.state_m[0])
            return true;
        else if (lhs.state_m[0] == rhs.state_m[0] && lhs.state_m[1] > rhs.state_m[1])
            return true;
        else if (
                lhs.state_m[0] == rhs.state_m[0] && lhs.state_m[1] == rhs.state_m[1]
                && lhs.state_m[2] > rhs.state_m[2])
            return true;
        return false;
    }
 
    // remove *map if it exists; delete this if there are no more VectorMaps
    void ThreeDGrid::remove(VectorMap* map) {
        std::vector<VectorMap*>::iterator it = std::find(maps_m.begin(), maps_m.end(), map);
        if (it < maps_m.end()) {
            maps_m.erase(it);
        }
        if (maps_m.begin() == maps_m.end()) {
            delete this;
        }
    }
 
    // add *map if it has not already been added
    void ThreeDGrid::add(VectorMap* map) {
        std::vector<VectorMap*>::iterator it = std::find(maps_m.begin(), maps_m.end(), map);
        if (it == maps_m.end()) {
            maps_m.push_back(map);
        }
    }
 
    void ThreeDGrid::setConstantSpacing() {
        constantSpacing_m = true;
        for (unsigned int i = 0; i < x_m.size() - 1; i++)
            if (std::abs(1 - (x_m[i + 1] - x_m[i]) / (x_m[1] - x_m[0])) > 1e-9) {
                constantSpacing_m = false;
                return;
            }
        for (unsigned int i = 0; i < y_m.size() - 1; i++)
            if (std::abs(1 - (y_m[i + 1] - y_m[i]) / (y_m[1] - y_m[0])) > 1e-9) {
                constantSpacing_m = false;
                return;
            }
        for (unsigned int i = 0; i < z_m.size() - 1; i++)
            if (std::abs(1 - (z_m[i + 1] - z_m[i]) / (z_m[1] - z_m[0])) > 1e-9) {
                constantSpacing_m = false;
                return;
            }
    }
 
    Mesh::Iterator ThreeDGrid::getNearest(const double* position) const {
        std::vector<int> index(3);
        lowerBound(position[0], index[0], position[1], index[1], position[2], index[2]);
        if (index[0] < xSize_m - 1 && index[0] >= 0)
            index[0] +=
                    (2 * (position[0] - x_m[index[0]]) > x_m[index[0] + 1] - x_m[index[0]] ? 2 : 1);
        else
            index[0]++;
        if (index[1] < ySize_m - 1 && index[1] >= 0)
            index[1] +=
                    (2 * (position[1] - y_m[index[1]]) > y_m[index[1] + 1] - y_m[index[1]] ? 2 : 1);
        else
            index[1]++;
        if (index[2] < zSize_m - 1 && index[2] >= 0)
            index[2] +=
                    (2 * (position[2] - z_m[index[2]]) > z_m[index[2] + 1] - z_m[index[2]] ? 2 : 1);
        else
            index[2]++;
        if (index[0] < 1) index[0] = 1;
        if (index[1] < 1) index[1] = 1;
        if (index[2] < 1) index[2] = 1;
        if (index[0] > xSize_m) index[0] = xSize_m;
        if (index[1] > ySize_m) index[1] = ySize_m;
        if (index[2] > zSize_m) index[2] = zSize_m;
        return Mesh::Iterator(index, this);
    }
}  // namespace interpolation