cartesian_grid.h++

Go to the documentation of this file.

/*
 *

 Copyright (C) 2009 Joseph Pingenot

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU Affero General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.

 This program 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 Affero General Public License for more details.

 You should have received a copy of the GNU Affero General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

#ifndef CARTESIAN_GRID_HXX___
#define CARTESIAN_GRID_HXX___

#include <stdlib.h>
#include <stdio.h>
#include <glib.h>
#include <math.h>
#include <libmodelxx/shapes.h++>
#include <libmodelxx/function_specification.h++>

namespace kdotp {
  namespace libmodelxx {
    namespace grid {

#define CARTESIAN_GRID_GERROR_DOMAIN 309

      char* long_unsigned_to_string(long unsigned* val);

      /**
   \var componets: number of pieces of data per grid site
   \var T: type of data to store in each grid site.
   \var D number of dimensions for the grid.
      */
      template <typename T, unsigned int D, unsigned int Nc=1>
      class cartesian_grid {
      public:
  /**This function pointer is used to do something at each point,
   *  e.g. set a member of the struct at location.
   */
  typedef GError* (*set_points_callback)(T* location, void* data);

  /**This is a function pointer used whne using a builtin callback, in order to determine the actual value of T from a double-precision value (i.e. the value of the function at that point.
   */
  typedef GError* (*get_point_value_from_double_callback)(T* location, double* value, int dim, void* data);

  /**This is used to set a location based upon a position-dependent function.
   * x is a vector to the location in the global coordinate system
   * location is a pointer to the storage at this point
   */
  typedef GError* (*set_point_function_callback_d)(double* x, int dim, T* location, void* data);
  
  /**This is used to set a location based upon a position-dependent function.
   * X is a vector to the location in grid sites
   * location is a pointer to the storage at this point
   */
  typedef GError* (*set_point_function_callback_i)(unsigned int* X, int dim, T* location, void* data);

  /**This is used when restoring from a buffer (e.g. HDF5)*/
  cartesian_grid(const unsigned* box_size, const double* site_length, long unsigned bufflen, const T* data_buffer) {
    catesian_grid_init(box_size, site_length, NULL, data_buffer, NULL, 0, NULL);
  }
  /**This is used when restoring from an indexed grid*/
  cartesian_grid(const T* indices, long unsigned n_indices, const cartesian_grid<long unsigned, D, Nc>* indexed_grid) {
    catesian_grid_init(indexed_grid->L, indexed_grid->a, NULL, NULL, indices, n_indices, indexed_grid);
  }
  
  cartesian_grid(GList* box_size, GList* site_length, T def_value) {
    unsigned ell[D];
    double ay[D];
    unsigned int d;
    GList *b;
    GList *a_l;
    for(d=0, b=(box_size), a_l=(site_length); d<D; d++, a_l=(a_l->next), b=(b->next)) {
      if(b == NULL) {fprintf(stderr, "kdotp::libmodelxx::grid::cartesian_grid::cartesian_grid(const struct grid*): Too few box sizes given for the dimension (gave %d, dimension is =%d). Fix your grid(){} information!\n", d, D); exit(1);}
      if(a_l == NULL) {fprintf(stderr, "kdotp::libmodelxx::grid::cartesian_grid::cartesian_grid(const struct grid*): Too few site lengths given for the dimension (gave %d, dimension is =%d). Fix your grid(){} information!\n", d, D); exit(2);}
      ell[d] = *(int*)(b->data);
      ay[d] = *(double*)(a_l->data);
    }
    cartesian_grid_init(ell, ay, &def_value, NULL, NULL, 0, NULL);
  }

  cartesian_grid(const unsigned* box_size, const double* site_length, T def_value) {
    cartesian_grid_init(box_size, site_length, &def_value, NULL, NULL, 0, NULL);
  }

  /**This is the central constructor (but an init method) since constructors can't call each other*/
  void cartesian_grid_init(const unsigned* box_size, const double* site_length, T* def_value, const T* buffer, const T* indices, long unsigned n_indices, const cartesian_grid<long unsigned, D, Nc>* indexed_grid) {
    num_pts = 1;
    unsigned int d;
    for(d=0; d<D; d++) {
      L[d] = box_size[d];
      a[d] = site_length[d];
      if(a[d] < 0.0) {fprintf(stderr, "kdotp::libmodelxx::grid::cartesian_grid::cartesian_grid_init(const struct grid*): negative site length given (%f) Fix your grid(){} information!\n", a[d]); exit(3);}
      if(L[d] < 0) {fprintf(stderr, "kdotp::libmodelxx::grid::cartesian_grid::cartesian_grid_init(const struct grid*): negative box size given (%d) Fix your grid(){} information!\n", L[d]); exit(4);}
      num_pts *= L[d];
      //fprintf(stderr, "L[%d]=%d num_pts=%lu\n", d, L[d], num_pts);
    }
    if((indexed_grid != NULL) && (indexed_grid->num_pts != num_pts)) {
      fprintf(stderr, "kdotp::libmodelxx::grid::cartesian_grid::cartesian_grid_init(const struct grid*): indexed grid data size (%lu) not the same as the target's size (%lu)\n", indexed_grid->num_pts, num_pts); 
      exit(4);
    }
    data = (T*)malloc(num_pts*Nc*sizeof(T));
    /*The grid is now allocated; we now need to set up the structure.*/
    for(long unsigned i=0; i<num_pts; i++) {
      if(def_value != NULL) {
        data[i] = *def_value;
      }else if(buffer != NULL) {
        data[i] = buffer[i];
      }else{
        //if(unlikely(indexed_grid->data[i] >= n_indices)) {
        if(indexed_grid->data[i] >= n_indices) {
    fprintf(stderr, "kdotp::libmodelxx::grid::cartesian_grid::cartesian_grid_init(const struct grid*): invalid index in the index array (%lu)\n", indexed_grid->data[i]); exit(4);
        }
        data[i] = indices[indexed_grid->data[i]];
      }
    }
    /*Set up ste_x, ste_X, and the respective local symbol tables to save several operations per unit.*/
    switch(D) {
    case 1:
      ste_x.dt=::kdotp::libmodelxx::postprocessing::SCALAR;
      break;
    case 2:
      ste_x.dt=::kdotp::libmodelxx::postprocessing::VECTOR2;
      break;
    case 3:
      ste_x.dt=::kdotp::libmodelxx::postprocessing::VECTOR3;
      break;
    default:
      ste_x.dt=::kdotp::libmodelxx::postprocessing::INVALID;
    };
    ste_X.dt=ste_x.dt;
    ste_x.value = ste_X.value = NULL;
    ste_X.st=ste_x.st=::kdotp::libmodelxx::postprocessing::CONSTANT_SYMBOL;
    st_x = ::kdotp::libmodelxx::postprocessing::new_symbol_table();
    st_X = ::kdotp::libmodelxx::postprocessing::new_symbol_table();
    add_symbol_table_entry(st_x, "x", (struct ::kdotp::libmodelxx::postprocessing::symbol_table_entry*)&ste_x);
    add_symbol_table_entry(st_X, "X", (struct ::kdotp::libmodelxx::postprocessing::symbol_table_entry*)&ste_X);
  }

  /**
   *Set all components at the points that are in the shape
   *defined by the tree.
   *\param st* pointer to the shape tree
   *\param spcb callback called at each point.
   *\param privdat private data for the callback function.
   *\return GError containing the error, if any.
   */
  struct ezstruct {
    set_points_callback spcb;
    void* privdat;
  };
  static GError* set_point_function_callback_d_imp(double* x, int dim, T* location, void* data) {
    struct ezstruct* ezs = (struct ezstruct*)data;
    return ezs->spcb(location, ezs->privdat);
  }
  GError* set_points(const struct shape_tree *st, set_points_callback spcb, GHashTable* global_symbols, void* privdat=NULL) {
    struct ezstruct ezs = {spcb, privdat};
    return set_points(st, &set_point_function_callback_d_imp, NULL, global_symbols, &ezs);
  }

  /**This is implemented as a set_point_function callback
   * x is the floating-point position
   * T is the location to store the value.
   * data is a pointer to the callback data struct
   */
  struct cartesian_set_function_callback_data {
    /**The expression (function) to evaluate*/
    struct postprocessing::function_exp *exp;
    /**This gets called with the return value of the function.
     */
    get_point_value_from_double_callback cb;
    /*Holds the global symbol table*/
    GHashTable* global_symbol_table;
    GHashTable* local_symbol_table;
    struct ::kdotp::libmodelxx::postprocessing::symbol_table_entry_const *ste_x;
    /**Private data for the callback function above*/
    void* privdat;
  };
  static GError* set_point_function(double* x, int dee, T* location, void* privdat) {
    struct cartesian_set_function_callback_data* d = (struct cartesian_set_function_callback_data*)privdat;
    get_point_value_from_double_callback cb = d->cb;
    int dim = get_function_datatype(d->exp, ((struct cartesian_set_function_callback_data*)privdat)->global_symbol_table, ((struct cartesian_set_function_callback_data*)privdat)->local_symbol_table);
    if(dim < 0) return g_error_new(CARTESIAN_GRID_GERROR_DOMAIN, 200, "kdotp::libkdotp::libmodelxx::grid::cartesian_grid::set_point_function: got invalid dimension from function datatype: %d", dim);
    if(dim < 1) dim = 1;
    double* value = new double[dim];
    GError *e = NULL;
    /*Set the x value.*/
    ((struct cartesian_set_function_callback_data*)privdat)->ste_x->value = x;
    if((e=evaluate_function(value, d->exp, ((struct cartesian_set_function_callback_data*)privdat)->global_symbol_table, ((struct cartesian_set_function_callback_data*)privdat)->local_symbol_table)) != NULL) {
      GError *err = NULL;
      g_propagate_prefixed_error(&err, e, "kdotp::libkdotp::libmodelxx::grid::cartesian_grid::set_point_function: got error when evaluating function.");
      delete[] value;
      return err;
    }
    e= NULL;
    e = (*cb)(location, value, dim, d->privdat);
    if(e != NULL) {
      GError *err = NULL;
      g_propagate_prefixed_error(&err, e, "kdotp::libkdotp::libmodelxx::grid::cartesian_grid::set_point_function: got error when setting material.");
      delete[] value;
      return err;
    }
    fprintf(stderr, "x=%g->location=%p\n", x[0], (void*)location);
    delete[] value;
    return NULL;
  }
  /**
   *Set all components at the points that are in the shape
   *defined by the tree.
   *\param st* pointer to the shape tree
   *\param T* pointer to an Nc-size array of values to be set
   *\return int with the error code. Needs upgraded to GError.
   */
  GError* set_points(const struct shape_tree *st, postprocessing::function_exp *exp, GHashTable* global_symbol_table, get_point_value_from_double_callback cb, void* privdat=NULL) {
    struct cartesian_set_function_callback_data csf = {exp, cb, global_symbol_table, st_x, &ste_x, privdat};
    return set_points(st, &set_point_function, NULL, global_symbol_table, &csf);
  }

  GError* set_points(const struct shape_tree *st, set_point_function_callback_d dcb, set_point_function_callback_i icb, GHashTable* global_symbols, void* privdat = NULL) {
    /*X is the integer location*/
    unsigned int X[D];
    /*x is the floating-point location*/
    double x[D];
    unsigned int d;
    for(d=0; d<D; d++) X[d]=0;
    /*X[0] is the lowest dimension; X[D-1] is the lowest
     *Here's the idea: X can be viewed as an integer number with
     *    a digit-dependent basis.  The algorithm:
     *1) carry is set to 0 going into the loop
     *2) When the counter is updated, carry will be set to 0,
     *     unless the last "digit" overflows, in which case carry stays 1.
     *3) If carry is still 1 coming out of the X update, we're done.
     */
    GError *e = NULL;
    unsigned int carry=0;
    unsigned int c;
    while(carry == 0) {
      point_index_to_realspace(X, x);
      if(point_is_in_tree(x, st, global_symbols)) {
        //fprintf(stderr, "Y");
        /*Note that this for loop can likely be unrolled.
         *The if statement should be compiled out; it'll never
         *change at runtime.
         */
        for(c=0; c<Nc; c++) {
    if(Nc == 1) {
      if(dcb != NULL) {
        e = dcb(x, D, &data[point_index_to_storage_index(X)], privdat);
      }else{
        e = icb(X, D, &data[point_index_to_storage_index(X)], privdat);
      }
    }else{
      if(dcb != NULL) {
        e = dcb(x, D, &data[point_index_to_storage_index(X, c)], privdat);
      }else{
        e = icb(X, D, &data[point_index_to_storage_index(X, c)], privdat);
      }
    }
    if(e != NULL) {
      GError *err = NULL;
      GString *s = g_string_new("");
      GString *S = g_string_new("");
      for(unsigned int i=0; i<D; i++) {
        g_string_append_printf(s, " %g", x[i]);
        g_string_append_printf(S, " %u", X[i]);
      }
      g_propagate_prefixed_error(&err, e, "kdotp::libmodeling::grid::cartesian_grid::set_points<full>: Got error evaluating point %s (grid site %s) c=%u: ", s->str, S->str, c);
      g_string_free(s, TRUE);
      g_string_free(S, TRUE);
      return err;
    }
        }
        //}else{
        //fprintf(stderr, "N");
      }
      /*Now update our point.*/
      carry=1;
      for(d=0; ((d<D)&&(carry==1)); d++) {
        //fprintf(stderr, "%d->", X[d]);
        X[d] += carry;
        //fprintf(stderr, "%d/", X[d]);
        if(X[d] >= L[d]) {
    X[d] = 0;
    carry = 1;
        }else{
    carry = 0;
        }
      }
      //fprintf(stderr, "\n");
    }
    return 0;
  }

  /*Does something at each point.
   *\param callback. Takes unsigned int* (location in grid space) unsigned int* (location in real space), component ,dimensions, location, and private data.
   *\param privdat private data for the callback.
   *\return GError if error occurred.
   */
  GError* for_each_point(GError* (*thing_to_do)(unsigned* X, double* x, unsigned c, unsigned dim, T* location, void* privdat), void* privdat=NULL) {
    /*X is the integer location*/
    unsigned int X[D];
    /*x is the floating-point location*/
    double x[D];
    unsigned int d;
    for(d=0; d<D; d++) X[d]=0;
    /*X[0] is the lowest dimension; X[D-1] is the lowest
     *Here's the idea: X can be viewed as an integer number with
     *    a digit-dependent basis.  The algorithm:
     *1) carry is set to 0 going into the loop
     *2) When the counter is updated, carry will be set to 0,
     *     unless the last "digit" overflows, in which case carry stays 1.
     *3) If carry is still 1 coming out of the X update, we're done.
     */
    GError *e = NULL;
    unsigned int carry=0;
    unsigned int c;
    while(carry == 0) {
      point_index_to_realspace(X, x);
      //fprintf(stderr, "Y");
      /*Note that this for loop can likely be unrolled.
       *The if statement should be compiled out; it'll never
       *change at runtime.
       */
      for(c=0; c<Nc; c++) {
        if(Nc == 1) {
    e = thing_to_do(X, x, 0, D, &data[point_index_to_storage_index(X)], privdat);
        }else{
    e = thing_to_do(X, x, c, D, &data[point_index_to_storage_index(X, c)], privdat);
        }
        if(e != NULL) {
    GError *err = NULL;
    GString *s = g_string_new("");
    GString *S = g_string_new("");
    for(unsigned int i=0; i<D; i++) {
      g_string_append_printf(s, " %g", x[i]);
      g_string_append_printf(S, " %u", X[i]);
    }
    g_propagate_prefixed_error(&err, e, "kdotp::libmodeling::grid::cartesian_grid::set_points<full>: Got error evaluating point %s (grid site %s) c=%u: ", s->str, S->str, c);
    g_string_free(s, TRUE);
    g_string_free(S, TRUE);
    return err;
        }
      }
      /*Now update our point.*/
      carry=1;
      /*We should exit this loop with carry=0 unless we incremented past the end of the final dimension (i.e. we'd need to carry into the D+1 dimension*/
      for(d=0; ((d<D)&&(carry==1)); d++) {
        //fprintf(stderr, "%d->", X[d]);
        X[d] += carry;
        //fprintf(stderr, "%d/", X[d]);
        if(X[d] >= L[d]) {
    X[d] = 0;
    carry = 1;
        }else{
    carry = 0;
        }
      }
      //fprintf(stderr, "\n");
    }
    return NULL;
  }
      
  /**
   *Set all components at the points that are in the shape
   *defined by the tree.
   *\param st* pointer to the shape tree
   *\param T* pointer to an Nc-size array of values to be set
   *\return int with the error code. Needs upgraded to GError.
   */
  int set_points(const struct shape_tree *st, T* val, GHashTable* global_symbols) {
    /*X is the integer location*/
    unsigned int X[D];
    /*x is the floating-point location*/
    double x[D];
    unsigned int d;
    for(d=0; d<D; d++) X[d]=0;
    /*X[0] is the lowest dimension; X[D-1] is the lowest
     *Here's the idea: X can be viewed as an integer number with
     *    a digit-dependent basis.  The algorithm:
     *1) carry is set to 0 going into the loop
     *2) When the counter is updated, carry will be set to 0,
     *     unless the last "digit" overflows, in which case carry stays 1.
     *3) If carry is still 1 coming out of the X update, we're done.
     */
    unsigned int carry=0;
    unsigned int c;
    while(carry == 0) {
      point_index_to_realspace(X, x);
      if(point_is_in_tree(x, st, global_symbols)) {
        //fprintf(stderr, "Y");
        /*Note that this for loop can likely be unrolled.
         *The if statement should be compiled out; it'll never
         *change at runtime.
         */
        for(c=0; c<Nc; c++) {
    if(Nc == 1) {
      data[point_index_to_storage_index(X)] = val[c];
    }else{
      data[point_index_to_storage_index(X, c)] = val[c];
    }
        }
        //}else{
        //fprintf(stderr, "N");
      }
      /*Now update our point.*/
      carry=1;
      for(d=0; ((d<D)&&(carry==1)); d++) {
        //fprintf(stderr, "%d->", X[d]);
        X[d] += carry;
        //fprintf(stderr, "%d/", X[d]);
        if(X[d] >= L[d]) {
    X[d] = 0;
    carry = 1;
        }else{
    carry = 0;
        }
      }
      //fprintf(stderr, "\n");
    }
    return 0;
  }

  /**
   *Covert a point within the grid from index to realspace
   *\param X integer array of size D, holding the grid indices
   *\param x double array of size D, to hold the position
   */
  void point_index_to_realspace(const unsigned int* X, double* x) const {
    unsigned int d;
    for(d=0; d<D; d++) {
      /*Formulated to put 0 in the center of the grid*/
      x[d] = (double(X[d]) - double(L[d])/2.0)*a[d];
    }
  }
  /**
   *Covert a point (index-based) to a location within storage
   *\param X integer array of size D, holding the grid indices.
   *\param c component number (default is 0; saves a save/load
   *if there's just one!)
   *\return long unsigned int, index into the storage.
   **/
  inline long unsigned int point_index_to_storage_index(const unsigned int* X, unsigned int c=0) const {
    unsigned int d;
    /*We initialize with c, the component, already done.*/
    unsigned long int block_offset = Nc;
    unsigned long int index = c;
    for(d=0; d<D; d++) {
      index += X[d]*block_offset;
      block_offset *= L[d];
    }
    return index;
  }

  bool point_is_in_tree(const double* x, const struct shape_tree *st, GHashTable* global_syms) const {
    bool l, r;
    switch(st->op) {
    case OR:
      return point_is_in_tree(x, st->l, global_syms) || point_is_in_tree(x, st->r, global_syms);
    case XOR:
      l = point_is_in_tree(x, st->l, global_syms);
      r = point_is_in_tree(x, st->r, global_syms);
      return (l || r) && (!(l && r));
    case AND:
      return point_is_in_tree(x, st->l, global_syms) && point_is_in_tree(x, st->r, global_syms);
    case NOT:
      return !point_is_in_tree(x, st->l, global_syms);
    case LEAF:
      break;
    default:
      fprintf(stderr, "kdotp::libmodelxx::grid::cartesian_grid::point_is_in_tree: unknown shape tree operation %d. Assuming point NOT in tree.\n", st->op);
      return false;
    }
    /*This code only executes if we are in a leaf*/
    struct shape* s = (struct shape*)st;
    switch(s->type) {
    case GLOBAL:
      /*A "global" shape is everywhere.*/
      return true;
    case BLOCK:
      return point_is_in_shape(x, (struct shape_block*)(s), global_syms);
    case ABSBLOCK:
      return point_is_in_shape(x, (struct shape_block*)(s), global_syms);
    default:
      fprintf(stderr, "kdotp::libmodelxx::grid::cartesian_grid::point_is_in_tree: unknown shape type %d. Assuming point NOT in tree.\n", s->type);
      return false;
    }
    /*If we've hit this point, we've gone too far.*/
    return false;
  }
  /**
   *Determines if a point is within the given shape.
   *\note TODO_FOR_SHAPE
   *\note Overload point_is_in_shape to implement a new shape.*/
  bool point_is_in_shape(const double* x, const struct shape_block* b, GHashTable* global_syms) const {
      /*Check if the point lies within the box (dist to center
       * less than half edge length, for each edge.*/
      //fprintf(stderr, "piis: x[%d]=%g center[%d]=%g L[%d]=%g\n", d, x[d], d, (b->center)[d], d, (b->L)[d]);
      /*fprintf(stderr, "point_is_in_shape: %g: center=%g L=%g x-center=%g L/2=%g: ", x[d], b->center[d], b->L[d], fabs(x[d] - b->center[d]), fabs(b->L[d]/2.0));*/
      double center[D];
      double L[D];
      if(get_function_datatype(b->center, global_syms, NULL) != D) {
        return false;
      }
      if(get_function_datatype(b->L, global_syms, NULL) != D) {
        return false;
      }
      GError *e = evaluate_function(center, b->center, global_syms, NULL);
      if(e != NULL) {
        fprintf(stderr, "ERROR evaluating center function for block shape: %s", e->message);
        g_error_free(e);
        return false;
      }
      e = evaluate_function(L, b->L, global_syms, NULL);
      if(e != NULL) {
        fprintf(stderr, "ERROR evaluating size function for block shape: %s", e->message);
        g_error_free(e);
        return false;
      }
    for(unsigned int d = 0; d<D; d++) {
      switch (b->type) {
      case BLOCK:
        if(!(fabs(x[d] - (center)[d]) <= fabs((L)[d]/2.0))) {
    return false;
        }
        break;
      case ABSBLOCK:
        if(center[d] > L[d]) {
    if(!((x[d] >= L[d]) && (x[d] <= center[d]))) return false;
        }else{
    if(!((x[d] <= L[d]) && (x[d] >= center[d]))) return false;
        }
        break;
      default:
        fprintf(stderr, "ERROR in point_is_in_shape (block argument): block isn't a block or an absblock?! (is %d).\n", b->type);
        return false;
      }
    }
    /*fprintf(stderr, "TRUE\n");*/
    return true;
  }

  /**helper function for create_indexed_unique_component_grid*/
  long unsigned get_index_for_component_and_store_if_new(GList** comps, bool (*is_same)(T* a, T* b), T* comp) {
    long unsigned index=0;
    for(GList* l=*comps; l != NULL; l=l->next, index++) {
      if(is_same(comp, (T*)l->data)) {
        return index;
      }
    }
    *comps = g_list_append(*comps, comp);
    return index;
  }
  /**Create a new grid which contains unique component indices.  The initial use case for this is creating an integer-based material map.  We then have an array of materials and a grid of material assignments.
   *\param indices pointer to a place to store an array of type T (the final array; the index stored in the grid is the array index to get the component
   *\param n_indices number of T's stored in the indices array
   *\param indexed_grid pointer to a new, identically-spec'ed Cartesian grid, but with indices into the component array instead of component.
   *\param is_same callback to compare to items (return true if same, false otherwise)
   *\returns pointer to a GError containing the error, or NULL if no error occurred.
   */
  GError* create_indexed_unique_component_grid(T** indices, long unsigned *n_indices, cartesian_grid<long unsigned, D, Nc> **indexed_grid, bool (*is_same)(T* a, T* b)) {
    *indexed_grid = new cartesian_grid<long unsigned, D, Nc>(L, a, 0);
    /*The easy way to do this is not to care about our location (it's really irrelevant; we only care about whether or not this component is unique or not.*/
    /*List of the unique components.*/
    GList *components = NULL;
    for(long unsigned i=0; i<num_pts; i++) {
      ((*indexed_grid)->data)[i] = get_index_for_component_and_store_if_new(&components, is_same, &(data[i]));
    }
    *n_indices = g_list_length(components);
    *indices = (T*)malloc(sizeof(T)*(*n_indices));
    if(*indices == NULL) {
      return g_error_new(CARTESIAN_GRID_GERROR_DOMAIN, 200, "kdotp::libkdotp::libmodelxx::grid::cartesian_grid::create_indexed_unique_component_grid: failed to allocate storage for the index array");
    }
    GList *l;
    long unsigned i=0;
    for(l=components, i=0; l!=NULL; l=l->next, i++) {
      (*indices)[i] = *(T*)(l->data);
    }
    g_list_free(components);
    return NULL;
  }


  /*Callback to implement the dump.
   */
  struct dump_foreach_callback_privdat {
    GOutputStream *gos;
    char* (*t_to_string)(T*);
    GString *s;
  };
  static GError* dump_foreach_callback(unsigned* X, double* x, unsigned c, unsigned dim, T* location, void* privdat) {
    struct dump_foreach_callback_privdat *cbp = (struct dump_foreach_callback_privdat*)privdat;
    g_string_printf(cbp->s, "%u\t%g\t", X[0], x[0]);
    for(unsigned d=1; d<dim; d++) {
      g_string_append_printf(cbp->s, "%u\t%g\t", X[d], x[d]);
    }
    char* data = cbp->t_to_string(location);
    g_string_append_printf(cbp->s, "%u\t%s\n", c, data);
    free(data);
    GError *e=NULL;
    gsize bytes_written;
    //fprintf(stderr, "str=%s\n", cbp->s->str);
    if(!g_output_stream_write_all(cbp->gos, cbp->s->str, cbp->s->len, &bytes_written, NULL, &e)) {
      GError *err=NULL;
      g_propagate_prefixed_error(&err, e, "kdotp::libmodelxx::grid::cartesian_grid.h++::dump: error writing line (%s). ", cbp->s->str);
      return err;
    }
    return NULL;
  }
  /**Dumps the grid to a file (gfos) in ASCII.
   *\param gos GOutputStream to write to.
   *\param t_to_string callback to convert T to a string.
   *\return the GError of the error that occurred, or NULL if no error.
   */
  GError* dump(GOutputStream *gos, char* (*t_to_string)(T*)) {
    struct dump_foreach_callback_privdat dcp = {gos, t_to_string, NULL};
    dcp.s = g_string_new("");
    /*Write the header.*/
    GError *e = NULL;
    gsize bytes_written;
    GString *s = g_string_new("#Grid dimensions:\n");
    if(!g_output_stream_write_all(gos, s->str, s->len, &bytes_written, NULL, &e)) {
      GError *err=NULL;
      g_propagate_prefixed_error(&err, e, "kdotp::libmodelxx::grid::cartesian_grid.h++::dump: error writing header.\n");
      return err;
    }
    g_string_printf(s, "%u", D);
    for(unsigned d=0; d<D; d++) {
      g_string_append_printf(s, "\t%u\t%g", L[d], a[d]);
    }
    g_string_append_printf(s, "\t%u\n", Nc);
    if(!g_output_stream_write_all(gos, s->str, s->len, &bytes_written, NULL, &e)) {
      GError *err=NULL;
      g_propagate_prefixed_error(&err, e, "kdotp::libmodelxx::grid::cartesian_grid.h++::dump: error writing header.\n");
      return err;
    }
    s = g_string_assign(s, "#");
    for(unsigned d=0; d<D; d++) {
      g_string_append_printf(s, "X%u\tx%u\t", d, d);
    }
    s = g_string_append(s, "\n");
    if(!g_output_stream_write_all(gos, s->str, s->len, &bytes_written, NULL, &e)) {
      GError *err=NULL;
      g_propagate_prefixed_error(&err, e, "kdotp::libmodelxx::grid::cartesian_grid.h++::dump: error writing header.\n");
      return err;
    }
    e = for_each_point(dump_foreach_callback, (void*)(&dcp));
    g_string_free(s, TRUE);
    if(e != NULL) {
      GError *err = NULL;
      g_propagate_prefixed_error(&err, e, "kdotp::libmodelxx::grid::cartesian_grid.h++::dump: error writing grid data.\n");
      return err;
    }
    return NULL;
  }

  /**Like dump, but dumps an *indexed* version.*/
  GError* dump_indexed(GOutputStream *gos, bool (*is_same)(T* a, T* b), char* (*t_to_string)(T*)) {
    T* indices;
    long unsigned n_indices;
    cartesian_grid<long unsigned, D, Nc>* indexed_grid;
    GError *e=NULL;
    e = create_indexed_unique_component_grid(&indices, &n_indices, &indexed_grid, is_same);
    if(e != NULL) {
      GError *err = NULL;
      g_propagate_prefixed_error(&err, e, "kdotp::libmodelxx::grid::cartesian_grid.h++::dump_indexed: error getting indexed grid.\n");
      return err;
    }
    GString *s = g_string_new("#material assignments:\n");
    g_string_append_printf(s, "%lu\n", n_indices);
    gsize bytes_written;
    if(!g_output_stream_write_all(gos, s->str, s->len, &bytes_written, NULL, &e)) {
      GError *err=NULL;
      g_propagate_prefixed_error(&err, e, "kdotp::libmodelxx::grid::cartesian_grid.h++::dump_indexed: error writing header.\n");
      free(indices);
      delete indexed_grid;
      return err;
    }
    for(long unsigned i=0; i<n_indices; i++) {
      char* rstr;
      g_string_printf(s, "%s\n", rstr=t_to_string(&(indices[i])));
      if(!g_output_stream_write_all(gos, s->str, s->len, &bytes_written, NULL, &e)) {
        GError *err=NULL;
        g_propagate_prefixed_error(&err, e, "kdotp::libmodelxx::grid::cartesian_grid.h++::dump_indexed: error writing material %lu (%s).\n", i, rstr);
        free(indices);
        free(rstr);
        delete indexed_grid;
        return err;
      }
      //fprintf(stderr, "index %lu: %s\n", i, rstr);
      free(rstr);
    }
    e=indexed_grid->dump(gos, long_unsigned_to_string);
    if(e != NULL) {
      GError *err=NULL;
      g_propagate_prefixed_error(&err, e, "kdotp::libmodelxx::grid::cartesian_grid.h++::dump_indexed: error dumping index grid.\n");
      free(indices);
      delete indexed_grid;
      return err;
    }
    free(indices);
    delete indexed_grid;
    return NULL;
  }

  /**Sets a point in the grid, by grid site
   *\note this is pretty much just a wrapper around using point_index_to_storage_index.
   */
  inline void set_point(unsigned* X, T* val, unsigned c=0) {
    data[point_index_to_storage_index(X, c)] = *val;
  }
  inline void set_point(unsigned* X, T& val, unsigned c=0) {
    data[point_index_to_storage_index(X, c)] = val;
  }
  inline void set_point_(unsigned* X, T& val, unsigned c=0) {
    data[point_index_to_storage_index(X, c)] = val;
  }

  /**Sets a point in the grid, by grid site
   *\note this is pretty much just a wrapper around using point_index_to_storage_index.
   *\note that we use store_pointer to get out a *pointer* to the datum, *not* a copy!
   *\note THIS DOES NOT MAKE A COPY!
   */
  inline void get_point(unsigned* X, T** store_pointer, unsigned c=0) {
    *store_pointer = &data[point_index_to_storage_index(X, c)];
  }
  /**same as get_point, but returns the pointer instead
   *\note THIS DOES NOT MAKE A COPY!
   */
  inline T* get_point_r(unsigned* X, unsigned c=0) {
    return &data[point_index_to_storage_index(X, c)];
  }

  /*\var Length (in grid sites) of each edge of the compuational box*/
  unsigned int L[D];
  /*\var distance between grid points*/
  double a[D];
  /*Data storage*/
  long unsigned int num_pts;
  T* data;
  struct ::kdotp::libmodelxx::postprocessing::symbol_table_entry_const ste_x;
  GHashTable* st_x;
  struct ::kdotp::libmodelxx::postprocessing::symbol_table_entry_const ste_X;
  GHashTable* st_X;
      private:
  /*Make these unsafe defaults not exist*/
  cartesian_grid() {}
  cartesian_grid(const cartesian_grid& g) {}
  cartesian_grid operator=(const cartesian_grid& g) {return *this;}
      };
      

    }
  }
}
#endif