/* interaction.c */

#include "interaction.h"
#include "pools.h"
#include "error.h"
#include "basis.h"

#ifndef NOPOOLS

static POOL *interactionPool = (POOL *)NULL;
#define NEWINTERACTION() (INTERACTION *)NewPoolCell(sizeof(INTERACTION), 0, "galerkin interactions", &interactionPool)
#define DISPOSEINTERACTION(interaction) Dispose(interaction, &interactionPool)

static POOL *ffPool[MAXBASISSIZE*MAXBASISSIZE];	/* automatically initialized to NULL */
#define NEWFORMFACTORS(nr)	(float *)NewPoolCell(nr*sizeof(float), 0, "galerkin form factors", &ffPool[nr-1])
#define DISPOSEFORMFACTORS(ff, nr) Dispose(ff, &ffPool[nr-1])

#else /*NOPOOLS*/

#define NEWINTERACTION() (INTERACTION *)Alloc(sizeof(INTERACTION))
#define DISPOSEINTERACTION(interaction) Free((char *)interaction, sizeof(INTERACTION))

#define NEWFORMFACTORS(nr)	(float *)Alloc(nr * sizeof(float))
#define DISPOSEFORMFACTORS(ff, nr) Free((char *)ff, nr * sizeof(float))

#endif /*NOPOOLS*/

static int TotalInteractions = 0, 
  CCInteractions = 0, CSInteractions = 0, SCInteractions = 0, SSInteractions = 0
;

int GetNumberOfInteractions(void)
{
  return TotalInteractions;
}

int GetNumberOfClusterToClusterInteractions(void)
{
  return CCInteractions;
}

int GetNumberOfClusterToSurfaceInteractions(void)
{
  return CSInteractions;
}

int GetNumberOfSurfaceToClusterInteractions(void)
{
  return SCInteractions;
}

int GetNumberOfSurfaceToSurfaceInteractions(void)
{
  return SSInteractions;
}

INTERACTION *InteractionCreate(ELEMENT *rcv, ELEMENT *src,
			       FloatOrPointer K,
			       FloatOrPointer deltaK,
			       unsigned char nrcv, 
			       unsigned char nsrc,
			       unsigned char crcv,
			       unsigned char vis
			       /* #ident <<< WMP */
#ifdef WMP_WEIGHTS
			       , float *FK
#endif
			       /* #ident >>> WMP */
			       )
{
  INTERACTION *interaction;
  int i;

  interaction = NEWINTERACTION();
  interaction->rcv = rcv;
  interaction->src = src;
  interaction->nrcv = nrcv;
  interaction->nsrc = nsrc;
  interaction->crcv = crcv;
  interaction->vis = vis;

  if (nrcv == 1 && nsrc == 1)
    interaction->K.f = K.f;
  else {
    interaction->K.p = NEWFORMFACTORS(nrcv*nsrc);
    for (i=0; i<nrcv*nsrc; i++) interaction->K.p[i] = K.p[i];
  }

  /* #ident <<< WMP */
#ifdef WMP_WEIGHTS
  interaction->FK = NEWFORMFACTORS(nrcv);
  for (i=0; i<nrcv; i++) { interaction->FK[i] = FK[i]; }
#endif
  /* #ident >>> WMP */

  if (crcv>1) {
    Fatal(2, "InteractionCreate", "Not yet implemented for higher order approximations");
  }
  interaction->deltaK.f = deltaK.f;

  TotalInteractions++;
  if (IsCluster(rcv)) {
    if (IsCluster(src))
      CCInteractions++;
    else
      SCInteractions++;
  } else {
    if (IsCluster(src))
      CSInteractions++;
    else
      SSInteractions++;
  }

  return interaction;
}

INTERACTION *InteractionDuplicate(INTERACTION *interaction)
{
  return InteractionCreate(interaction->rcv, interaction->src,
			   interaction->K, interaction->deltaK,
			   interaction->nrcv, interaction->nsrc,
			   interaction->crcv, interaction->vis
			   /* #ident <<< WMP */
#ifdef WMP_WEIGHTS
			   , interaction->FK
#endif
			   /* #ident >>> WMP */
			   );
}

void InteractionDestroy(INTERACTION *interaction)
{
  ELEMENT *src=interaction->src, *rcv=interaction->rcv;

  if (interaction->nrcv>1 || interaction->nsrc>1) {
    DISPOSEFORMFACTORS(interaction->K.p, interaction->nrcv*interaction->nsrc);
  }

  /* #ident <<< WMP */
#ifdef WMP_WEIGHTS
  DISPOSEFORMFACTORS(interaction->FK, interaction->nrcv);
#endif
  /* #ident >>> WMP */

  DISPOSEINTERACTION(interaction);

  TotalInteractions--;
  if (IsCluster(rcv)) {
    if (IsCluster(src))
      CCInteractions--;
    else
      SCInteractions--;
  } else {
    if (IsCluster(src))
      CSInteractions--;
    else
      SSInteractions--;
  }
}

void InteractionPrint(FILE *out, INTERACTION *link)
{
  int a, b, c, alpha, beta, i;

  fprintf(out, "%d (", link->rcv->id);
  PrintElementId(out, link->rcv);
  fprintf(out, ") <- %d (", link->src->id);
  PrintElementId(out, link->src);
  fprintf(out, "): vis=%d, nrcv=%d, nsrc=%d, crcv=%d,\n",
	  link->vis,
	  link->nrcv,
	  link->nsrc,
	  link->crcv);
  if (link->vis==0)
    return;

  a = link->nrcv;
  b = link->nsrc;
  if (a==1 && b==1) 
    fprintf(out, "\tK = %f\n", link->K.f);
  else {
    fprintf(out, "\tK[%d][%d] = \n", a, b);
    for (alpha=0; alpha<a; alpha++) {
      fprintf(out, "\t\t");
      for (beta=0; beta<b; beta++) {
	fprintf(out, "%7.7f\t", link->K.p[alpha*b + beta]);
      }
      fprintf(out, "\n");
    }
  }
    
  c = link->crcv;
  if (c==1)
    fprintf(out, "\tdeltaK = %f\n", link->deltaK.f);
  else {
    fprintf(out, "\tdeltaK[%d] = ", c);
    for (i=0; i<c; i++)
      fprintf(out, "%7.7f ", link->deltaK.p[i]);
    fprintf(out, "\n");
  }
}