align.cpp

算断裂的

  // that is not on the box face.  We check whether it is OK by testing
  // whether the volume of the simplices induced with previously protected
  // boxes diminishes too much.

  bool warp_ok(Real dist,
    const Point& warp_point,
    const ActiveBox& b,
    Base3 warp_rindex,
    const BoxFaceDag& dag,
    BoxFaceDag::Index linkparent,
    const SimpComplex_Under_Construction& sc,
    Real asp_degrade,
    Logstream& logstr,
    double scfac,
    double tol,
    // workspace
    vector<SimpComplex::VertexOrdinalIndex>& volist) {
    if (logstr.verbosity() >= 4) {
      logstr.str() << "In warp_ok for box " << b << " linkparent is bfdi#" 
        << b.link_parent() << "# ";
    }
    if (dist == 0.0) {
      if (logstr.verbosity() >= 4) {
        logstr.str() << " dist=0 return true\n";
      }
      return true;
    }

    int di = sc.embedded_dim();
    int bdim = b.dim();
    int bdimpower3 = Base3::power3(bdim);
    Real realwidth = b.real_width();
    Point unwarped_point = compute_l1_closepoint(warp_point, b, warp_rindex);

#ifdef DEBUGGING
    Real computed_dist = Point::l1dist(warp_point, unwarped_point); 
    if (fabs(computed_dist - dist) > 1e3 * scfac * tol) {
      throw_error("Problem computing unwarped point");
    }
#endif

    Point currentpoint;
    Real matrix1_space[2 * MAXDIM * MAXDIM];
    Matrix matrix1(di, di, matrix1_space);
    Matrix matrix2(di, di, &matrix1_space[MAXDIM * MAXDIM]);



    for (int ri = 0; ri < bdimpower3; ++ri) {
      Base3 face_rindex = ri;
      if (bdim > 0) {
        if (face_rindex.dim() != bdim - 1 ||
          Base3::merge(face_rindex, warp_rindex) == face_rindex)
          continue;
      }
      currentpoint = b.face_real_lowerleft(face_rindex);

      // Row 0 gets overwritten below if bdim = 0.

      for (int jj = 0; jj < di; ++jj) {
        matrix1(0,jj) = currentpoint[jj] - warp_point[jj];
        matrix2(0,jj) = currentpoint[jj] - unwarped_point[jj];
      }
      int matrixrow = 1;
      {
        int actual_dim = -1;
        for (int ii = 0; ii < bdim; ++ii) {
          while (b.flatdim(++actual_dim));
          if (face_rindex[ii] < 2) continue;
          currentpoint[actual_dim] += realwidth;
          for (int jj = 0; jj < di; ++jj) {
            matrix1(matrixrow, jj) = currentpoint[jj] - warp_point[jj];
            matrix2(matrixrow, jj) = currentpoint[jj] - unwarped_point[jj];
          }
          ++matrixrow;
        }
      }
#ifdef DEBUGGING
      if (bdim > 0 && matrixrow != bdim)
        throw_error("Inconsistent matrixrow in align");
#endif

      if (linkparent < 0) {
#ifdef DEBUGGING
        if (bdim != di)
          throw_error("Inconsistent dimensions in align");
#endif
        Real dt1 = matrix1.determinant();
        Real dt2 = matrix2.determinant();
        if (logstr.verbosity() >= 4) {
          logstr.str() << " rind = ";
          warp_rindex.output(logstr.str(), bdim);
          logstr.str() << " facerindex = ";
          face_rindex.output(logstr.str(), bdim);
          if (logstr.verbosity() >= 7) {
            logstr.str() << '\n' << matrix1 << '\n' << matrix2;
          }
          logstr.str() << " dets = " << dt1 << ',' << dt2 << " asp_degrade = " << asp_degrade;
        }

        if (dt1 * dt2 <= 0.0 || fabs(dt1) < asp_degrade * fabs(dt2)) {
          if (logstr.verbosity() >= 4)
            logstr.str() << " return false 1 \n";
          return false;
        }
      }
      else {
        
        int scrap1, scrap2;
        dag.trace_path_to_root(linkparent, volist, scrap1, scrap2);
        
#ifdef DEBUGGING
        if (volist.size() != di - bdim) 
          throw_error("Inconsistent dimensions in align 2");
#endif

        for (int ii = 0; ii < di - bdim; ++ii) {
          SimpComplex::VertexOrdinalIndex thisvert = volist[ii];
          for (int jj = 0; jj < di; ++jj) {
            Real parentcoord = sc.real_coord_o(thisvert, jj);
            matrix1(ii + bdim, jj) = parentcoord - warp_point[jj];
            matrix2(ii + bdim, jj) = parentcoord - unwarped_point[jj];
          }
        }
        Real dt1 = matrix1.determinant();
        Real dt2 = matrix2.determinant();
        if (logstr.verbosity() >= 4) {
          logstr.str() << " rind = ";
          warp_rindex.output(logstr.str(), bdim);
          logstr.str() << " face_rindex = ";
          warp_rindex.output(logstr.str(), bdim);
          logstr.str() << " anc vtx = ";
          for (int l1 = 0; l1 < di - bdim; ++l1)
            logstr.str() << "ord#" << volist[l1] << "# ";
          if (logstr.verbosity() >= 7) {
            logstr.str() << '\n' << matrix1 << '\n' << matrix2;
          }
          logstr.str() << " dets = " << dt1 << ',' << dt2 << " asp_degrade = " << asp_degrade;
        }            
        if (dt1 * dt2 <= 0.0 || fabs(dt1) < asp_degrade * fabs(dt2)) {
          if (logstr.verbosity() >= 4) 
            logstr.str() << " return false\n";
          return false;
        }
      }
    }
    if (logstr.verbosity() >= 4) {
      logstr.str() << " return true\n";
    }
    return true;
  }

  
  // ------------------------------------------------------------------
  // launch_get_subbox_weight
  // When a facet of a box is launched as a new box, this routine computes 
  // the weight of its faces based on the weight of the parent.

  Weight_t launch_get_subbox_weight(const ActiveBox& ab,
    int parentbdim,
    Base3 closeface,
    Base3 rindex,
    int reldim,
    int facetval) {
    
    int newind = 0;
    Base3 parent_rindex = 0;      
    int rindex_dim = 0;
    int overlap_count = 0;
    for (int parentind = 0; parentind < parentbdim; ++parentind) {
      int thisdigit = (parentind == reldim)? facetval : rindex[newind++];
      parent_rindex.set_digit(parentind, thisdigit);
      if (thisdigit < 2) {
        if (thisdigit == closeface[parentind])
          ++overlap_count;
      }
      else {
        ++rindex_dim;
      }
    }

    Weight_t newwt = ab.weight(parent_rindex) / (parentbdim - rindex_dim - overlap_count);
    return newwt;
  }

  
  // ------------------------------------------------------------------
  // class ActiveBox_For_Launch
  // This is an active box whose facets are about to be launched as
  // new active boxes.  This class contains extra information
  // about computing weights, etc.
  // Specialized postprocessing inserts new box close faces into priority queue.
    

    
  class ActiveBox_For_Launch : public ActiveBox_As_Parent {
  private:
    int bdim_;
    Base3 closeface_;
    int reldim_, facetval_;
    mutable Brep::Face_Spec chambers[MAXDIM][2];
    mutable bool is_front_[MAXDIM][2];
    Brep::Face_Spec alignfspec_;
    IncidencePriorityQueue_ABC& ipq_;
  public:
    ActiveBox_For_Launch(const ActiveBox& ab, 
      Base3 closeface,
      const Brep::Face_Spec& alignfspec,
      IncidencePriorityQueue_ABC& ipq) :
    ActiveBox_As_Parent(ab), 
    bdim_(ab.dim()),
    closeface_(closeface),
    alignfspec_(alignfspec),
    ipq_(ipq) {
      for (int i = 0; i < bdim_; ++i) {
        for (int j = 0; j < 2; ++j)
          chambers[i][j] = Brep::Face_Spec(-2,-2);
      }
    }
    void set_reldim_facetval(int r, int f) {
      reldim_ = r; facetval_ = f;
    }
    Weight_t get_subbox_weight(Base3 ri) const {
      return launch_get_subbox_weight(*this, bdim_, closeface_, ri, reldim_, facetval_);
    }
    pair<Brep::Face_Spec,bool> chamber_containing_subbox(const ActiveBox& testbox,
      ActiveBoxVec::Create_Subbox_Workspace& wkspa) const;
    virtual void specialized_postprocess_subbox(const ActiveBox& newb,
      int lowest_patch_dim) const {
      if (lowest_patch_dim <= alignfspec_.fdim())
        ipq_.insert_box_closefaces(newb);
    }
  };
  
  
  pair<Brep::Face_Spec, bool> 
    ActiveBox_For_Launch::
    chamber_containing_subbox(const ActiveBox& testbox,
    ActiveBoxVec::Create_Subbox_Workspace& wkspa)  const {
    bool found = false;
    for (int i = 0; i < bdim_; ++i) {
      if (i == reldim_) continue;
      
      for (int j = 0; j < 2; ++j) {
        if (chambers[i][j].faceind() > -2) {
          found = true;
          chambers[reldim_][facetval_] = chambers[i][j];
          is_front_[reldim_][facetval_] = is_front_[i][j];
          break;
        }
      }
      if (found) break;
    }
    if (!found) {
      pair<Brep::Face_Spec,bool> rval = 
        ActiveBox_As_Parent::chamber_containing_subbox(testbox, wkspa);
      chambers[reldim_][facetval_] = rval.first;
      is_front_[reldim_][facetval_] = rval.second;
    }
    return pair<Brep::Face_Spec,bool>(chambers[reldim_][facetval_], 
      is_front_[reldim_][facetval_]);
  }

  
  // ------------------------------------------------------------------
  // compute_facet_address
  // Computes a box-address for a facet of an active box; used to
  // launch the facets as new boxes.


  void compute_facet_address(BoxAddress& bxa, 
    const ActiveBox& b,
    Base3 subfaceindex) {
    bxa = b.boxaddress();
    int di = b.embedded_dim();
    int reldim = -1;
    for (int i = 0; i < di; ++i) {
      if (!b.flatdim(i)) {
        ++reldim;
        if (subfaceindex[reldim] < 2) {
          bxa.flatdim().set_bit(i);
          if (subfaceindex[reldim] == 1) {
            bxa.lowerleft().increment(i, b.iwidth());
          }
        }
      }
    }
  }



  // ------------------------------------------------------------------
  // protect_box_and_launch_subfaces
  // This routine protects a box and launches its facets as new boxes.



  void protect_box_and_launch_subfaces(const ActiveBox& b,
    SimpComplex::VertexOrdinalIndex vertnum_o,
    Base3 closeface_rindex,
    const Brep::Face_Spec& apex,
    ActiveBoxVec& orbit,
    ActiveBoxVec& idle_vec,
    BoxFaceDag& dag,
    IncidencePriorityQueue& ipq,
    ActiveBoxVec::Create_Subbox_Workspace& wkspa) {

    BoxFaceDag::Index newind = dag.add_protected_box(b, vertnum_o);
    if (wkspa.logstr.verbosity() >= 2) { ////
      wkspa.logstr.str() << "New protected box " << b << " face ";
      closeface_rindex.output(wkspa.logstr.str(), b.dim());
      wkspa.logstr.str() << " at bfdi#" << newind << "#  parent is bfdi#"
        << b.link_parent() << "#\n";
    }
    int bwidth = b.iwidth();
    int bdim = b.dim();
    int actualdim = -1;
    Base3 fullindex = Base3::power3(bdim) - 1;
    BoxAddress bxa;
    

    ActiveBox_For_Launch abl(b, closeface_rindex, apex, ipq);
    
    for (int reldim = 0; reldim < bdim; ++reldim) {
      while (b.flatdim(++actualdim));
      for (int facetval = 0; facetval < 2; ++facetval) {
        Base3 subface_rindex = fullindex.modify_digit(reldim, facetval);
        if (Base3::merge(subface_rindex, closeface_rindex) != subface_rindex) {
          compute_facet_address(bxa, b, subface_rindex);
          abl.set_reldim_facetval(reldim, facetval);
          ActiveBoxVec::create_subbox(abl, bxa, bwidth, apex.fdim(),
            newind, facetval, actualdim, orbit, idle_vec, wkspa);
        }
      }
    }
  }

  
  // ------------------------------------------------------------------
  // protect_box_and_launch_subfaces_lastphase
  // Same as the above, except we are in the last phase so
  // the create_subbox routine is simplified.


  void protect_box_and_launch_subfaces_lastphase(const ActiveBox& b,
    SimpComplex::VertexOrdinalIndex vertnum_o,
    Base3 closepoint_rindex,
    const Brep::Face_Spec& apex,
    ActiveBoxVec& orbit,
    BoxFaceDag& dag,
    IncidencePriorityQueue_LastPhase& ipq,
    Logstream& logstr) {
    
    BoxFaceDag::Index newind = dag.add_protected_box(b,vertnum_o);
    if (logstr.verbosity() >= 2) { ////
      logstr.str() << "New protected interior box ";
      if (logstr.verbosity() >= 4)
        b.full_dump(logstr.str());
      else
        logstr.str() << b;
      logstr.str() << " closeface = ";
      closepoint_rindex.output(logstr.str(),b.dim()); 
      logstr.str() << " at bfdi#" << newind << "#  parent is bfdi#"
        << b.link_parent() << "#\n";
    }
    int bwidth = b.iwidth();
    int bdim = b.dim();
    int actualdim = -1;
    Base3 fullindex = Base3::power3(bdim) - 1;
    BoxAddress bxa;

    
    ActiveBox_For_Launch abl(b, closepoint_rindex, apex, ipq);
    
    for (int reldim = 0; reldim < bdim; ++reldim) {
      while (b.flatdim(++actualdim));
      for (int facetval = 0; facetval < 2; ++facetval) {
        Base3 subface_rindex = fullindex.modify_digit(reldim, facetval);
        if (Base3::merge(subface_rindex, closepoint_rindex) != subface_rindex) {
          compute_facet_address(bxa, b, subface_rindex);          
          abl.set_reldim_facetval(reldim, facetval);
          ActiveBoxVec::create_subbox_completely_interior(abl, bxa, 
            bwidth, orbit,  newind, facetval, actualdim, logstr);
        }
      }
    }
  }
  
  // ------------------------------------------------------------------
  // struct DeduplicationMapKey
  // Used as a key for a map to deduplicate a box stack.  Boxes are
  // duplicates if they have the same address and the same highest
  // incidence number.

    struct DeduplicationMapKey {
    BoxAddress bxa;
    PatchTable::Index highest_inc;
    bool operator<(const DeduplicationMapKey& other) const;
  };

  bool DeduplicationMapKey::operator<(const DeduplicationMapKey& other) const {
    if (highest_inc < other.highest_inc)
      return true;
    if (highest_inc > other.highest_inc)
      return false;
    if ((UInt32) (bxa.flatdim()) < (UInt32) (other.bxa.flatdim()))
      return true;
    if ((UInt32) (bxa.flatdim()) > (UInt32) (other.bxa.flatdim()))
      return false;
    int sgndiff = IntCoord::compare(bxa.lowerleft(), other.bxa.lowerleft());
    if (sgndiff < 0)
      return true;
    return false;
  }