/* See {SpreadEnergy.h} */
#include <SpreadEnergy.h>


#include <r4.h>
#include <Triangulation.h>
#include <Energy.h>

REVEAL
  T == Public BRANDED OBJECT
      ElemTableRec_t *top;     /* The map elements. */
      double K;                     /* Energy normalization constant */
      /* Defined by "defVar: */
      termVar: REF ARRAY OF bool_t;   /* Tells which terms are variable */
    OVERRIDES
      init = Init;
      defTop = DefTop;
      defVar = DefVar;
      eval = Eval;
      name = Name;
    }
    
SELF_Init(SELF_T *erg)
  { return erg;
  } /* END Init */

void SELF_DefTop(SELF_T *erg, ElemTableRec_t *top)
  { erg->top = top;
    erg->K = 1.0/((double)top->node.nel); 
    erg->termVar = NEW(REF ARRAY OF bool_t, top->node.nel);
    /* In case the client forgets to call "defVar": */
    for (i = 0 TO LAST(erg->termVar^)){ erg->termVar^[i] = FALSE; }
    
  } /* END DefTop */
  
void SELF_DefVar(SELF_T *erg, *variable: ARRAY OF bool_t)
  { /*
      There is a term for each existing node. The term of "v"
      is variable (termVar[v] == TRUE) if "v" is variable (variable[v] ==
      TRUE).
     
      A node contributes to the energy only if "termVar[v] == TRUE", i.e.
      if "v" exists and "v" is variable. 
    */
    int NV = erg->top->node.nel;
    ??? node = erg->top->node^;
    ??? termVar = erg->termVar^;
    {
      assert((variable.nel) == NV);
      for (v = 0; v < NV; v++)
      {
        termVar[v] = variable[v]) && (node[v]->exists
      }    
    }
  } /* END DefVar */
  
void SELF_Eval(
    SELF_T erg; 
    Coords_t *c; 
    double *e; 
    bool_t grad;
    Gradient *eDc;
 )
  { int NV = erg->top->node.nel;
    ??? termVar = erg->termVar^;
    ??? K = erg->K;
    {

      void AccumTerm(*cv: r4_t; VAR evDcv: r4_t)
        /*
          Adds to "e" the energy term corresponding to a node at "cv".
          Returns also the gradient "evDcv" of that term relative to "cv".
        */
        {
          e = e + K * r4_NormSqr(cv);
          if (grad)
          {
            evDcv = r4_Scale(2.0 * K, cv)
          }
          else
          {
            evDcv = (r4_t){0.0, 0.0, 0.0, 0.0};
          }
        } AccumTerm;

      {
        e = 0.0;
        for (v = 0; v < NV; v++)
        {
          if (termVar[v])
          {
            AccumTerm(c[v], eDc[v]);
          }
          else
          {
            eDc[v] = (r4_t){0.0, 0.0, 0.0, 0.0};
          }
        }
      }
    }
  } /* END Eval */

PROCEDURE Name(<*UNUSED); erg: T): char *==
  { return "Spread()";
  } /* END Name */

{
} SpreadEnergy.