Line data    Source code

       1             : ///////////////////////////////////////////////////////////////////////////////
       2             : // Copyright 2019 Lakshman Anumolu, Raunak Bardia.
       3             : //
       4             : // Redistribution and use in source and binary forms, with or without
       5             : // modification, are permitted provided that the following conditions are
       6             : // met:
       7             : //
       8             : // 1. Redistributions of source code must retain the above copyright notice,
       9             : // this list of conditions and the following disclaimer.
      10             : //
      11             : // 2. Redistributions in binary form must reproduce the above copyright notice,
      12             : // this list of conditions and the following disclaimer in the documentation
      13             : // and/or other materials provided with the distribution.
      14             : //
      15             : // 3. Neither the name of the copyright holder nor the names of its contributors
      16             : // may be used to endorse or promote products derived from this software without
      17             : // specific prior written permission.
      18             : //
      19             : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
      20             : // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
      21             : // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
      22             : // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
      23             : // HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
      24             : // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
      25             : // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
      26             : // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
      27             : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
      28             : // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
      29             : // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
      30             : ///////////////////////////////////////////////////////////////////////////////
      31             : 
      32             : #include "gals/cpu/interpolation/linear.h"
      33             : 
      34             : #include <math.h>
      35             : 
      36             : template <typename T>
      37         115 : T GALS::INTERPOLATION::Linear<T, GALS::CPU::Grid<T, 1>>::linearInterpolation(
      38             :     const GALS::CPU::Grid<typename GALS::CPU::Grid<T, 1>::value_type, GALS::CPU::Grid<T, 1>::dim> &grid,
      39             :     const typename GALS::CPU::Grid<T, 1>::position_type &x_interp,
      40             :     const GALS::CPU::Array<GALS::CPU::Grid<T, 1>, T> &alpha)
      41             : {
      42             :   typedef GALS::CPU::Grid<T, 1> T_GRID;
      43             : 
      44             :   const int dim = T_GRID::dim;
      45             :   const auto &axis_vectors = GALS::CPU::Grid<typename T_GRID::value_type, T_GRID::dim>::axis_vectors;
      46         230 :   const GALS::CPU::Vec3<int> base_node_id = grid.baseNodeId(x_interp);
      47         230 :   const GALS::CPU::Vec3<int> base_node_id_p1 = GALS::CPU::Vec3<int>(base_node_id[0] + axis_vectors(0, 0), 0, 0);
      48             : 
      49         230 :   const typename T_GRID::position_type x_base = grid(base_node_id);
      50         115 :   const auto &one_over_dx = grid.oneOverDX();
      51             : 
      52         115 :   T eta = (x_interp[0] - x_base[0]) * one_over_dx[0];
      53             : 
      54         115 :   T alpha_interpolated = (1. - eta) * alpha(base_node_id) + eta * alpha(base_node_id_p1);
      55             : 
      56             :   // TODO (lakshman): DELETE
      57             :   // std::cout << "base_node_id = " << base_node_id << "; eta = " << eta
      58             :   //<< "; alpha_interpolation = " << alpha_interpolated << std::endl;
      59             :   // T xo = 0., ro = 0.5;
      60             :   // T alpha_exact = (x_interp[0] - xo) * (x_interp[0] - xo) - ro * ro;
      61             :   // std::cout << "error = " << fabs(alpha_interpolated - alpha_exact) << std::endl;
      62             : 
      63         115 :   return alpha_interpolated;
      64             : }
      65             : 
      66             : template <typename T>
      67           5 : void GALS::INTERPOLATION::Linear<T, GALS::CPU::Grid<T, 1>>::compute(
      68             :     const GALS::CPU::Array<GALS::CPU::Grid<T, 1>, typename GALS::CPU::Grid<T, 1>::position_type> &x_interp,
      69             :     const GALS::CPU::Array<GALS::CPU::Grid<T, 1>, T> &alpha,
      70             :     GALS::CPU::Array<GALS::CPU::Grid<T, 1>, T> &alpha_interpolated)
      71             : {
      72             :   typedef GALS::CPU::Grid<T, 1> T_GRID;
      73             : 
      74          10 :   const GALS::CPU::Vec3<int> num_cells = alpha.numCells();
      75           5 :   const T_GRID &grid = alpha.grid();
      76          10 :   const GALS::CPU::Vec3<typename T_GRID::value_type> dx = grid.dX();
      77             :   const auto &axis_vectors = GALS::CPU::Grid<typename T_GRID::value_type, T_GRID::dim>::axis_vectors;
      78             : 
      79             :   // for (int i = 0; i < num_cells[0]; ++i)
      80             :   // for (int j = 0; j < num_cells[1]; ++j)
      81             :   // for (int k = 0; k < num_cells[2]; ++k)
      82             :   // alpha_interpolated(i, j, k) = linearInterpolation(grid, x_interp(i, j, k), alpha);
      83         235 :   for (int i = 0; i < x_interp.grid().numCells()[0]; ++i)
      84         345 :     for (int j = 0; j < x_interp.grid().numCells()[1]; ++j)
      85         345 :       for (int k = 0; k < x_interp.grid().numCells()[2]; ++k)
      86         115 :         alpha_interpolated(i, j, k) = linearInterpolation(grid, x_interp(i, j, k), alpha);
      87           5 : }
      88             : 
      89             : template <typename T>
      90           5 : void GALS::INTERPOLATION::Linear<T, GALS::CPU::Grid<T, 1>>::compute(
      91             :     const GALS::CPU::Array<GALS::CPU::Grid<T, 1>, typename GALS::CPU::Grid<T, 1>::position_type> &x_interp,
      92             :     GALS::CPU::Levelset<GALS::CPU::Grid<T, 1>, T> &levelset)
      93             : {
      94             :   GALS_FUNCTION_NOT_IMPLEMENTED("linear-1d.cc: compute(...) with levelset.")
      95           5 : }
      96             : 
      97           3 : template class GALS::INTERPOLATION::Linear<double, GALS::CPU::Grid<double, 1>>;