crtsubbox.cpp

来自「算断裂的」· C++ 代码 · 共 1,123 行 · 第 1/3 页

// ------------------------------------------------------------------
// crtsubbox.cpp
//
// This file contains member functions from qboxstack.h for creating
// subboxes. 
// ------------------------------------------------------------------
// 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"
#include "qunionfind.h"

namespace QMG {
  namespace MG {
    using namespace QMG;

    
    //global function in this file


    void boxstack_init_static_data(int di); 
  }
}

// File scope functions.


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


namespace {

  using namespace QMG;
  using namespace QMG::MG;


  // ------------------------------------------------------------------
  // Combinatorial merge: merge each patch with its child patches.  This leaves
  // all bounding faces in the same connected component with the same label.
  // The return value is left in workspace variable wkspa.uf.

  void combinatorial_merge(const ActiveBox& b,
    ActiveBoxVec::Create_Subbox_Workspace& wkspa) {

 
    wkspa.map_to_relnum.reset(-1);
  
    for (int relnum = 0; relnum < b.num_patch(); relnum++) {
      PatchTable::Index p = b.patch(relnum).patchind();
      wkspa.map_to_relnum[p] = relnum;
      // Loop over the children patches.
      for (PatchTable::Loop_over_children loop(wkspa.patchtable, p);
      loop.notdone();
      ++loop) {
        PatchTable::Index c = loop.child_index();
        int relnum2 = wkspa.map_to_relnum[c];
        if (relnum2 >= 0) {
          wkspa.uf.merge(relnum, relnum2);
        }
      }
    }
  }  
  
  // Merge patches that belong to the same brep face.  Necessary else
  // alignment won't work correctly. Leave results in uf.

  void combinatorial_merge2(const ActiveBox& b,
    ActiveBoxVec::Create_Subbox_Workspace& wkspa) {
    {
      for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
        PatchTable::Index p = b.patch(relnum).patchind();
        wkspa.example_patchrelnum[wkspa.patchtable.owner(p)] = relnum;
      }
    }

    {
      for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
        PatchTable::Index p = b.patch(relnum).patchind();
        int relnum2 = wkspa.example_patchrelnum[wkspa.patchtable.owner(p)];
        wkspa.uf.merge(relnum, relnum2);
      }
    }
  }




    // ------------------------------------------------------------------
  // Draw a segment through a box; return the sorted list of hits on patches.
  // Used for geometric merging of connected components.  Direct return
  // value is true if there is a degeneracy else false.  The sorted list
  // is returned in workspace variable wkspa.hitlist.

  
  bool draw_segment_thru_box(const ActiveBox& b,
    IncidenceTable::Index startpoint_incind,
    int startpoint_relnum,
    const Point& endpoint,
    ActiveBoxVec::Create_Subbox_Workspace& wkspa) {

    wkspa.hitlist.resize(0, ActiveBoxVec::RayHitRecord());

    PatchTable::Index startpoint_patchind = wkspa.inc_table.patchind(startpoint_incind);
    Point startpoint = wkspa.inc_table.real_coord(startpoint_incind);
    int di = wkspa.patchtable.embedded_dim();

    Point normal;
    Point segment_tangent = Point::subtract(endpoint, startpoint);
    
    Real farthest_hit_dist = 0.0;
    pair<PatchTable::Index, PatchTable::ConeIntersectCode> rval1 =
      wkspa.patchtable.determine_ray_front_face(startpoint_patchind, 
      wkspa.inc_table.param_coord(startpoint_incind),
      segment_tangent, wkspa.tol);
    if (rval1.second == PatchTable::DEGEN) {
      if (wkspa.logstr.verbosity() >= 8) {
        wkspa.logstr.str() << " Degeneracy in determine_ray_front_face\n";
      }
      return true;
    }

    PatchTable::Index patchind_s = rval1.first;
    Brep::Face_Spec base_fspec = (rval1.second == PatchTable::HITS_PATCH_FRONT)? 
      wkspa.patchtable.front_side(patchind_s) : wkspa.patchtable.back_side(patchind_s);

    if (wkspa.logstr.verbosity() >= 8)
      wkspa.logstr.str() <<  "reached shootray tangent = " << segment_tangent 
      << " start point = " << startpoint << "= incd#" << wkspa.inc_table.seqno(startpoint_incind)
      << "#"
      <<   "\n basefspec = " << base_fspec << " startpoint_patchind = pat#" 
      << startpoint_patchind << "#\n"
      << " endpoint = " << endpoint << '\n'; 

    ActiveBoxVec::RayHitRecord rayhit;
    rayhit.relnum = startpoint_relnum;
    rayhit.front_hit = base_fspec;
    rayhit.dist = 0.0; 
    rayhit.is_front = rval1.second == PatchTable::HITS_PATCH_FRONT;
    wkspa.hitlist.push_back(rayhit);

    for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
      PatchTable::Index patchind = b.patch(relnum).patchind();      
      if (wkspa.patchtable.gdim(patchind) < di - 1)
        continue;

      if (wkspa.logstr.verbosity() >= 8) {
        wkspa.logstr.str() << "loop relnum = " << relnum
          << "pat#" << patchind << "#" << " orig = " 
          << wkspa.patchtable.owner(patchind) << " # " 
          << wkspa.patchtable.orig_index(patchind) << '\n';
      }

      bool require_startpoint_degen = 
        wkspa.patchtable.is_ancestor(startpoint_patchind, patchind);

      wkspa.output_stack.resize(0, PatchMath::IntersectRecord());
      wkspa.patchtable.patchmath(patchind).intersect_seg(wkspa.output_stack, wkspa.rworkspace, 
        startpoint, endpoint, wkspa.scfac, wkspa.tol, false);      
      if (wkspa.logstr.verbosity() >= 8)
        wkspa.logstr.str() << "require startpoint degen = " << require_startpoint_degen 
        << " size = " << wkspa.output_stack.size() << '\n';

      int sz = wkspa.output_stack.size();
      int degen_count = 0;
      for (int ii = 0; ii < sz; ++ii) {
        if (wkspa.logstr.verbosity() >= 8) {
          wkspa.logstr.str() << "hit " << ii << " at real point " 
            << wkspa.output_stack[ii].realcoord << " param " << wkspa.output_stack[ii].paramcoord
            << "\n degen = " << wkspa.output_stack[ii].degen
            << " uncer = " << wkspa.output_stack[ii].param_uncertainty
            << " seguncer = " << wkspa.output_stack[ii].param_uncertainty_seg
            << " dist = " << wkspa.output_stack[ii].dist << '\n';
        }
        if (wkspa.output_stack[ii].degen) {
          ++degen_count;
          continue;
        }
        normal = 
          wkspa.patchtable.patchmath(patchind).normal(wkspa.output_stack[ii].paramcoord);
        Real ip = Point::inner_product(normal, segment_tangent);
        rayhit.relnum = relnum;
        rayhit.dist = wkspa.output_stack[ii].dist;        
        if (ip > 0) {
          rayhit.front_hit = wkspa.patchtable.front_side(patchind);
          rayhit.back_hit = wkspa.patchtable.back_side(patchind);
          rayhit.is_front = true;
        }
        else {
          rayhit.front_hit = wkspa.patchtable.back_side(patchind);
          rayhit.back_hit = wkspa.patchtable.front_side(patchind);
          rayhit.is_front = false;
        }
        if (wkspa.logstr.verbosity() >= 8) {
          wkspa.logstr.str() << " front = " << rayhit.front_hit << " back = " << rayhit.back_hit << '\n';
        }
        wkspa.hitlist.push_back(rayhit);
      }
      
      if (degen_count > 1 || (degen_count > 0 && !require_startpoint_degen))
        return true;
      if (degen_count == 0 && require_startpoint_degen) {

        wkspa.logstr.str() << "Error! Output stack size = " << sz << '\n';
        {
          for (int jj = 0; jj < sz; ++jj) {
            wkspa.logstr.str() << "output_stack[" << jj << "] dist = "
              << wkspa.output_stack[jj].dist << " degen = "
              << wkspa.output_stack[jj].degen << " paramcoord = "
              << wkspa.output_stack[jj].paramcoord << '\n';
          }
        }
        QMG::Local_to_PatchMath_CPP::debug_str = &wkspa.logstr.str();
        QMG::Local_to_PatchMath_CPP::dump_everything = true;
        wkspa.output_stack.resize(0, PatchMath::IntersectRecord());
        wkspa.patchtable.patchmath(patchind).intersect_seg(wkspa.output_stack, wkspa.rworkspace,
          startpoint, endpoint,
          wkspa.scfac, wkspa.tol, false);
        
        int sz = wkspa.output_stack.size();
        wkspa.logstr.str() << " Second time Output stack size = " << sz << '\n';
        {
          for (int jj = 0; jj < sz; ++jj) {
            wkspa.logstr.str() << "wkspa.output_stack[" << jj << "] dist = "
              << wkspa.output_stack[jj].dist << " degen = "
              << wkspa.output_stack[jj].degen << " paramcoord = "
              << wkspa.output_stack[jj].paramcoord << '\n';
          }
        }
        throw_error("Ray missed starting patch?  This message may indicate a tolerance error");
      }
    }
    sort(wkspa.hitlist.begin(), wkspa.hitlist.end());


#ifdef DEBUGGING
    for (int hitlistind = 0; 
    hitlistind < wkspa.hitlist.size() - 1; 
    hitlistind++) {
      if (wkspa.hitlist[hitlistind].front_hit != wkspa.hitlist[hitlistind+1].back_hit) {
        ostringstream os;
        os << "Problem 2 in draw_segment.\n"
          << "This error may indicate a problem with illegally crossing faces\n"
          << "in the input brep. The two surfaces causing the error are\n"
          << wkspa.patchtable.owner(b.patch(wkspa.hitlist[hitlistind].relnum).patchind())
          << " and "
          << wkspa.patchtable.owner(b.patch(wkspa.hitlist[hitlistind+1].relnum).patchind());
        throw_error(os.str());
      }
    }
#endif
    return false;
  }

    // ------------------------------------------------------------------
  // Merge connected components of the content of a box by shooting rays.
  // Return variable is wkspa.uf (union find).
  
  void segshoot_merge(const ActiveBox& b,
    ActiveBoxVec::Create_Subbox_Workspace& wkspa) {
    
    int di = wkspa.patchtable.embedded_dim();
    int ufsize = wkspa.uf.size();     
    wkspa.ray_shot_from_label.resize(ufsize);
    for (int jj = 0; jj < ufsize; ++jj)
      wkspa.ray_shot_from_label[jj] = false;

    bool degenhit = true;
    int trycount = 0;
    if (wkspa.logstr.verbosity() >= 8) {
      wkspa.logstr.str() << " reached segshoot merge box = " << b << '\n';
    }

    int oldverb = wkspa.logstr.verbosity();

    while (degenhit && ++trycount < 400) {
      Point randpoint = b.select_rand_point();
      if (trycount == 398) {        
        wkspa.logstr.change_verbosity(8);
        wkspa.logstr.str() << "trycount at 398 b = ";
        b.full_dump(wkspa.logstr.str());
      }
      if (wkspa.logstr.verbosity() >= 8) {
        wkspa.logstr.str() << "trycount = " << trycount 
          << " randpoint = " << randpoint << "\n";
      }

      degenhit = false;
      int root_component_patchrelnum = -1;

      for (int relnum1 = 0; relnum1 < b.num_patch(); ++relnum1) 
        wkspa.ray_shot_from_label[relnum1] = false;


      for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
        PatchTable::Index patchind = b.patch(relnum).patchind();
        int c = wkspa.uf.get_label(relnum);
        if (wkspa.logstr.verbosity() >= 8) {
          wkspa.logstr.str() << " in loop relnum = " << relnum << " c = " << c
            << " rsfl = " << wkspa.ray_shot_from_label[c] << '\n';
        }      

        if (!wkspa.ray_shot_from_label[c]) {

          // Select starting point on patch relnum.

          Patch_in_ActiveBox ptc = b.patch(relnum);
          if (ptc.num_incidence() == 0) 
            throw_error("No incidences on patch of component?");
          IncidenceTable::Index incind = ptc.incidence(0);
          if (wkspa.logstr.verbosity() >= 8) {
            wkspa.logstr.str() << " calling draw seg from inc incd#" 
              << wkspa.inc_table.seqno(incind) << "#\n";
          }

          degenhit = draw_segment_thru_box(b, incind, relnum, randpoint, wkspa);
          
          // Last call puts return list in wkspa.hitlist.

          if (degenhit) {
            if (wkspa.logstr.verbosity() >= 8) {
               wkspa.logstr.str() << "degenhit = true\n";
            }
            break;
          }
#ifdef DEBUGGING
          // First hit should be on p.
          if (wkspa.hitlist.size() == 0 || wkspa.hitlist[0].relnum != relnum) {
            throw_error("Problem 1 in rayshoot_merge");
            return;
          }
#endif        

          if (wkspa.logstr.verbosity() >= 8) {
            wkspa.logstr.str() << " hitlist size = " << wkspa.hitlist.size() << '\n';
            for (int iii = 0; iii < wkspa.hitlist.size(); ++iii) {
              wkspa.logstr.str() << "hit list " << iii << " dist = " << wkspa.hitlist[iii].dist 
                << " relnum = " <<  wkspa.hitlist[iii].relnum
                << " front = " << wkspa.hitlist[iii].front_hit << " back = " 
                << wkspa.hitlist[iii].back_hit << '\n';
              
            }
          }
          // Merge each item in the hit list with its successor.

          for (int hitlistind = 0; 
          hitlistind < wkspa.hitlist.size() - 1; 
          hitlistind++) {

            int relnum1 = wkspa.hitlist[hitlistind].relnum;
            int relnum2 = wkspa.hitlist[hitlistind + 1].relnum;