qpatchmath-mini.cpp

算断裂的

    }
  };

  // Class for math with bezier triangles.


  class Solve_BezTri {
  public:
    virtual Real get_bez_coeff(int eqnum, int j) const = 0;
    Solve_BezTri() { }
    virtual ~Solve_BezTri() { }



    bool solve(int degree,
      vector<PatchMath::IntersectRecord>& output_stack,
      PatchMath::Workspace& workspace,
      double scfac, 
      double tol,
      bool newton_improve) const;
    Real eval(const R2Point& param,
      int eqnum,
      int degree) const;
    Real eval_direc_deriv(const R2Point& param,
      const R2Point& direc,
      int eqnum,
      int degree) const;
  };

  // Derived class for intersecting a line with a triangle.

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

  // Derived class for evaluating directional deriv.

  class Solve_BezTri2 : public Solve_BezTri {
  private:
    const Point& direc_;
    int di_;
    const PatchMath& pm_;
  public:
    Solve_BezTri2(const Point& direc, 
      int di,
      const PatchMath& pm) : direc_(direc), di_(di), pm_(pm) 
    { }
    Real get_bez_coeff(int eqnum, int cpnum) const;
  };
  
  

  // Derived class for evaluating one coordinate of a bezier triangle.

  class Solve_BezTri3 : public Solve_BezTri {
  private:
    int select_dim_;
    const PatchMath& pm_;
  public:
    Solve_BezTri3(int select_dim, 
      const PatchMath& pm) : select_dim_(select_dim), pm_(pm) 
    { }
    Real get_bez_coeff(int, int cpnum) const {
      return PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, select_dim_);
    }
  };
  
  
  Real Solve_BezTri2::get_bez_coeff(int eqnum, int cpnum) const {
    Real t = 0.0;
    int offset = static_cast<int>(sqrt(2 * static_cast<double>(cpnum)));
    if (offset * (offset + 1) / 2 > cpnum)
      --offset;
    int bcp = cpnum + offset + 1;
    int ecp = (eqnum == 0)? bcp + 1 : cpnum;
    for (int j = 0; j < di_; ++j) {
      t += (PatchMath::Bezier_Helper::control_point_coord(pm_, ecp, j) -
        PatchMath::Bezier_Helper::control_point_coord(pm_, bcp, j)) *
        direc_[j];
    }
    return t;
  }
    
  // Compute real coord from param coord using Horner's algorithm.

  Real Solve_BezTri::eval(const R2Point& param,
    int eqnum,
    int degree) const {
    Real p0 = param.coord[0];
    Real p1 = param.coord[1];
    Real p2 = 1.0 - p0 - p1;
    int outer_binom = 1;
    Real pwrv = 1.0;
    Real result = 0.0;
    for (int row = degree; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      int inner_binom = 1;
      Real pwru = 1.0;
      Real row_result = 0.0;
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += inner_binom * 
          get_bez_coeff(eqnum, base_idx + col) * pwru;
        pwru *= p0;
        inner_binom *= (row - col);
        inner_binom /= (col + 1);
      }
      result += row_result * pwrv * outer_binom;
      outer_binom *= row;
      outer_binom /= (degree - row + 1);
      pwrv *= p1;
    }
    return result;
  }

    
  // Compute directional derivative using Horner's rule.
  
  Real Solve_BezTri::eval_direc_deriv(const R2Point& param,
    const R2Point& direc,
    int eqnum,
    int degree) const {   
    Real p0 = param.coord[0];
    Real p1 = param.coord[1];
    Real p2 = 1.0 - p0 - p1;
    Real d0 = direc.coord[0];
    Real d1 = direc.coord[1];
    Real d2 = -d0 - d1;    
    int outer_binom = 1;
    Real pwrv = 1.0;
    Real result = 0.0;
    for (int row = degree - 1; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      int inner_binom = 1;
      Real pwru = 1.0;
      Real row_result = 0.0;
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += inner_binom * 
          (d1 * get_bez_coeff(eqnum, base_idx + col) +
          d2 * get_bez_coeff(eqnum, base_idx + col + row + 1) +
          d0 * get_bez_coeff(eqnum, base_idx + col + row + 2))
          * pwru;
        pwru *= p0;
        inner_binom *= (row - col);
        inner_binom /= (col + 1);
      }
      result += row_result * pwrv * outer_binom;
      outer_binom *= row;
      outer_binom /= (degree - row);
      pwrv *= p1;
    }
    return degree * result;
  }



  // Solve a system of two polynomial equations
  // in Bezier triangle form.  Use Macaulay-type
  // result and random u-direction.  Reduce
  // to small system and invoke LAPACK for
  // generalized eigenvalues.
  // Returns true if there was a global degeneracy
  // during the computation.

  
  bool Solve_BezTri::solve(int degree,
    vector<PatchMath::IntersectRecord>& output_stack,
    PatchMath::Workspace& workspace, 
    double scfac, 
    double tol,
    bool newton_improve) const {

    throw_error1();
    return false;
  }

  // Helper class for math with bezier quad patches.

  class Solve_BezQuad {
  public:
    virtual Real get_bez_coeff(int eqnum, int j) const = 0;
    Solve_BezQuad() { }
    virtual ~Solve_BezQuad() { }
    bool solve(int degree1[2],
      int degree2[2],
      vector<PatchMath::IntersectRecord>& output_stack,
      PatchMath::Workspace& workspace, 
      double scfac, 
      double tol,
      bool newton_improve) const;
    Real eval(const R2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Real eval_partialu(const R2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Real eval_partialv(const R2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
  };

  // Derived class for intersecting a patch with a line.
  
  class Solve_BezQuad1 : public Solve_BezQuad {
  private: 
    const Matrix& lineeq_;
    const Point& linerhs_;
    const PatchMath& pm_;
  public:
    Solve_BezQuad1(const Matrix& lineeq,
      const Point& linerhs,
      const PatchMath& pm) : lineeq_(lineeq), linerhs_(linerhs), pm_(pm) { }
    Real get_bez_coeff(int eqnum, int cpnum) const {
      return lineeq_(eqnum,0) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 0) +
        lineeq_(eqnum,1) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 1) +
        lineeq_(eqnum,2) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 2) -
        linerhs_[eqnum];
    }
  };

  // derived class for directional derivative.

  class Solve_BezQuad2 : public Solve_BezQuad {
  private:
    const Point& direc_;
    int di_;
    const PatchMath& pm_;
    int degree1_;
    int degree2_;
  public:
    Solve_BezQuad2(const Point& direc, 
      int di,
      const PatchMath& pm,
      int degree1,
      int degree2) : direc_(direc), di_(di), pm_(pm), degree1_(degree1), degree2_(degree2)
    { }
    Real get_bez_coeff(int eqnum, int cpnum) const;
  };

  // derived class for evaluating one dimension.
    
  class Solve_BezQuad3 : public Solve_BezQuad {
  private:
    int select_dim_;
    const PatchMath& pm_;
  public:
    Solve_BezQuad3(int select_dim, 
      const PatchMath& pm) : select_dim_(select_dim), pm_(pm) 
    { }
    Real get_bez_coeff(int, int cpnum) const {
      return PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, select_dim_);
    }
  };

    
  Real Solve_BezQuad2::get_bez_coeff(int eqnum, int cpnum) const {
    Real t = 0.0;
    int ocp1, ocp2, fac;
    if (eqnum == 0) {
      int d2 = cpnum / degree1_;
      int d1 = cpnum % degree1_;
      ocp1 = d1 + 1 + d2 * (degree1_ + 1);
      ocp2 = d1 + d2 * (degree1_ + 1);
      fac = degree1_;
    }
    else {
      int d2 = cpnum / (degree1_ + 1);
      int d1 = cpnum % (degree1_ + 1);
      ocp1 = d1 + (d2 + 1) * (degree1_ + 1);
      ocp2 = d1 + d2 * (degree1_ + 1);
      fac = degree2_;
    }
    for (int j = 0; j < di_; ++j) {
        t += (PatchMath::Bezier_Helper::control_point_coord(pm_, ocp1, j) -
          PatchMath::Bezier_Helper::control_point_coord(pm_, ocp2, j)) * fac 
          * direc_[j];
    }
    return t;
  }

  // Horner algorithm for evaluating bezier quad (parametric
  // to real conversion)


  Real Solve_BezQuad::eval(const R2Point& param,
    int eqnum,
    int degree1,
    int degree2) const {

    Real u = param.coord[0];
    Real oppu = 1.0 - u;
    Real v = param.coord[1];
    Real oppv = 1.0 - v;

    int outer_binom = 1;
    Real outer_result = 0.0;
    Real pwrv = 1.0;
    for (int d2 = 0; d2 <= degree2; ++d2) {
      // Compute the inner coefficient.

      int inner_binom = 1;
      Real pwru = 1.0;
      Real inner_result = 0.0;
      for (int d1 = 0; d1 <= degree1; ++d1) {
        inner_result *= oppu;
        inner_result += get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1) *
          inner_binom * pwru;
        pwru *= u;
        inner_binom *= (degree1 - d1);
        inner_binom /= (d1 + 1);
      }
      outer_result *= oppv;
      outer_result += inner_result * outer_binom * pwrv;
      pwrv *= v;
      outer_binom *= (degree2 - d2);
      outer_binom /= (d2 + 1);
    }
    return outer_result;
  }

  
  
  
  Real Solve_BezQuad::eval_partialu(const R2Point& param,
    int eqnum,
    int degree1,
    int degree2) const { 

    Real u = param.coord[0];
    Real oppu = 1.0 - u;
    Real v = param.coord[1];
    Real oppv = 1.0 - v;
    int outer_binom = 1;
    Real outer_result = 0.0;
    Real pwrv = 1.0;
    for (int d2 = 0; d2 <= degree2; ++d2) {
      // Compute the inner coefficient.

      int inner_binom = 1;
      Real pwru = 1.0;
      Real inner_result = 0.0;
      Real prev_coef = get_bez_coeff(eqnum, d2 * (degree1 + 1));
      for (int d1 = 0; d1 < degree1; ++d1) {
        inner_result *= oppu;
        Real next_coef = get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1 + 1);
        inner_result += (next_coef - prev_coef) * inner_binom * pwru;
        prev_coef = next_coef;
        pwru *= u;
        inner_binom *= (degree1 - d1 - 1);
        inner_binom /= (d1 + 1);
      }
      outer_result *= oppv;
      outer_result += inner_result * outer_binom * pwrv;
      pwrv *= v;
      outer_binom *= (degree2 - d2);
      outer_binom /= (d2 + 1);
    }
    return outer_result * degree1;
  }



  Real Solve_BezQuad::eval_partialv(const R2Point& param,
    int eqnum,
    int degree1,
    int degree2) const { 

    Real u = param.coord[0];
    Real oppu = 1.0 - u;
    Real v = param.coord[1];
    Real oppv = 1.0 - v;
    int outer_binom = 1;
    Real outer_result = 0.0;
    Real pwrv = 1.0;
    for (int d2 = 0; d2 < degree2; ++d2) {
      // Compute the inner coefficient.

      int inner_binom = 1;
      Real pwru = 1.0;
      Real inner_result = 0.0;
      for (int d1 = 0; d1 <= degree1; ++d1) {
        inner_result *= oppu;
        Real coef =  get_bez_coeff(eqnum, (d2 + 1) * (degree1 + 1) + d1) -
          get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1);
        inner_result += coef * inner_binom * pwru;
        pwru *= u;
        inner_binom *= (degree1 - d1);
        inner_binom /= (d1 + 1);
      }
      outer_result *= oppv;
      outer_result += inner_result * outer_binom * pwrv;
      pwrv *= v;
      outer_binom *= (degree2 - d2 - 1);
      outer_binom /= (d2 + 1);
    }
    return outer_result * degree2;
  }