qboxstack.cpp

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

// ------------------------------------------------------------------
// qboxstack.cpp
//
// This file contains member functions qboxstack.h.  These are classes
// for boxes in the mesh generator, and vectors of boxes.
// ------------------------------------------------------------------
// 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 "qboxstack.h"


namespace QMG {
  namespace MG {
    using namespace QMG;

    
    //global function in this file


    void boxstack_init_static_data(int di); 
  }
}

// File scope functions.

namespace {

  using namespace QMG;
  using namespace QMG::MG;

#ifdef DEBUGGING
  ostream* debug_ostream;  
  bool dump_everything;
  int create_subbox_invocation_count;
#endif



  // ------------------------------------------------------------------
  // Routine to compute a random point in the unit octahedron.
  // Selection is not uniform -- I don't know an easy way to do this uniformly.
  Point select_rand_point_in_unit_octahedron(int di) {
    Point p;
    {
      for (int i = 0; i < di; ++i)
        p[i] = random_real() * 2.0 - 1.0;
    }
    Real t = random_real();
    {
      for (int i = 0; i < di; ++i)
        t += fabs(p[i]);
    }
    {
      for (int i = 0; i < di; ++i)
        p[i] /= t;
    }
    return p;
  }

}


// ------------------------------------------------------------------
// contains_point_in_ex
// This routine determines whether a point is contained in the expanded
// version of a box.

bool 
QMG::MG::ActiveBox::contains_point_in_ex(const Point& pt, 
                                         Real param_uncertainty,
                                         Real scfac) const { 

  Real l1dist = 0.0;
  Real mult = BoxAddress::multiplier();
  Real rwidth = mult / (1<<iwidth());
  for (int actualdim = 0; actualdim < di_; ++actualdim) {
    Real lower = ((Real) this -> boxaddress().lowerleft(actualdim)) * mult / 
      (1 << MAXLEV) + BoxAddress::global_coordinate_origin(actualdim);
    Real upper = lower + ( (flatdim(actualdim)) ? 0 : rwidth);
    if (pt[actualdim] < lower)
      l1dist += lower - pt[actualdim];
    else if (pt[actualdim] > upper)
      l1dist += pt[actualdim] - upper;
  }
  return (l1dist <= rwidth * expa() - scfac * param_uncertainty);
}


// ------------------------------------------------------------------
// This routine selects a random point in the expanded neighborhood
// of a box.

QMG::MG::Point 
QMG::MG::ActiveBox::select_rand_point() const {

  Real mult = BoxAddress::multiplier();
  Real rwidth = mult / (1<<iwidth());
  Point pt =  select_rand_point_in_unit_octahedron(di_);

#ifdef DEBUGGING
  if (dump_everything) {
    *debug_ostream << " in select rand initalpt = " << pt << "expa = " << expa() 
      << " mult = " << mult << " rwidth = " << rwidth << '\n';
  }
#endif

  for (int actualdim = 0; actualdim < di_; ++actualdim) {
    Real lower = ((Real) this -> boxaddress().lowerleft(actualdim)) * mult / 
      (1 << MAXLEV) + BoxAddress::global_coordinate_origin(actualdim);
    Real thisw =  (flatdim(actualdim)) ? 0.0: rwidth;
    pt[actualdim] = pt[actualdim] * rwidth * expa() + 
      random_real() * thisw + lower;
  }
  return pt;
}

// ------------------------------------------------------------------
// This routine computes the real coordinates of the lower left corner
// of the box.

QMG::MG::Point
QMG::MG::BoxAddress::real_lowerleft() const {
  Point p;
  for (int ii = 0; ii < di_; ++ii) {
    p[ii] = lowerleft_[ii] * multiplier_ / (1 << MAXLEV) +
      global_coordinate_origin_[ii];
  }
  return p;
}



// ------------------------------------------------------------------
// This routine computes the real coordinates of the lower left
// corner of a face of a box.

QMG::MG::Point 
QMG::MG::ActiveBox::face_real_lowerleft(Base3 rindex) const {
  int bdim = this -> dim();
  int actual_dim = -1;
  BoxAddress bxa = this -> boxaddress();
  int iwidth = this -> iwidth();

  for (int ii = 0; ii < bdim; ++ii) {
    while (this -> boxaddress().flatdim(++actual_dim));
    if (rindex[ii] == 1)
      bxa.lowerleft().increment(actual_dim,iwidth);
  }
  return bxa.real_lowerleft();
}



// ------------------------------------------------------------------
// This routine computes the weight of a box face.  Use the static
// weights for fulldimensional boxes, properly adjusted if the box
// lies at the boundary of the top box.  Else use the weights stored
// in the box vector.

QMG::MG::Weight_t
QMG::MG::ActiveBox::weight(Base3 rindex) const {
#ifdef RANGECHECK
  if ((UInt32) rindex >= Base3::power3(this -> dim())) {
    throw_error("Arg out of range in weight()");
  }
#endif
  if (this -> dim() == di_) {
    Weight_t w = fulldim_weights_[rindex];
    UInt32 rightmost = (1 << MAXLEV) - (1 << (MAXLEV - this -> iwidth()));
    for (int j = 0; j < di_; ++j) {
      int ri = rindex[j];
      if (ri == 0 && this -> boxaddress().lowerleft(j) == 0)
        w *= 2;
      else if (ri == 1 && this -> boxaddress().lowerleft(j) == rightmost)
        w *= 2;
    }
    return w;

  }
  else {
    return (this -> owner_).weights_[(this -> owner_).boxvec_[this -> idx_].weights_begin + 
      (UInt32) rindex];
  }
}


// ------------------------------------------------------------------
// Initialize static data in this file.


void QMG::MG::ActiveBox::init_static_data(int di) {
  di_ = di;
  active_box_count_ = 0;
  for (Base3 rindex = 0; rindex < Base3::power3(di); ++rindex) {
    int d = rindex.dim();
    fulldim_weights_[rindex] = 1.0 / ((Weight_t) (1 << (di - d)));
  }
}



// ------------------------------------------------------------------
// Returns the weight of a box face as a nonconstant reference (so
// that it can be updated).  Valid only for non-full-dimensional boxes.

QMG::MG::Weight_t& 
QMG::MG::NonConstActiveBox::weight_nonconst(Base3 rindex) {
#ifdef RANGECHECK
  if ((UInt32) rindex >= Base3::power3(this -> dim())) {
    throw_error("Arg out of range in weight_nonconst()");
  }
#endif
#ifdef DEBUGGING
  if (this -> dim() == di_) {
    throw_error("Attempt to set weight for fulldimensional box");
  }
#endif
  return (const_cast<ActiveBoxVec&>(this -> owner_)).
    weights_[(this -> owner_).boxvec_[this -> idx_].weights_begin + 
    (UInt32) rindex];
}

 

// ------------------------------------------------------------------
// Loop_over_vertices: class for looping over the vertices of a box.

void 
QMG::MG::ActiveBox::Loop_over_vertices::operator++() {
  this -> ind_++;
  int reldim = 0;
  for (int j = 0; j < this -> di_; ++j) {
    int bit = (this -> bflatdim_[j])? 0 : ((ind_ >> (reldim++)) & 1);
    this -> intcoord_[j] = this -> baseintcoord_[j];
    this -> realcoord_[j] = this -> baserealcoord_[j];
    if (bit) {
      this -> intcoord_.increment(j, this -> bwidth_);
      this -> realcoord_[j] += brealwidth_;
    }
  }
}

Base3 
QMG::MG::ActiveBox::Loop_over_vertices::rindex() const {
  UInt32 rind = 0;
  for (int j = 0; j < this -> bdim_; ++j) {
    if ((this -> ind_ >> j) & 1)
      rind += Base3::power3(j);
  }
  return (Base3) rind;
}


// ------------------------------------------------------------------
// Loop_over_incidences
// Class for looping over incidences of a box.



QMG::MG::ActiveBox::Loop_over_incidences::Loop_over_incidences(const ActiveBox& ab,
                                                               const PatchTable& patchtable,
                                                               const Brep::Face_Spec& fspec) :

ab_(ab),
patchtable_(patchtable),
fspec_(fspec) {
  for (this -> relnum_ = 0; 
  this -> relnum_ < (this -> ab_).num_patch(); 
  ++(this -> relnum_)) {
    PatchTable::Index p = (this -> ab_).patch(this -> relnum_).patchind();
    if ((this ->patchtable_).owner(p) == (this -> fspec_) && 
      (this -> ab_).patch(this -> relnum_).num_incidence() > 0)
      break;
  }
  this -> ii_ = 0;
}


void 
QMG::MG::ActiveBox::Loop_over_incidences::operator++() {
  this -> ii_++;
  if (this -> ii_ < (this -> ab_).patch(this -> relnum_).num_incidence())
    return;
  for (++(this ->relnum_); 
  (this -> relnum_) < (this -> ab_).num_patch(); 
  ++(this -> relnum_)) {
    PatchTable::Index p = (this -> ab_).patch((this -> relnum_)).patchind();
    if ((this -> patchtable_).owner(p) == (this -> fspec_) &&
      (this -> ab_).patch(this -> relnum_).num_incidence() > 0)
      break;
  }
  this -> ii_ = 0;
}

QMG::MG::ActiveBox 
QMG::MG::ActiveBoxVec::get_last_box_() const {
  return ActiveBox(*this, this -> boxvec_.size() - 1);
}



 

void QMG::MG::boxstack_init_static_data(int di) {
  ActiveBoxVec::init_static_data(di);
  ActiveBox::init_static_data(di);
  BoxAddress::init_static_data(di);
}




#ifdef NSF_DEMO
namespace QMG {
  namespace Nsf_Demo {
    using namespace QMG;
    using namespace QMG::MG;
    extern void show_box(const ActiveBox& b, int code);
  }
}
#endif


// ------------------------------------------------------------------
// push_top_box
// This routine computes the initial (top) box of the quadtree and
// pushes it into an ActiveBoxVec.  This is how mesh generation
// is initiated.

void QMG::MG::ActiveBoxVec::push_top_box(const Brep& g,
                                         IncidenceTable& inc_table,