qinctab.cpp

来自「算断裂的」· C++ 代码 · 共 559 行 · 第 1/2 页

// ------------------------------------------------------------------
// qinctab.cpp
//
// This file contains member functions for IncidenceTable.
// ------------------------------------------------------------------
// Author: Stephen A. Vavasis
// 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.
// ------------------------------------------------------------------

#ifdef _MSC_VER
#if _MSC_VER == 1200
#include "qmatvec.h"
#endif
#endif

#include "qinctab.h"
#include "qmatvec.h"
#include "qmg_gui.h"
#include <cstddef>


#ifdef DEBUGGING
namespace {
  int get_incidence_invocation_count;
}
#endif


// ------------------------------------------------------------------
// PatchMath_for_curve_on_triangle
// This routine allows us to treat a curve that sits on a patch as
// a separate Bezier curve.  The patch is described by two parameters
// alpha and beta which are the slopes of the curve.

class QMG::MG::IncidenceTable::PatchMath_for_curve_on_triangle : public QMG::PatchMath {
private:
  const PatchTable& patchtab_;
  PatchTable::Index patchind_;
  Real alpha_,beta_;
  int deg_;
  Real control_point_coord_(int cpnum, int d) const;
public:
  PatchMath_for_curve_on_triangle(const PatchTable& patchtab,
    PatchTable::Index patchind,
    Real alpha,
    Real beta) :
  PatchMath(BEZIER_TRIANGLE, IncidenceTable::get_patch_degree1_(patchtab, patchind), 
    0, 1, 3, true),
    patchtab_(patchtab),
    patchind_(patchind),
    alpha_(alpha),
    beta_(beta),
    deg_(IncidenceTable::get_patch_degree1_(patchtab,patchind))
  { }
  ~PatchMath_for_curve_on_triangle() { }
};


// ------------------------------------------------------------------
// The control point coordinates of the curve are interpolated
// from the patch's control points.
    
    
QMG::Real 
QMG::MG::IncidenceTable::PatchMath_for_curve_on_triangle::
control_point_coord_(int cpnum, int d) const {
  // Should use decasteljau here, but no obvious way to allocate workspace
  Real coef = 1.0;
  for (int jj = 0; jj < cpnum; ++jj)
    coef *= beta_;
  Real t = 0.0;
  int liftcp = (deg_ - cpnum) * (deg_ - cpnum + 1) / 2;
  int incliftcp = deg_ - cpnum + 2;
  for (int j = 0; j <= cpnum; ++j) {
    t += IncidenceTable::get_patch_control_point_coord_(patchtab_, patchind_, liftcp, d) 
      * coef;
    coef *= alpha_ / beta_ * (cpnum - j) / (j + 1);
    liftcp += incliftcp;
    ++incliftcp;
  }
  return t;
}




// ------------------------------------------------------------------
// PatchMath_for_curve_on_quad
// This routine allows us to treat a curve that sits on a quad patch as
// a separate Bezier curve.  The curve has a fixed v and varying u.
// The fixed value of v is parameter beta.

class QMG::MG::IncidenceTable::PatchMath_for_curve_on_quad : public QMG::PatchMath {
private:
  const PatchTable& patchtab_;
  PatchTable::Index patchind_;
  Real beta_;
  int deg1_, deg2_;
  Real control_point_coord_(int cpnum, int d) const;
public:
  PatchMath_for_curve_on_quad(const PatchTable& patchtab,
    PatchTable::Index patchind,
    Real beta) :
  PatchMath(BEZIER_TRIANGLE, IncidenceTable::get_patch_degree1_(patchtab, patchind), 
    0, 1, 3, true),
    patchtab_(patchtab),
    patchind_(patchind),
    beta_(beta),
    deg1_(IncidenceTable::get_patch_degree1_(patchtab,patchind)),    
    deg2_(IncidenceTable::get_patch_degree2_(patchtab,patchind))
  { }
  ~PatchMath_for_curve_on_quad() { }
};

    
// The control points of the curve are interpolated from
// control points of the patch.
    
QMG::Real 
QMG::MG::IncidenceTable::PatchMath_for_curve_on_quad::
control_point_coord_(int cpnum, int d) const {
  // Should use decasteljau here, but no obvious way to allocate workspace
  Real coef = 1.0;
  for (int jj = 0; jj < deg2_; ++jj)
    coef *= (1 - beta_);
  Real t = 0.0;
  int liftcp = cpnum;
  for (int j = 0; j <= deg2_; ++j) {
    t += IncidenceTable::get_patch_control_point_coord_(patchtab_, patchind_, liftcp, d)
      * coef;
    coef *= beta_ / (1 - beta_) * (deg2_ - j) / (j + 1);
    liftcp += deg1_ + 1;
  }
  return t;
}


// Functions to expose otherwise hidden patchtable fields.

int
QMG::MG::IncidenceTable::get_patch_degree1_(const PatchTable& patchtable,
                                           PatchTable::Index patchind) {
  return patchtable.patches_[patchind].degree1;
}


int
QMG::MG::IncidenceTable::get_patch_degree2_(const PatchTable& patchtable,
                                           PatchTable::Index patchind) {
  return patchtable.patches_[patchind].degree2;
}



QMG::PatchType
QMG::MG::IncidenceTable::get_patch_type_(const PatchTable& patchtable,
                                           PatchTable::Index patchind) {
  return patchtable.patches_[patchind].ptype;
}



QMG::Real 
QMG::MG::IncidenceTable::get_patch_control_point_coord_(const PatchTable& patchtable,
                                                    PatchTable::Index patchind,
                                                    int cpnum,
                                                    int d)  {
  int cpstart = patchtable.patches_[patchind].control_point_start;
  Brep::ControlPointIndex cpi = patchtable.control_point_list_[cpstart + cpnum];
  return patchtable.brep_.control_point(cpi, d);
}



bool QMG::MG::IncidenceTable::Key::operator<(const Key& other) const {
  if (patchind < other.patchind)
    return true;
  if (patchind > other.patchind)
    return false;
  if ((UInt32) flatdim < (UInt32) other.flatdim)
    return true;
  if ((UInt32) flatdim > (UInt32) other.flatdim)
    return false;
  int sgndif = IntCoord::compare(intcoord, other.intcoord);
  if (sgndif < 0)
    return true;
  if (sgndif > 0)
    return false;
  if (seqno < other.seqno)
    return true;
  return false;
}


int QMG::MG::IncidenceTable::compare_index(const Index& i1, const Index& i2) {
  return (i1 -> first.seqno) - (i2 -> first.seqno);
}


QMG::MG::IncidenceTable::IncidenceTable() : sequence_number_(0) {
#ifdef DEBUGGING
  get_incidence_invocation_count = 0;
#endif
}

// This the main routine for IncidenceTable: it computes
// the incidence (intersection) between a box entity
// and a brep entity.  The box entity is described by
// flatdim and midpoint coord, and midpointcoord_real
// is midpoint_coord converted to real coordinates.
// The box entity is patchind.
// The result is placed in output_vec.


void     
QMG::MG::IncidenceTable::get_incidence(Mask flatdim, 
                                       const IntCoord& midpointcoord, 
                                       const Point& midpointcoord_real,
                                       const PatchTable& patchtable,
                                       PatchTable::Index patchind,
                                       vector<Index>& output_vec, // output variable
                                       Logstream& logstr, 
                                       double scfac,
                                       double tol,
                                       PatchMath::Workspace& workspace,
                                       vector<PatchMath::IntersectRecord>& os_workspace,
                                       Meshgen_gui& gui) {
#ifdef DEBUGGING
  ++get_incidence_invocation_count;
  int invoc_count = get_incidence_invocation_count;
#endif


  int di = patchtable.embedded_dim();
  int patchdim = patchtable.patches_[patchind].gdim;
  int affinedim = flatdim.dim();
  const int GUI_UPDATE_FREQUENCY = 20;

  // If the dimension is too low, no intersection.

  if (affinedim + patchdim < di)
    return;

  // Look up the incidence; see if already computed.

  Key key;
  key.intcoord = midpointcoord;
  key.flatdim = flatdim;
  key.patchind = patchind;
  key.seqno = 0;
  for (int j = 0; j < di; ++j) {
    if (!key.flatdim[j]) {
      key.intcoord[j] = 0;
    }
  }
  typedef map<Key,Data> TTTYPE;
  typedef TTTYPE::iterator TTTYPE_IT;

  Index rangestart = table_.lower_bound(key);
  key.seqno = sequence_number_;
  Index rangeend = table_.upper_bound(key);
  if (rangestart != rangeend) {
    // Check if the first record indicates "no intersection."
    if (rangestart -> second.param_coord[0] == code_for_no_intersection()) 
      return;

    // Otherwise return all the incidence records.

    for (Index i1 = rangestart; i1 != rangeend; ++i1) {
      output_vec.push_back(i1);
    }
    return;
  }