split1.cpp

来自「算断裂的」· C++ 代码 · 共 324 行

// ------------------------------------------------------------------
// split1.cpp
//
// This file contains the routine for splitting a box.
// ------------------------------------------------------------------
// 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 {

  using namespace QMG;
  using namespace QMG::MG;


  // overrides virtual function: compute subbox weights
  // based on how the subbox overlaps with the parent box.

  Weight_t split_get_subbox_weight(const ActiveBox& b,
    Base3 rindex,
    int bdim,
    int newboxind) {
    Weight_t factor = 1.0;
    Base3 source_rindex = rindex;
    for (int j = 0; j < bdim; ++j) {
      int bit = (newboxind >> j) & 1;
      if ((bit && rindex[j] == 0) || (!bit && rindex[j] == 1)) {
        factor *= 0.5;
        source_rindex = source_rindex.modify_digit(j, 2);
      }
    }
    return b.weight(source_rindex) * factor;
  }

  class ActiveBox_For_Split : public ActiveBox_As_Parent {
  private:
    int bdim_;
    int newboxind_;
    const ActiveBox& auxbox_;
    mutable int numinc0_;
    mutable int numinc1_;
    mutable Brep::Face_Spec interior_determined_;
    mutable bool is_front_;

  public:
    ActiveBox_For_Split(const ActiveBox& ab,
      const ActiveBox& auxbox) :
        ActiveBox_As_Parent(ab),
        bdim_(ab.dim()),
        auxbox_(auxbox),
        interior_determined_(-2,-2) { }
    ~ActiveBox_For_Split() { }
    void set_boxind(int boxind) {newboxind_ = boxind;}
    pair<Brep::Face_Spec, bool>
      chamber_containing_subbox(const ActiveBox& testbox,
      ActiveBoxVec::Create_Subbox_Workspace& wkspa) const;

    Weight_t get_subbox_weight(Base3 ri) const {
      return split_get_subbox_weight(*this, ri, bdim_, newboxind_);
    }
    void incloop_init(int relnum) const {
      relnum_ = relnum;
      numinc0_ = this -> patch(relnum).num_incidence();
      numinc1_ = auxbox_.patch(relnum).num_incidence();
    }
    int incloop_numinc() const { return numinc0_ + numinc1_; }
    IncidenceTable::Index incloop_getinc(int ind) const {
      return (ind < numinc0_)?
        this -> patch(relnum_).incidence(ind) :
        auxbox_.patch(relnum_).incidence(ind - numinc0_);
    }
  };

  // This member function overrides the base class
  // chamber_containing subbox.  The idea is that
  // whatever chamber contains one of the split boxes
  // contains them all since they have a common
  // point.  Therefore, we need call the base class
  // function just once, and then cache the result.

  pair<Brep::Face_Spec, bool>
    ActiveBox_For_Split::chamber_containing_subbox(const ActiveBox& testbox,
    ActiveBoxVec::Create_Subbox_Workspace& wkspa) const {
    
    if (interior_determined_.faceind() == -2) {
      pair<Brep::Face_Spec, bool> rval = 
      ActiveBox_As_Parent::chamber_containing_subbox(testbox, wkspa);
      interior_determined_ = rval.first;
      is_front_ = rval.second;
    }
    return pair<Brep::Face_Spec, bool>(interior_determined_, is_front_);
  }
  

  class ActiveBox_For_Split_Completely_Interior : public ActiveBox_As_Parent {
  private:
    int bdim_;
    int newboxind_;
  public:
    explicit ActiveBox_For_Split_Completely_Interior(const ActiveBox& ab) :
        ActiveBox_As_Parent(ab),
        bdim_(ab.dim())
        { }
    ~ActiveBox_For_Split_Completely_Interior() { }
    void set_boxind(int boxind) {newboxind_ = boxind;}
    Weight_t get_subbox_weight(Base3 ri) const {
      return split_get_subbox_weight(*this, ri, bdim_, newboxind_);
    }
  };
}



#ifdef DEBUGGING
namespace QMG {
  namespace Local_to_PatchMath_CPP {
    extern ostream* debug_str;
    extern bool dump_everything;
  }
}
#endif


void 
QMG::MG::ActiveBoxVec::
split_box_and_push_children(int cutoffdim,
                            ActiveBoxVec& tempvec,
                            ActiveBoxVec& smallboxvec,
                            ActiveBoxVec& idlevec,
                            const ActiveBox& b,
                            IncidenceTable& inc_table,
                            Create_Subbox_Workspace& wkspa) {
#ifdef DEBUGGING
  using QMG::Local_to_PatchMath_CPP::debug_str;
  using QMG::Local_to_PatchMath_CPP::dump_everything;
#endif
  
  
  // Compute intersection of the splitting faces with the
  // new box faces.  Store these intersections in a temporary box.
  
  
  // Save the initial sizes so that we can restore the vector when done.
  
  int tempvec_patch_initsize = tempvec.patchvec_.size();
  int tempvec_incidence_initsize = tempvec.incidencevec_.size();
  int tempvec_box_initsize = tempvec.boxvec_.size();
  
  int di = b.embedded_dim();
  int diexp = 1 << di;
  int iwidth = b.iwidth();
  Real rwidth = b.real_width();
  
  if (wkspa.logstr.verbosity() >= 2) {
    wkspa.logstr.str() << "Splitting " << b << '\n';
  }
  
  
  
  // Create a new temporary box, initially with no patches.
  
  tempvec.boxvec_.push_back(Box_Data());
  Box_Data& tempboxd = tempvec.boxvec_.back();
  tempboxd.box_completely_interior = false;
  tempboxd.patchvec_begin = tempvec.patchvec_.size();
  tempboxd.num_patch = 0;
  
  Patch_Data newptcd;

  int itable1[18] = {1,2,4,  3,3,3,3,3, 5,5,5,5,5,     6,6,6,6,6};
  int itable2[18] = {2,6,18, 2,6,5,7,8, 2,11,18,19,20, 6,15,18,21,24};

  int itablesize = (di == 2)? 2 : 18;

  Mask bflatdim = b.flatdim();
  BoxAddress bxa;
  
  for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
    // Create a patch in the new box.
    
    Patch_in_ActiveBox ptc = b.patch(relnum);
    
    newptcd.inc_begin = tempvec.incidencevec_.size();
    newptcd.patchnum = ptc.patchind(); 
    // Loop over all flatdims/masks in the table,
    
    for (int jj = 0; jj < itablesize; ++jj) {
      Mask flatdim = (Mask) itable1[jj];
      Base3 incr = (Base3) itable2[jj];
      
      
      // Get the box midpoint coordinates.
      
      bxa = b.boxaddress();
      bool validi = true;
      for (int kk = 0; kk < di; ++kk) {
        if (incr[kk] == 2) {
          bxa.lowerleft().increment(kk, iwidth + 1);
        }
        else if (incr[kk] == 1) {
          bxa.lowerleft().increment(kk, iwidth);
        }
        if (bflatdim[kk] && incr[kk] > 0)
          validi = false;
      }
      if (validi) {
        Point midpointcoord_real = bxa.real_lowerleft();
        
        inc_table.get_incidence(flatdim, bxa.lowerleft(), midpointcoord_real, 
          wkspa.patchtable, newptcd.patchnum, tempvec.incidencevec_, wkspa.logstr, 
          wkspa.scfac, wkspa.tol, wkspa.rworkspace, wkspa.output_stack, wkspa.gui);
      }
    }

    newptcd.num_inc = tempvec.incidencevec_.size() - newptcd.inc_begin;
    ++tempboxd.num_patch;
    tempvec.patchvec_.push_back(newptcd);
  }
  
  
  // Create the new boxes.
  
  int bdim = b.dim();
  int bdimexp = 1<<bdim;
  
  ActiveBox auxbox(tempvec, tempvec_box_initsize);
  ActiveBox_For_Split bspl(b, auxbox);
  
  BoxFaceDag::Index parentind = b.link_parent();
  unsigned int parentside = b.link_parent_side();
  unsigned int parentdim = b.link_parent_dim();
  
  for (int binboxind = 0; binboxind < bdimexp; ++binboxind) {
    // compute the new box address.
    
    bxa = b.boxaddress();
    {
      int k = 0;
      for (int i = 0; i < di; ++i) {
        if (!b.flatdim(i) && (binboxind & (1 << (k++)))) {
          bxa.lowerleft().increment(i, iwidth+1);
        }
      }
    }
    
    
    bspl.set_boxind(binboxind);
    
    // Create the new box.
    
    create_subbox(bspl, bxa, iwidth + 1, cutoffdim, parentind, parentside, parentdim,
      smallboxvec, idlevec,wkspa);
 
  }
  
  // Free up the temporary space.
  
  tempvec.delete_last_box_();                                           
}


void 
QMG::MG::ActiveBoxVec::split_completely_interior_box_and_push_children(const ActiveBox& b, 
    ActiveBoxVec& smallboxvec,
    Logstream& logstr) {

  int di = b.embedded_dim();
  int iwidth = b.iwidth();
  int bdim = b.dim();
  int bdimexp = 1<<bdim;

  BoxAddress bxa;

  if (logstr.verbosity() >= 2)
    logstr.str() << "splitting completely interior box " << b << '\n';


  ActiveBox_For_Split_Completely_Interior bspl(b);


  for (int binboxind = 0; binboxind < bdimexp; ++binboxind) {
    // compute the new box address.

    bxa = b.boxaddress();
    {
      int k = 0;
      for (int i = 0; i < di; ++i) {
        if (!b.flatdim(i) && (binboxind & (1 << (k++)))) {
          bxa.lowerleft().increment(i, iwidth+1);
        }
      }
    }
  
     // Create the new box.

    bspl.set_boxind(binboxind);
    create_subbox_completely_interior(bspl, bxa, 
      iwidth + 1, smallboxvec, b.link_parent(), b.link_parent_side(),
      b.link_parent_dim(), logstr);
  }
}