qpatchmath-mini.cpp

算断裂的

  // Routine to solve a 2-var polynomial system in 
  // tensor product Bezier form.  Use Macaulay
  // resultant combined with u-resultant with randomized direction.
  // Return true if degenerate.
  
  bool Solve_BezQuad::solve(int udegree[2],
    int vdegree[2],
    vector<PatchMath::IntersectRecord>& output_stack,
    PatchMath::Workspace& workspace,
    double scfac, 
    double tol,
    bool newton_improve) const {

    throw_error1();
    return false;
  }
}


// Intersect a segment with a curve (2D) or patch (3D).

void QMG::PatchMath::intersect_seg(vector<IntersectRecord>& output_stack,
                                   Workspace& workspace,
                                   const Point& startpoint,
                                   const Point& endpoint,
                                   double scfac,
                                   double tol,
                                   bool newton_improve) const {
  
#ifdef DEBUGGING
  using Local_to_PatchMath_CPP::dump_everything;
  using Local_to_PatchMath_CPP::debug_str;
  if (dump_everything)
    *debug_str << " reached intersect_seg\n";
#endif

  if (gdim_ != di_ - 1)
    return;

  Point tangent = Point::subtract(endpoint,startpoint);
  Real lth = tangent.l2norm();
  if (lth == 0.0)
    throw_error("Zero-length segment in intersect_seg");


  // Convert the segment from implicit to parametric form.

  Real lineeqspace[MAXDIM * MAXDIM];
  Matrix lineeq(di_ - 1, di_, lineeqspace);

  if (di_ == 2) {
    lineeq(0,0) = tangent[1] / lth;
    lineeq(0,1) = -tangent[0] / lth;
  }
  else {
    Point hhtransform;
    Real beta = 
      Matrix::compute_hh_transform(&hhtransform[0], &tangent[0], 1, 3);
    for (int i = 0; i < 2; ++i) {
      Real mult = -beta * hhtransform[i + 1];
      for (int j = 0; j < 3; ++j) {
        lineeq(i,j) = mult * hhtransform[j];
      }
      lineeq(i, i + 1) += 1.0;
    }
  }

  Point linerhs;
  {
    for (int i = 0; i < di_ - 1; ++i) {
      Real t = 0.0;
      for (int j = 0; j < di_; ++j) {
        t += lineeq(i,j) * startpoint[j];
      }
      linerhs[i] = t;
    }
  }
  int init_stacksize = output_stack.size();
  intersect_line(output_stack, workspace, lineeq, linerhs, scfac, tol, newton_improve);
  int end_stacksize = output_stack.size();

  int output_pos = init_stacksize;
  int degree_fac;
  if (gdim_ == 1 || ptype_ == BEZIER_TRIANGLE) {
    degree_fac = (degree1_ + 1) * (degree1_ + 1) * 2;
  }
  else {
    degree_fac = (degree1_ + 1) * (degree2_ + 1) * 2;
  }



  Real ip0 = Point::inner_product(startpoint, tangent);
  Real ip1 = lth * lth;

  for (int input_pos = init_stacksize; input_pos < end_stacksize; ++input_pos) {
    Real ip2 = Point::inner_product(output_stack[input_pos].realcoord, tangent);
    Real t = (ip2 - ip0) / ip1;

    Real seg_unc = output_stack[input_pos].param_uncertainty * degree_fac *
      scfac / lth;
    
#ifdef DEBUGGING
    if (dump_everything) {
      *debug_str << "input_pos = " << input_pos << " output_pos = " << output_pos
        << " ips = " << ip0 << " " << ip1 << " " << ip2
        << "\n t = " << std::setprecision(17) << t <<  " segunc = " << seg_unc << '\n';
    }
#endif
    if (t > -seg_unc && t < 1.0 + seg_unc) {
      output_stack[output_pos] = output_stack[input_pos];
      output_stack[output_pos].dist = t;
      if (fabs(t) <= seg_unc || fabs(t - 1.0) <= seg_unc) {
        output_stack[output_pos].degen = true;
      }
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "keeping, degen = " << output_stack[output_pos].degen << '\n';
      }
#endif
      output_stack[output_pos].param_uncertainty_seg = seg_unc;
      ++output_pos;
    }
  }
  output_stack.resize(output_pos, IntersectRecord());
}


// Intersect a curve (in 2D) or patch (in 3D) with a line.


void QMG::PatchMath::intersect_line(vector<IntersectRecord>& output_stack,
                                    Workspace& workspace,
                                    const Matrix& lineeq,
                                    const Point& linerhs,
                                    double scfac,
                                    double tol,
                                    bool newton_improve) const {

  if (gdim_ != di_ - 1)
    return;

#ifdef DEBUGGING
  using Local_to_PatchMath_CPP::dump_everything;
  using Local_to_PatchMath_CPP::debug_str;
  if (dump_everything)
    *debug_str << " reached intersect_line\n";
#endif

  int initsize = output_stack.size();

  if (di_ == 2) {
    Solve_OneVar_Bezier1 helper(lineeq, linerhs, *this);
    helper.solve(degree1_, output_stack, workspace, scfac, tol);
  }
  else {
    if (ptype_ ==  BEZIER_TRIANGLE) {
      Solve_BezTri1 helper(lineeq, linerhs, *this);
      helper.solve(degree1_, output_stack, workspace, scfac, tol, newton_improve);
    }
    else {
      Solve_BezQuad1 helper(lineeq, linerhs, *this);
      int deg1[2];
      int deg2[2];
      deg1[0] = deg1[1] = degree1_;
      deg2[0] = deg2[1] = degree2_;
      helper.solve(deg1, deg2, output_stack, workspace, scfac, tol, newton_improve);
    }
  }
  for (int j = initsize; j < output_stack.size(); ++j) {
    output_stack[j].realcoord = 
      get_real_point(output_stack[j].paramcoord);
#ifdef DEBUGGING
    if (dump_everything)
      *debug_str << " real point = " << output_stack[j].realcoord << '\n';
#endif
    
  }
  return;
}


// Intersect a plane with a bezier curve.


void QMG::PatchMath::intersect_plane(vector<IntersectRecord>& output_stack,
                                     Workspace& workspace,
                                     const Point& plane_normal,
                                     Real rhs,
                                     double scfac,
                                     double tol,
                                     bool newton_improve) const {

#ifdef DEBUGGING
  if (di_ != 3 || gdim_ != 1)
    throw_error("Intersect_plane called on wrong dimensions");
#endif

  int initsize = output_stack.size();
  Solve_OneVar_Bezier2 helper(plane_normal, rhs, *this);
  helper.solve(degree1_, output_stack, workspace, scfac, tol);
  for (int j = initsize; j < output_stack.size(); ++j) {
    output_stack[j].realcoord = 
      get_real_point(output_stack[j].paramcoord);
  }
}

// Get extreme points of a Bezier patch or curve.
// This means extreme in a certain direc direc.
// Do this by finding roots of the directional derivative.


void QMG::PatchMath::get_extreme_points(vector<IntersectRecord>& output_stack,
                                        Workspace& workspace,
                                        const Point& direc,
                                        double scfac,
                                        double tol) const {
  int initsize = output_stack.size();
  if (gdim_ == 0) {
    IntersectRecord result;
    for (int j = 0; j < di_; ++j) 
      result.realcoord[j] = control_point_coord_(0,j);
    result.degen = false;
    result.param_uncertainty = tol;
    result.param_uncertainty_seg = 0.0;
    result.dist = 0.0;
    output_stack.push_back(result);
  }
  else {
    if (degree1_ == 1) return;
    if (gdim_ == 1) {
      Solve_OneVar_Bezier4 helper(direc, di_, *this);
      helper.solve(degree1_ - 1, output_stack, workspace, scfac, tol);
    }
    else { // get here if gdim == 2, di = 3.
      if (ptype_ == BEZIER_TRIANGLE) {
        Solve_BezTri2 helper(direc, di_, *this);
        helper.solve(degree1_ - 1, output_stack, workspace, scfac, tol, true);
      }
      else {
        if (degree2_ == 1) return;
        Solve_BezQuad2 helper(direc, di_, *this, degree1_, degree2_);
        int deg1[2];
        int deg2[2];
        deg1[0] = degree1_ - 1;
        deg2[0] = degree2_;
        deg1[1] = degree1_;
        deg2[1] = degree2_ - 1;
        helper.solve(deg1, deg2, output_stack, workspace, scfac, tol, true);
      }
    }
    for (int k = initsize; k < output_stack.size(); ++k) {
      output_stack[k].realcoord = 
        get_real_point(output_stack[k].paramcoord);
    }
  }
}
      

QMG::Point 
QMG::PatchMath::get_real_point(const Point& param_coord) const {
  Point real_coord;

  if (gdim_ == 0) {
    for (int j = 0; j < di_; ++j)
      real_coord[j] = control_point_coord_(0,j);  
  }
  else if (gdim_ == 1) {
    for (int k = 0; k < di_; ++k) {
      Solve_OneVar_Bezier3 e1(k, *this);
      real_coord[k] = e1.eval(param_coord[0], degree1_);
    }
  }
  else  {//  gdim == 2

    if (ptype_ == BEZIER_TRIANGLE) {
      
      // Use deCasteljau algorithm for a Bezier triangle

      for (int k = 0; k < di_; ++k) {
        Solve_BezTri3 e1(k, *this);
        real_coord[k] = e1.eval(R2Point(param_coord[0], param_coord[1]),
          0, degree1_);
      }
    }
    else {
      for (int k = 0; k < di_; ++k) {
        Solve_BezQuad3 e1(k, *this);
        real_coord[k] = e1.eval(R2Point(param_coord[0], param_coord[1]),
          0, degree1_, degree2_);
      }
    }
  }
  return real_coord;
}

QMG::Point 
QMG::PatchMath::direc_deriv(const Point& param_coord, 
                            const Point& direc) const {

  Point real_coord;

  if (gdim_ == 0) {
    throw_error("Directional derivative of 0D-patch not allowed.");
  }
  else if (gdim_ == 1) {
    for (int k = 0; k < di_; ++k) {
      Solve_OneVar_Bezier3 e1(k, *this);
      real_coord[k] = direc[0] * e1.eval_deriv(param_coord[0], 
        degree1_);
    }
  }
  else if (gdim_ == 2) {

    if (ptype_ == BEZIER_TRIANGLE) {
      
      for (int k = 0; k < di_; ++k) {
        Solve_BezTri3 e1(k, *this);
        real_coord[k] = e1.eval_direc_deriv(R2Point(param_coord[0], param_coord[1]),
          R2Point(direc[0], direc[1]), 0, degree1_);
      }
    }
    else {
      Real partial;
      for (int k = 0; k < di_; ++k) {
        real_coord[k] = 0.0;
        Solve_BezQuad3 e1(k, *this);
        if (direc[0] != 0) {
          partial = e1.eval_partialu(R2Point(param_coord[0], param_coord[1]),
            0, degree1_, degree2_);
          real_coord[k] += partial * direc[0];
        }
        if (direc[1] != 0) {
          partial = e1.eval_partialv(R2Point(param_coord[0], param_coord[1]),
            0, degree1_, degree2_);
          real_coord[k] += partial * direc[1];
        }
      }
    }
  }
  return real_coord;
}


QMG::Point 
QMG::PatchMath::normal(const Point& paramcoord) const {
  if (gdim_ != di_ - 1)
    throw_error("Normal called for non-full-dimensional patch");
  Point direc;
  Point returnval;
  if (gdim_ == 1) {
    direc[0] = 1;
    Point tan = direc_deriv(paramcoord, direc);
    returnval[0] = tan[1];
    returnval[1] = -tan[0];
  }
  else {
    direc[0] = 1;
    direc[1] = 0;
    Point tan1 = direc_deriv(paramcoord, direc);
    direc[0] = 0;
    direc[1] = 1;
    Point tan2 = direc_deriv(paramcoord, direc);
    returnval = Point::cross_product(tan1, tan2);
  }
  if (returnval.l2norm() == 0)
    throw_error("zero normal");
  returnval.normalize();
  if (!is_forward_) {
    for (int i = 0; i < di_; ++i) 
      returnval[i] = -returnval[i];
  }
  return returnval;
}

QMG::Point QMG::PatchMath::normal_at_barycenter() const {
  return normal(barycenter_param());
}

QMG::Point QMG::PatchMath::barycenter_param() const {
  Point returnval;
  if (gdim_ == 1)
    returnval[0] = 0.5;
  else if (gdim_ == 2) {
    if (ptype_ == BEZIER_TRIANGLE) {
      returnval[0] = 0.33333333333333333;
      returnval[1] = 0.33333333333333333;
    }
    else {
      returnval[0] = 0.5;
      returnval[1] = 0.5;
    }
  }
  return returnval;
}

 
int QMG::PatchMath::num_control_points() const { 
  if (gdim_ == 0) {
    return 1;
  }
  else if (gdim_ == 1) {
    return degree1_ + 1;
  }
  else if (ptype_ == BEZIER_TRIANGLE) {
    return (degree1_ + 1) * (degree1_ + 2) / 2;
  }
  else {
    return (degree1_ + 1) * (degree2_ + 1);
  }
}