#include <iostream>

#define DEBUG 0


struct point 
{
    float x1;
    float x2;

    point() {}
    point(float x, float y) :x1(x), x2(y) {}
};

struct _face;
struct _cell;
typedef _face face;
typedef _cell cell;

// double edge?
struct _face
{
    point* p1;
    point* p2;
    struct _cell *c0, *c1;
    int face_id;
    _face() {
        p1 = NULL;
        p2 = NULL;
        c0 = NULL;
        c1 = NULL;
    }

    _face(point *pt1, point *pt2): p1(pt1), p2(pt2){
        c1 = c0 = NULL;
#if DEBUG
        printf("face: (%.0f,%.0f)->(%.0f,%.0f)\n", p1->x1, p1->x2, p2->x1, p2->x2);
#endif
    }
};

struct _cell
{
    struct _face *fw;
    struct _face *fn;
    struct _face *fe;
    struct _face *fs;
    int cell_id;
    _cell() {
        fs = fe = fn = fw = NULL;
    }
    _cell(face *f1, face *f2, face *f3, face *f4)
    : fw(f1), fn(f2), fe(f3), fs(f4) {
#if DEBUG
        printf("cell: (%.0f,%.0f)->(%.0f,%.0f)\n   ->(%.0f,%.0f)->(%.0f,%.0f)\n",
         fw->p1->x1, fw->p1->x2, fw->p2->x1, fw->p2->x2, fe->p1->x1, fe->p1->x2, fe->p2->x1, fe->p2->x2);
#endif
    }
};


#define NXMAX 7
#define NYMAX 4

struct Grid
{
    float xmin, xmax;
    float ymin, ymax;

    point points[NXMAX+1][NYMAX+1];
    face faces_row[NXMAX+1][NYMAX+1];
    face faces_col[NXMAX+1][NYMAX+1];
    cell cells[NXMAX][NYMAX];

    Grid() {
        xmin = 0.0;
        ymin = 0.0;
        xmax = 7.0;
        ymax = 4.0;
        create(1.0, 1.0);
    }

    void create(float dx, float dy);

};

void Grid::create(float dx, float dy) {
    float x = xmin; 
    for(int i=0; i < NXMAX + 1; ++i) {
        float y = ymin;
        for(int j=0; j < NYMAX + 1; ++j) {
            points[i][j] = point(x, y);
            y += dy;
            // std::cout << i << "," << j << ":" << x << "," << y << std::endl;
        }
        x += dx;
    }

    int count = 0;
    for(int i=0; i < NXMAX; ++i) {
        for(int j=0; j < NYMAX+1; ++j) {
            faces_row[i][j] = face(&points[i][j], &points[i+1][j]);
            faces_row[i][j].face_id = count++;
            if(j == 0) {
                faces_row[i][j].c0 = &cells[i][j];
            } else if (j == NYMAX) {
                faces_row[i][j].c0 = &cells[i][j-1];
            } else {
                faces_row[i][j].c0 = &cells[i][j-1];
                faces_row[i][j].c1 = &cells[i][j];
            }
        }
    }

    for(int j=0; j < NYMAX; ++j) {
        for(int i=0; i < NXMAX+1; ++i) {
            // std::cout << i << "," << j ;
            faces_col[i][j] = face(&points[i][j], &points[i][j+1]); 
            faces_row[i][j].face_id = count++;
            if(i == 0) {
                faces_col[i][j].c0 = &cells[i][j];
            } else if (i == NXMAX) {
                faces_col[i][j].c0 = &cells[i-1][j];
            } else {
                faces_col[i][j].c0 = &cells[i-1][j];
                faces_col[i][j].c1 = &cells[i][j];
            }
        }
    }

    count = 0;
    for(int i=0; i < NXMAX ; ++i) {
        for(int j=0; j < NYMAX ; ++j) {
            // std::cout << i <<"," << j <<" ";
            cells[i][j] = cell(&faces_col[i][j], &faces_row[i][j+1], &faces_col[i+1][j], &faces_row[i][j]);
            cells[i][j].cell_id = count++;
        }
    }

    printf("mesh generated: \n %d points, %d faces, %d cells\n",  
        (NXMAX+1)*(NYMAX+1), (NXMAX+1)*NYMAX+NXMAX*(NYMAX+1), NXMAX*NYMAX);
}

#include <iostream>
#include "grid.cpp"

struct face_node 
{
    face* f;
    face_node *prev, *next;
    face_node(face* fc) {
        f = fc;
        prev = NULL;
        next = NULL;
    }
};


struct face_t_iterator
{
    face_node *current;
    face_t_iterator(face_node *f) {
        current = f;
    }

    face_t_iterator(const face_t_iterator& other) {
        current = other.current;
    }

    face_t_iterator operator++(int) {
        current = current->next;
        return *this;
    }

    face_t_iterator operator--(int) {
        current = current->prev;
        return *this;
    }

    face* operator->() {
        return current->f;
    }

    bool operator==(const face_t_iterator & rhs) const
        { return current == rhs.current; }

    bool operator!=(const face_t_iterator & rhs) const
        { return !( *this == rhs ); }
};


struct face_t
{
    int _size;
    face_node *head;
    face_node *tail;

    face_t(){
        init();
    }

    face_t(face_node *f){
        init();
        push_back(f);
    }


    void init() {
        _size = 0;
        point p1(0,0);
        face f1(&p1, &p1);
        head = new face_node(&f1);
        tail = new face_node(&f1);
        head->next = tail;
        tail->prev = head;
    }

    face_t_iterator begin() const {
        return face_t_iterator(head->next);
    }

    face_t_iterator end() const {
        return face_t_iterator(tail);
    }

    void push_back(face_node* f) {
        face_node* p = end().current;
        _size++;
        f->prev = p->prev;  
        f->next = p;
        p->prev->next = f;
        p->prev = f;
    }
};

/***********************************************/

struct cell_node 
{
    cell *c;
    cell_node *prev, *next;
    cell_node(cell *cl) {
        c = cl;
        prev = NULL;
        next = NULL;
    }
};

struct cell_t_iterator
{
    cell_node *current;
    cell_t_iterator(cell_node* c) {
        current = c;
    }

    cell_t_iterator(const cell_t_iterator& other) {
        current = other.current;
    }

    cell_t_iterator operator++(int) {
        current = current->next;
        return *this;
    }

    cell_t_iterator operator--(int) {
        current = current->prev;
        return *this;
    }

    cell* operator->() {
        return current->c;
    }

    bool operator==(const cell_t_iterator & rhs) const
        { return current == rhs.current; }

    bool operator!=(const cell_t_iterator & rhs) const
        { return !( *this == rhs ); }
};

struct cell_t
{
    int _size;
    cell_node *head;
    cell_node *tail;

    cell_t(){
        init();
    }

    cell_t(cell_node *f){
        init();
        push_back(f);
    };


    void init() {
        _size = 0;
        point p1(0,0);
        face f1(&p1, &p1);
        cell c1(&f1, &f1, &f1, &f1);
        head = new cell_node(&c1);
        tail = new cell_node(&c1);
        head->next = tail;
        tail->prev = head;
    }

    cell_t_iterator begin() const {
        return cell_t_iterator(head->next);
    }

    cell_t_iterator end() const {
        return cell_t_iterator(tail);
    }

    void push_back(cell_node* f) {
        cell_node* p = end().current;
        _size++;
        f->prev = p->prev;  
        f->next = p;
        p->prev->next = f;
        p->prev = f;
    }
};

/***********************************************/
#define Nfaces (NXMAX+1)*NYMAX + NYMAX*(NXMAX+1)
#define Ncells NXMAX * NYMAX

struct domain 
{
    face_t f_thread[NYMAX+1 + NXMAX+1];
    cell_t c_thread[NYMAX];

    enum variable{
        u, v, rho, p,
        F0, F1, F2, F3, 
        G0, G1, G2, G3,
        NumberUDS
    };

    float SV_F[NumberUDS][Nfaces];
    float SV_C[NumberUDS][Ncells];

    domain(Grid* G);
    void set_domain(variable U, float value);
    void set_boundary(variable U, float value, int i);

};


domain::domain(Grid* G) {
    int i, j;
    for(i=0; i < NXMAX; ++i) {
        for(j=0; j < NYMAX+1; ++j) {
            face_node* f1 = new face_node(&G->faces_row[i][j]);
            f_thread[j].push_back(f1);
        }
    }

    for(j=0; j < NYMAX; ++j) {
        for(i=0; i < NXMAX+1; ++i) {
            face_node* f1 = new face_node(&G->faces_col[i][j]);
            f_thread[(NYMAX+1) + i].push_back(f1);
        }
    }

    for(i=0; i<NXMAX; ++i) {
        for(j=0; j < NYMAX; ++j) {
            cell_node* c1 = new cell_node(&G->cells[i][j]);
            c_thread[i].push_back(c1);
        }
    }

    // for(j=0; j < NYMAX; ++j) {
    //     for(i=0, i < NXMAX, ++i) {
    //         cell_node* c1 = new cell_node(cell[i][j]);
    //         c_thread[j].push_back(c1);
    //     }
    // }
}

void domain::set_domain(variable U, float value) {
    for(int i = 0; i < NYMAX+1; ++i) {
        cell_t t = c_thread[i];
        cell_t_iterator itr = t.begin();
        while(itr != t.end()) {
            int index = itr->cell_id;
            SV_C[U][index] = value;
            itr++;
        }
    }
}

void domain::set_boundary(variable U, float value, int i) {
    face_t t = f_thread[i];
    face_t_iterator itr = t.begin();
    if(itr->c1 != NULL) {
        printf("error: set_boundary(): invalid thread id!\n"); 
    }
    while(itr != t.end()) {
        int index = itr->face_id;
        SV_F[U][index] = value;
        itr++;
    }
}


int main()
 {

 }