  * @param x X-coordinates of interpolation points
  * @param y Y-coordinates of interpolation points
  * @param g The Grid on which to operate
- * @param globalbcz NEU for common interpolation. DIR for zeros at Box
+ * @param boundary determines what to do when a point
+ lies exactly on the boundary:  DIR generates zeroes in the interpolation matrix, NEU and other boundaries interpolate the inner side polynomial.
+ * @attention all points (x,y) must lie within or on the boundaries of g.
  *
  * @return interpolation matrix
  */
-cusp::coo_matrix<int, double, cusp::host_memory> interpolation( const thrust::host_vector<double>& x, const thrust::host_vector<double>& y, const Grid2d& g , dg::bc globalbcz = dg::NEU)
+cusp::coo_matrix<int, double, cusp::host_memory> interpolation( const thrust::host_vector<double>& x, const thrust::host_vector<double>& y, const Grid2d& g , dg::bc boundary = dg::NEU)
 {
     assert( x.size() == y.size());
     std::vector<double> gauss_nodes = g.dlt().abscissas(); 

    cusp::coo_matrix<int, double, cusp::host_memory> A( x.size(), g.size(), values.size());
    A.row_indices = row_indices; A.column_indices = column_indices; A.values = values;

    if (globalbcz == DIR_NEU ) std::cerr << "DIR_NEU NOT IMPLEMENTED "<<std::endl;
    if (globalbcz == NEU_DIR ) std::cerr << "NEU_DIR NOT IMPLEMENTED "<<std::endl;
    if (globalbcz == dg::PER ) std::cerr << "PER NOT IMPLEMENTED "<<std::endl;
    
    return A;
}

 * @param y Y-coordinates of interpolation points
 * @param z Z-coordinates of interpolation points
 * @param g The Grid on which to operate
 * @param boundary determines what to do when a point
 lies exactly on the boundary:  DIR generates zeroes in the interpolation matrix, NEU and other boundaries interpolate the inner side polynomial.
 *
 * @return interpolation matrix
 * @attention all points (x, y, z) must lie within or on the boundaries of g
 */
cusp::coo_matrix<int, double, cusp::host_memory> interpolation( const thrust::host_vector<double>& x, const thrust::host_vector<double>& y, const thrust::host_vector<double>& z, const Grid3d& g, dg::bc boundary= dg::NEU)
{
    assert( x.size() == y.size());
    assert( y.size() == z.size());
