// Gnurbs - A curve and surface library
// Copyright (C) 2008-2026 Eric Bechet
//
// See the LICENSE file for contributions and license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//

#ifndef __NBSPLINESURFACE_H
#define __NBSPLINESURFACE_H


#include "nsurface.h"
#include "nbspline.h"
#include <vector>

class nbsplinesurface : public nparametricsurface
{
protected:
  std::vector<npoint> val; //!< Control points (ordered in u then in v).
  int nCPu; //!< Number of control points (u).
  int nCPv; //!< Number of control points (v).
  std::vector<double> knotsu; //!< Nodal sequence U.
  std::vector<double> knotsv; //!< Nodal sequence V.
  int nknotsu; //!< Number of knots (U).
  int nknotsv; //!< Number of knots (V).
  int degu; //!< Degree in u.
  int degv; //!< Degree in v.

public:

  /// \brief default constructor.
  /// \warning Control points and knots are not initialized.
  nbsplinesurface() : val(1),nCPu(1),degu(0),nCPv(1),degv(0),nknotsu(2),nknotsv(2),knotsu(2),knotsv(2)
  {
    knotsu[0]=0;
    knotsu[1]=1;
    knotsv[0]=0;
    knotsv[1]=1;
    val[0]=npoint();
  }

  /// \brief Constructor.
  /// \param[in] nCPu_ number of control points in u.
  /// \param[in] degreeu_ surface degree in u.
  /// \param[in] nCPu_ number of control points in v.
  /// \param[in] degreeu_ surface degree in v.
  /// \warning Control points and knots are not initialized.
  nbsplinesurface(int nCPu_,int degreeu_,int nCPv_,int degreev_) : val(nCPu_*nCPv_),
      nCPu(nCPu_),degu(degreeu_),nknotsu(nCPu_+degreeu_+1),knotsu(nCPu_+degreeu_+1),
      nCPv(nCPv_),degv(degreev_),nknotsv(nCPv_+degreev_+1),knotsv(nCPv_+degreev_+1) {}

  /// \brief Constructor.
  /// \param[in] nCPu_ number of control points in u.
  /// \param[in] degreeu_ surface degree in u.
  /// \param[in] nCPu_ number of control points in v.
  /// \param[in] degreeu_ surface degree in v.
  /// \param[in] knotsu_ knots vector U.
  /// \param[in] knotsv_ knots vector V.
  /// \param[in] pts_ control points vector (fastest varying index corresponds to U).
  nbsplinesurface(int nCPu_,int degreeu_,int nCPv_,int degreev_,double knotsu_[],double knotsv_[], npoint pts_[]) : val(nCPu_*nCPv_),
      nCPu(nCPu_),degu(degreeu_),nknotsu(nCPu_+degreeu_+1),knotsu(nCPu_+degreeu_+1),
      nCPv(nCPv_),degv(degreev_),nknotsv(nCPv_+degreev_+1),knotsv(nCPv_+degreev_+1) 
      {
        for (int i=0;i<nknotsu;++i) knotsu[i]=knotsu_[i];
        for (int i=0;i<nknotsv;++i) knotsv[i]=knotsv_[i];
        for (int i=0;i<nCPu*nCPv;++i) val[i]=pts_[i];
      }
  /// \brief Destructor.
  virtual ~nbsplinesurface() {}

  /// \brief Reserve memory data for the surface.
  /// \warning Control points and knots are not initialized.
  void init(int nCPu_,int degreeu_,int nCPv_,int degreev_)
  {
      val.resize(nCPu_*nCPv_);
      nCPu=nCPu_;
      degu=degreeu_;
      nknotsu=nCPu_+degreeu_+1;
      knotsu.resize(nCPu_+degreeu_+1);
      nCPv=nCPv_;
      degv=degreev_;
      nknotsv=nCPv_+degreev_+1;
      knotsv.resize(nCPv_+degreev_+1);
  }

  /// \brief Construct a Coons patch.
  /// \param[in] Cu0 first curve in u.
  /// \param[in] Cu1 second curve in u.
  /// \param[in] Cv0 first curve in v.
  /// \param[in] Cv1 second curve in v.
  /// \param[out] data data to fill-in.
  void makecoons(const nbspline &Cu0,const nbspline &Cu1,const nbspline &Cv0,const nbspline &Cv1,nbsplinesurface &S1,nbsplinesurface &S2,nbsplinesurface &S3);
  virtual int degree(void) const;
  virtual int degree_u(void) const
  {
    return degu;
  }
  virtual int degree_v(void) const
  {
    return degv;
  }
  virtual int nb_CP(void) const
  {
    return nCPu*nCPv;
  }
  virtual int nb_CP_u(int iv=0) const
  {
    return nCPu;
  }
  virtual int nb_CP_v(int iu=0) const
  {
    return nCPv;
  }

  /// \brief Number of knots in u.
  /// \return Number of knots.
  virtual int nb_knots_u(void) const
  {
    return nknotsu ;
  }

  /// \brief Number of knots in u.
  /// \return Number of knots.
  virtual int nb_knots_v(void) const
  {
    return nknotsv ;
  }
  virtual double min_u() const
  {
    return knotsu[degu];
  }
  virtual double max_u() const
  {
    return knotsu[nknotsu-degu-1];
  }
  virtual double min_v() const
  {
    return knotsv[degv];
  }
  virtual double max_v() const
  {
    return knotsv[nknotsv-degv-1];
  }
  virtual void P(double u,double v, npoint& ret) const ;
  
  
  virtual void dPdu(double u,double v, npoint& ret) const ;
  virtual void dPdv(double u,double v, npoint& ret) const ;
  
  virtual void d2Pdu2(double u,double v, npoint& ret) const ;
  virtual void d2Pdv2(double u,double v, npoint& ret) const ;
  
  virtual void d2Pdudv(double u,double v, npoint& ret) const ;
  
  
  /// \brief Gets control point (const version).
  /// \param[in] whichu control point number in u.
  /// \param[in] whichv control point number in v.
  /// \return Control point.
  virtual npoint CP(int whichu,int whichv) const
  {
    return val[whichu+nCPu*whichv];
  }
  virtual npoint CP(int which) const
  {
    return val[which];
  }

  /// \brief Gets control point (non-const version).
  /// \param[in] whichu control point number in u.
  /// \param[in] whichv control point number in v.
  /// \return Control point.
  virtual npoint& CP(int whichu,int whichv)
  {
    return val[whichu+nCPu*whichv];
  }

  /// \brief Gets control point (non-const version).
  /// \param[in] which control point number.
  /// \return Control point.
  virtual npoint& CP(int which)
  {
    return val[which];
  }

  /// \brief Sets control point.
  /// \param[in] whichu control point number in u.
  /// \param[in] whichv control point number in v.
  /// \param[in] pt new control point.
  virtual void set_CP(int whichu,int whichv,const npoint& pt)
  {
    val[whichu+nCPu*whichv]=pt;
  }
  virtual void set_CP(int which,const npoint& pt)
  {
    val[which]=pt;
  }

  /// \brief Gets a knot of U (const version).
  /// \param[in] which wanted knot number.
  /// \return Knot value.
  virtual double u(int which) const
  {
    return knotsu[which];
  }

  /// \brief Gets a knot of U (non-const version).
  /// \param[in] which wanted knot number.
  /// \return Knot value.
  virtual double& u(int which)
  {
    return knotsu[which];
  }

  /// \brief Gets a knot of V (const version).
  /// \param[in] which wanted knot number.
  /// \return Knot value.
  virtual double v(int which) const
  {
    return knotsv[which];
  }

  /// \brief Gets a knot of V (non-const version).
  /// \param[in] which wanted knot number.
  /// \return Knot value.
  virtual double& v(int which)
  {
    return knotsv[which];
  }

  /// \brief Inserts a knot r times in U.
  /// \param[in] u knot value.
  /// \param[in] r number of insertions.
  void insertknot_u(double u,int r);

  /// \brief Saturates nodal sequence U.
  void saturate_u();

  /// \brief Elevates surface degree by 1 in u.
  void degree_elevation_u();

  /// \brief Inserts a knot r times in V.
  /// \param[in] v knot value.
  /// \param[in] r number of insertions.
  void insertknot_v(double v,int r);

  /// \brief Saturates nodal sequence V.
  void saturate_v();

  /// \brief Elevates surface degree by 1 in v.
  void degree_elevation_v();
  
  int findspan(const std::vector<double> &knots,int nknots, int deg, double u) const ;
  void basisfuns(const std::vector<double> &knots,int nknots, int deg, int i, double para, std::vector<double> &N) const ;
  void gradbasisfuns(const std::vector<double> &knots,int nknots,int deg,int i, double u,int n,std::vector<std::vector<double> > &N) const ;
  double Nip(int p,int m,std::vector<double> &knots, int i, double u);
  /// \brief Compute ONE shape functions and its first derivatives at once. Takes weighting factor into account.
  void SF(double u, double v,int k,double &SF,double &dSFdu,double &dSFdv) const;
  virtual nsurface* clone() const
  {
    return new nbsplinesurface(*this);
  }
private:
  void insertknotu_(double u,int k,int s,int r);
  void insertknotv_(double v,int k,int s,int r);
  void findspanmult(const std::vector<double> &knots,int nknots, int deg,double u,int &k,int &s);
  void curvepoint(double u, double v, npoint &C) const ;
  void dercurvepoint(double u, double v, int nu,int nv,npoint &C) const ;
  
};



#endif // __NBSPLINESURFACE_H