#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "Memory.h"
#include "FileIO.h"
#include "Mesh.h"
#include "Progress.h"
#include "Solver.h"
#include "Triangle.h"
#include "Types.h"
#include "Implicit.h"

uint index2element(uint i,uint j,uint k){
    if(j >2 || j<0) printf("index2element-2nd is bad: %d is not in [0,2]\n",j);
    if(k >3 || k<0) printf("index2element-3rd is bad: %d is not in [0,3]\n",k);
    return 12*i+4*j+k;
}

void element2index(uint p, uint* i,uint* j,uint* k){
  *i= p/12;
  *j=(p-12* *i)/4;
  *k= p-12* *i-4* *j;
}

int Solver_Implicit_Init(tSolver_State *s){
tSolver_Matrix_Element *me, **row;
uint i,j,k,p,q;


i=2; j=2; k=3;
p = index2element(i,j,k);
printf("Test %d %d %d %d\n",i,j,k,p);
i=j=k=0;
p-=10;
element2index(p,&i,&j,&k);
printf("Test %d %d %d %d\n",i,j,k,p);

if(!Memory_GetM(s->matrix_row,tSolver_Matrix_Element*,s->nelem*s->DOF*4) )
  return 0;

for(i=0;i<s->nelem*s->DOF*4;i++) s->matrix_row[i]=0;
for(i=0;i<3;i++)
  for(j=0;j<3;j++)
    Matrix_ElementAdd(s,i*j,i,j);

for(i=0;i<3;i++){
  me=s->matrix_row[i];
  printf("A[%d][.]=",me->p);
  while(me->Last_row){
    printf("%7.5g(%d)",me->a,me->q);
    me=me->Next;
  }
  printf("\n");
}
}

tSolver_Matrix_Element *Matrix_ElementAdd(tSolver_State *s, double a,
					  uint p, uint q){
tSolver_Matrix_Element *me;
  if(!s||!Memory_GetM(me,tSolver_Matrix_Element,1))
    return 0;
  me->a = a;
  me->p = p;
  me->q = q;


  me->Last_row=FALSE;
  me->Prev=0;
  me->Next=s->matrix_row[p];
  if(me->Next)
    me->Next->Prev=me;
  else
    me->Last_row=TRUE;
  s->matrix_row[p]=me;
  s->nummatrix++;

  printf("Matrix %d %d %7.5g %i %d \n",p,q,a,me->Last_row,s->nummatrix);
  
  return me;
} 

void Solver_Implicit(tSolver_State *s){

}
