#ifndef MESHPATCH
#define MESHPATCH

#include <point.h>
#include <print.h>

#include <bernstein.h>

typedef point3<double> pt3;

/*!
\brief Rectangular bezier patch.

This is a generalization of a patch. The client evaluates
 patches with a functional object pij in eval.
*/
class meshpatch
{
public:

  /** The number of points on the u-axis. */
  uintc M;
  /** The number of points of the v-axis. */
  uintc N;

protected:
  bernsteinPoly * const vbu;
  bernsteinPoly * const vbv;

  double * const vbuval;
  double * const vbvval;

public:

  /** Build a mesh patch with bernstein polynomials. */
  meshpatch(uintc M_, uintc N_);

  /** Memory cleanup. */
  ~meshpatch();

  /** Evaluate an actual patch. */
  template< typename T, typename W >
  void eval
  ( 
    T & p, 
    double const u, 
    double const v,
    W const & pij 
  ) const
  {
    for (uint i=0; i<=M; ++i)
      vbuval[i] = vbu[i].eval(u);

    for (uint i=0; i<=N; ++i)
      vbvval[i] = vbv[i].eval(v);

    p = T();
    T p2;
    uint k;
    for (uint i=0; i<=M; ++i)
    {
      T p2;
      for (k=0; k<=N; ++k)
        p2 += pij(i,k)*vbvval[k];

      p += p2*vbuval[i];
    }
  }

};



#endif