qpatchmath-mini.cpp

算断裂的

 // ------------------------------------------------------------------
// qpatchmath-mini.cpp
//
// This file contains routines for math with patches including:
// converting parametric to real coordinates, computing directional
// derivatives, computing normals, and intersecting patches with
// lines.  In this version of the file, all the intersect routines
// are removed (to save space in mex files).
// ------------------------------------------------------------------
// Copyright (c) 1999 by Cornell University.  All rights reserved.
// 
// See the accompanying file 'Copyright' for authorship information,
// the terms of the license governing this software, and disclaimers
// concerning this software.
// ------------------------------------------------------------------
// This file is part of the QMG software.  
// Version 2.0 of QMG, release date September 3, 1999
// ------------------------------------------------------------------




#include "qpatchmath.h"
#include "qmatvec.h"



#ifdef DEBUGGING
namespace QMG {
  namespace Local_to_PatchMath_CPP {
    using namespace QMG;
    ostream* debug_str;
    bool dump_everything = false;
  }
}
#endif


// Types used throughout this file.


// Types passed to Lapack routines

#ifdef LAPACK_DOUBLECOMPLEX_T
typedef LAPACK_DOUBLECOMPLEX_T Lapack_doublecomplex;
#else
typedef double Lapack_doublecomplex;
#endif

#ifdef LAPACK_INT_T
typedef LAPACK_INT_T Lapack_int;
#else
typedef long int Lapack_int;
#endif

typedef std::complex<double> Complex;


// Helper class that makes private routines public for this file only.

class QMG::PatchMath::Bezier_Helper {
public:
  static Real control_point_coord(const PatchMath& pm, int cpnum, int d) {
    return pm.control_point_coord_(cpnum,d);
  }
};

// Class Workspace is a workspace for math operations.
// Used so that we don't have to keep calling 'new' to
// get workspace for Lapack.
// Used as follows: Set a marker.  Then allocate vectors
// and matrices.  When the marker is destroyed, the space
// is deallocated.


class QMG::PatchMath::Workspace::RealVector {
private:
  Workspace& wkspa_;
  int base_;
  int sz_;
  // no copying, no assignment.
  RealVector(const RealVector& o) : wkspa_(o.wkspa_)  { }
  void operator=(const RealVector&) { }
public:
  RealVector(Workspace& wkspa, int sz) : wkspa_(wkspa),
    base_(wkspa.rwkspa_.size()),
    sz_(sz) {
#ifdef DEBUGGING
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
    wkspa.rwkspa_.resize(base_ + sz_);
  }
  ~RealVector() { }
  Real& operator[](int i) {
#ifdef RANGECHECK
    if (i < 0 || i >= sz_)
      throw_error("Range err");
#endif
    return wkspa_.rwkspa_[base_ + i];
  }
  double* make_fortran_real() {
    return &(wkspa_.rwkspa_[base_]);
  }
};



// Matrices stored by column for fortran compatibility.

class QMG::PatchMath::Workspace::ComplexMatrix {
private:
  Workspace& wkspa_;
  int base_;
  int nr_;
  int nc_;
  // no copying, no assignment
  ComplexMatrix(const ComplexMatrix& o) : wkspa_(o.wkspa_) { }
  void operator=(const ComplexMatrix&) { }
public:
  ComplexMatrix(Workspace& wkspa, int nr, int nc) : wkspa_(wkspa),
    base_(wkspa.cwkspa_.size()),
    nr_(nr),
    nc_(nc) {
    wkspa.cwkspa_.resize(base_ + nr * nc);
#ifdef DEBUGGING
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
  }
  ~ComplexMatrix() { }
  Complex& operator()(int i, int j) {
#ifdef RANGECHECK
    if (i < 0 || i >= nr_ || j < 0 || j >= nc_)
      throw_error("Range err");
#endif
    return wkspa_.cwkspa_[base_ + i + j * nr_];
  }
  Lapack_doublecomplex* make_fortran_complex() {
    return reinterpret_cast<Lapack_doublecomplex*>(
      static_cast<Complex*>(&(wkspa_.cwkspa_[base_])));
  }
  int numrow() const { return nr_; }
  int numcol() const { return nc_; }
};



class QMG::PatchMath::Workspace::ComplexVector {
private:
  Workspace& wkspa_;
  int base_;
  int sz_;
  // no copying, no assignment
  ComplexVector(const ComplexVector& o) : wkspa_(o.wkspa_) { }
  void operator=(const ComplexVector&) { }
public:
  ComplexVector(Workspace& wkspa, int sz) : wkspa_(wkspa),
    base_(wkspa.cwkspa_.size()),
    sz_(sz) {
    wkspa.cwkspa_.resize(base_ + sz_, Complex(0.0,0.0));
#ifdef DEBUGGING
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
  }
  ~ComplexVector() { }
  Complex& operator[](int i) {
#ifdef RANGECHECK
    if (i < 0 || i >= sz_)
      throw_error("Range err");
#endif
    return wkspa_.cwkspa_[base_ + i];
  }
  Lapack_doublecomplex* make_fortran_complex() {
    return reinterpret_cast<Lapack_doublecomplex*>(
      static_cast<Complex*>(&(wkspa_.cwkspa_[base_])));
  }

};

class QMG::PatchMath::Workspace::RealMatrix {
private:
  Workspace& wkspa_;
  int base_;
  int nr_;
  int nc_;
  // no copying, no assignment
  RealMatrix(const RealMatrix& o) : wkspa_(o.wkspa_)  { }
  void operator=(const RealMatrix&) { }
public:
  RealMatrix(Workspace& wkspa, int nr, int nc) : wkspa_(wkspa),
    base_(wkspa.rwkspa_.size()),
    nr_(nr),
    nc_(nc) {
    wkspa.rwkspa_.resize(base_ + nr * nc);
#ifdef DEBUGGING
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
  }
  ~RealMatrix() { }
  Real& operator()(int i, int j) {
#ifdef RANGECHECK
    if (i < 0 || i >= nr_ || j < 0 || j >= nc_)
      throw_error("Range err");
#endif
    return wkspa_.rwkspa_[base_ + i + j * nr_];
  }
  Real* make_fortran_real() {
    return &(wkspa_.rwkspa_[base_]);
  }
  int numrow() const { return nr_; }
  int numcol() const { return nc_; }
};


class QMG::PatchMath::Workspace::StartPosMarker {
private:
  Workspace& wkspa_;
  int rpos_;
  int cpos_;
  StartPosMarker(const StartPosMarker& o) : wkspa_(o.wkspa_) { }
  void operator=(const StartPosMarker&) { }
public:
  explicit StartPosMarker(Workspace& wkspa) : wkspa_(wkspa),
    rpos_(wkspa.rwkspa_.size()),
    cpos_(wkspa.cwkspa_.size()) {
    if (wkspa.marker_active_)
      throw_error("startpos marker already active in workspace");
    wkspa.marker_active_ = true;
  }
  ~StartPosMarker() {
    wkspa_.rwkspa_.resize(rpos_);
    wkspa_.cwkspa_.resize(cpos_);
    wkspa_.marker_active_ = false;
  }
};



namespace {
  using namespace QMG;


        
  // This is a helper class for math on 
  // a one-variable Bezier.

  class Solve_OneVar_Bezier {
  public:
    virtual Real get_bez_coeff(int j) const = 0;
    Solve_OneVar_Bezier() { }
    virtual ~Solve_OneVar_Bezier() { }
    void solve(int degree,
      vector<PatchMath::IntersectRecord>& output_stack,
      PatchMath::Workspace& workspace, 
      double scfac, 
      double tol) const;
    Real eval(Real param,
      int degree) const;
    Real eval_deriv(Real param,
      int degree) const;
  };

  // Derived class for math on a one-variable Bezier that arises
  // from intersecting a line with a curve in 2D.

  class Solve_OneVar_Bezier1 : public Solve_OneVar_Bezier {
  private: 
    const Matrix& lineeq_;
    const Point& linerhs_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier1(const Matrix& lineeq,
      const Point& linerhs,
      const PatchMath& pm) : lineeq_(lineeq), linerhs_(linerhs), pm_(pm) { }
    Real get_bez_coeff(int cpnum) const {
      return lineeq_(0,0) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 0) +
        lineeq_(0,1) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 1) -
        linerhs_[0];
    }
  };

  // Derived class for math on a one-variable Bezier that
  // arises from intersecting a curve with a plane in 3D.

  class Solve_OneVar_Bezier2 : public Solve_OneVar_Bezier {
  private:
    const Point& plane_normal_;
    Real plane_rhs_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier2(const Point& plane_normal,
      Real plane_rhs,
      const PatchMath& pm) : plane_normal_(plane_normal), plane_rhs_(plane_rhs), 
      pm_(pm) { }
    Real get_bez_coeff(int cpnum) const {
      return plane_normal_[0] * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 0) +
        plane_normal_[1] * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum,1) +
        plane_normal_[2] * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum,2) -
        plane_rhs_;
    }
  };

  // Derived class for evaluating one dimension of a Bezier curve.

  class Solve_OneVar_Bezier3 : public Solve_OneVar_Bezier {
  private:
    int select_dim_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier3(int select_dim, 
      const PatchMath& pm) : select_dim_(select_dim), pm_(pm) 
    { }
    Real get_bez_coeff(int cpnum) const {
      return PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, select_dim_);
    }
  };


  // Derived class for evaluating the derivative of a Bezier curve.


  class Solve_OneVar_Bezier4 : public Solve_OneVar_Bezier {
  private:
    const Point& direc_;
    int di_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier4(const Point& direc, int di,
      const PatchMath& pm) : direc_(direc), di_(di), pm_(pm) 
    { }
    Real get_bez_coeff(int cpnum) const {
      Real t = 0.0;
      for (int j = 0; j < di_; ++j) {
        t += (PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum + 1, j) -
          PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, j)) *
          direc_[j];
      }
      return t;
    }
  };


  // Evaluate a one-variable Bezier curve using Horner algorithm.
  

  Real Solve_OneVar_Bezier::eval(Real param,
    int degree) const {
    Real u = 1.0 - param;
    Real pwrv = 1.0; 
    Real result = 0.0;
    int binom = 1;
    for (int j = 0; j < degree + 1; ++j) {
      result *= u;
      result += binom * pwrv * get_bez_coeff(j);      
      pwrv *= param;
      binom *= (degree - j);
      binom /= (j + 1);
    }
    return result;
  }
  
  // Evaluate derivative a one-variable Bezier curve using Horner algorithm.
  

  Real Solve_OneVar_Bezier::eval_deriv(Real param,
    int degree) const {
    Real u = 1.0 - param;
    Real pwrv = 1.0; 
    Real result = 0.0;
    int binom = 1;
    Real prevbez = get_bez_coeff(0);
    for (int j = 0; j < degree; ++j) {
      result *= u;
      Real nextbez = get_bez_coeff(j + 1);
      result += binom * pwrv * (nextbez - prevbez);
      prevbez = nextbez;
      pwrv *= param;
      binom *= (degree - j - 1);
      binom /= (j + 1);
    }
    return result * degree;
  }

  void throw_error1() {
    throw_error("QMG was not built correctly; qpatchmath-mini was linked to a mexfile instead of qpatchmath");
  }    

  // solve: find zeros of a one-variable Bezier.
  // Transforms to an eigenvalue problem and invokes LAPACK.
  // Doesn't set output_stack real coord fields -- caller must do this.
 
  void Solve_OneVar_Bezier::solve(int degree,
    vector<PatchMath::IntersectRecord>& output_stack,
    PatchMath::Workspace& workspace, 
    double scfac, 
    double tol) const {

    throw_error1();

  }

  // A point in R^2, used for representing Bezier parameters. 
  
  struct R2Point {
    Real coord[2];
    R2Point() { }
    R2Point(Real x, Real y) {
      coord[0] = x;
      coord[1] = y;
    }
    Real l2norm() const {
      return sqrt(coord[0] * coord[0] + coord[1] * coord[1]);