qseparation.cpp

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

  }

  // Loop over patches and evaluate the curve-control function
  // to set the rqvariance_ field.
  
  {
    Point incpoint_real;
    Point incpoint_par;
    for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
      Patch_in_ActiveBox patch = b.patch(relnum);
      PatchTable::Index patchind = patch.patchind();      
      Brep::Face_Spec owner = patchtable_.owner(patchind);
      if (owner.fdim() == 0)
        continue;      
      int numinc = patch.num_incidence();
      for (int ii = 0; ii < numinc; ii++) {
        IncidenceTable::Index inc = patch.incidence(ii);
        incpoint_real = inc_table_.real_coord(inc);
        incpoint_par = inc_table_.param_coord(inc);
        pair<Real,bool> rval = 
          curve_control_.eval(incpoint_real, incpoint_par, patchind);
        Real rqvariance1 = rval.first;
        if (rqvariance_(owner) > rqvariance1)
          rqvariance_(owner) = rqvariance1;
        if (logstr_.verbosity() >= 10) {
          logstr_.str() << " rqvariance(" << owner << ") = " << rqvariance_(owner) 
            << '\n';
        }
      }
    }
  }


  // Initialize the track_max_min map.
  
  {
    int old_serial_num = serial_num_;
    ++serial_num_;

    // check for overflow in serial_num_
    if (serial_num_ <= old_serial_num) {
      serial_num_ = 1;
      for (Brep::Face_Spec_Loop_Over_Faces fspec(brep_);
      fspec.notdone();
      ++fspec) {
        max_min_initialized_(fspec) = 0;
        curvature_is_checked_(fspec) = 0;
      }
    }

    for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
      PatchTable::Index patchind = b.patch(relnum).patchind();
      Brep::Face_Spec owner = patchtable_.owner(patchind);
      if (max_min_initialized_(owner) == serial_num_) continue;
      max_min_initialized_(owner) = serial_num_;
      for (int l = 0; l < num_test_directions_actual_; ++l) {
        track_max_min_(owner).lower_bound[l] = 100;
        track_max_min_(owner).upper_bound[l] = -100;
      }
    }
  }
 
  // Figure out the variation of the normal on each topological
  // face in the box.  Use a heuristics: find the maximum variation
  // in each direction in test_directions.
  
  {
    for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
      PatchTable::Index patchind = b.patch(relnum).patchind();
      Brep::Face_Spec owner = patchtable_.owner(patchind);
      int pdim = patchtable_.gdim(patchind);
      if (logstr_.verbosity() >= 9) {
        logstr_.str() << "relnum = " << relnum << " patchind = "
          << patchind << " owner = " << owner << '\n';      
      }
      int numinc = b.patch(relnum).num_incidence();
      for (int ii = 0; ii < numinc; ++ii) {
        IncidenceTable::Index incind = b.patch(relnum).incidence(ii);
        Point paramcoord = inc_table_.param_coord(incind);
        Point realcoord = inc_table_.real_coord(incind);
        process_incidence1_(patchind, paramcoord, realcoord);
        if (pdim == di_ - 1) continue;
        for (PatchTable::Loop_over_parents loop(patchtable_, patchind);
        loop.notdone();
        ++loop) {
          PatchTable::Index parentind = loop.parent_index();
          Point paramcoord1 = patchtable_.lift_parameters(paramcoord, 
            patchind, loop.listind());
          process_incidence1_(parentind, paramcoord1, realcoord);
          if (pdim == di_ - 2) continue;
          for (PatchTable::Loop_over_parents loop2(patchtable_, parentind);
          loop2.notdone();
          ++loop2) {
            PatchTable::Index grandparentind = loop2.parent_index();
            Point paramcoord2 = patchtable_.lift_parameters(paramcoord1, 
              parentind, loop2.listind());
            process_incidence1_(grandparentind, paramcoord2, realcoord);
          }
        }
      }
    }
  }
      
  // Loop over patches and compare the variation of the normal over
  // the patch to the curve-control function.  
  
  Size_Return ret;
  {
    Point incpoint_real;
    Point incpoint_par;
    for (int relnum = 0; relnum < b.num_patch(); relnum++) {
      Patch_in_ActiveBox patch = b.patch(relnum);
      PatchTable::Index patchind = patch.patchind();      
      Brep::Face_Spec owner = patchtable_.owner(patchind);
      if (owner.fdim() == 0)
        continue;      
      if (curvature_is_checked_(owner) == serial_num_) continue;
      curvature_is_checked_(owner) = serial_num_;
      Real rq1 = rqvariance_(owner);
      for (int l = 0; l < num_test_directions_actual_; ++l) {
        Real delta = track_max_min_(owner).upper_bound[l] -
          track_max_min_(owner).lower_bound[l];
        if (logstr_.verbosity() >= 14) {
          logstr_.str() << "pat#" << patchind << "# owner " << owner
            << " dir " << l << " delta = " << delta << '\n';
        }
        if (delta > rq1) {
          if (logstr_.verbosity() >= 8)
            logstr_.str() << " curvature not ok rqvariance = " << rq1 << " variation on "
            << owner << " is " << delta << " dir "  << l<< " crv-return false\n";
          ret.fspec = owner;
          ret.val = rq1;
          ret.ok = false;
          logstr_.change_verbosity(oldverb);
          return ret;
        }
      }
    }
  }

  // All the stuff in the box should be connected, either directly
  // or by using patches outside the box by a factor of 2.
  // The farther-away patches must satisfy the curvature rqt.

  patch_searched_.reset(My_bool(false));
  acceptable_for_search_cache_.reset(-1);
  map_to_relnum_.reset(-1);
  {
    for (int relnum = 0; relnum < b.num_patch(); ++relnum)
      map_to_relnum_[b.patch(relnum).patchind()] = relnum;
  }
  Point bmidpoint;
  for (int j = 0; j < di_; ++j) {
    if (b.flatdim(j)) {
      bmidpoint[j] = b.real_lowerleft()[j];
    }
    else {
      bmidpoint[j] = b.real_lowerleft()[j] + boxwidth;
    }
  }

  patchind_stack_.resize(0);
  patchind_stack_.push_back(b.patch(0).patchind());
  while (patchind_stack_.size() > 0) {
    PatchTable::Index newpatch = patchind_stack_.back();
    patchind_stack_.pop_back();
    if (patch_searched_[newpatch]) continue;
    if (logstr_.verbosity() >= 10) {
      logstr_.str() << " searching pat#" << newpatch << "#\n";
    }
    patch_searched_[newpatch] = true;
    for (PatchTable::Loop_over_children cloop(patchtable_, newpatch);
    cloop.notdone();
    ++cloop) {
      PatchTable::Index pind1 = cloop.child_index();
      if (acceptable_for_search_(pind1, bmidpoint, boxwidth))
        patchind_stack_.push_back(pind1);
    }
    for (PatchTable::Loop_over_parents ploop(patchtable_, newpatch);
    ploop.notdone();
    ++ploop) {
      PatchTable::Index pind1 = ploop.parent_index();
      if (acceptable_for_search_(pind1, bmidpoint, boxwidth))
        patchind_stack_.push_back(pind1);
    }
  }

  // Find if there is a patch that was not reached by the preceding search.
  for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
    if (!patch_searched_[b.patch(relnum).patchind()]) {      
      if (logstr_.verbosity() >= 8)
        logstr_.str() << " curvature not ok disconnection on patch pat#"
        << b.patch(relnum).patchind() << "# relnum " << relnum
        << " from root pat#" << b.patch(0).patchind() << "# relnum " << relnum
        << " crv-return false\n";
      ret.fspec = patchtable_.owner(b.patch(relnum).patchind());
      ret.val = rqvariance_(ret.fspec);
      ret.ok = false;
      logstr_.change_verbosity(oldverb);
      return ret;
    }
  }
  if (logstr_.verbosity() >= 8)
    logstr_.str() << " crv-return true";
  ret.ok = true;
  logstr_.change_verbosity(oldverb);
  return ret;
}



QMG::MG::Separation_Functor::Separation_Functor(const Brep& brep,
                                                const PatchTable& patchtable,
                                                IncidenceTable& inc_table,
                                                const SizeControl& size_control,
                                                const SizeControl& curve_control,
                                                Brep_Orbits& orbits,
                                                double scfac,
                                                double tol,
                                                Logstream& logstr,
                                                Meshgen_gui& gui) :
brep_(brep),
logstr_(logstr),
patchtable_(patchtable),
inc_table_(inc_table),
size_control_(size_control),
curve_control_(curve_control),
orbits_(orbits),
scfac_(scfac),
tol_(tol),
di_(brep.embedded_dim()),
lca_(brep),
track_max_min_(brep),
num_test_directions_actual_((di_ == 2)? NUM_TEST_DIRECTIONS_2D:NUM_CURVATURE_TEST_DIRECTION),
test_directions_(num_test_directions_actual_, Point()),
patch_searched_(patchtable, My_bool(false)),
acceptable_for_search_cache_(patchtable, 0),
map_to_relnum_(patchtable,0),
rqvariance_(brep),
serial_num_(1),
max_min_initialized_(brep),
curvature_is_checked_(brep),
crt_wkspa_(patchtable, inc_table, logstr, gui, scfac, tol)
{
  for (int i = 0; i < num_test_directions_actual_; ++i) {
    for (int j = 0; j < di_; ++j) {
      test_directions_[i][j] = test_direction_init[i*3+j];
    }
    test_directions_[i].normalize();
  }
}
  


// Process a completely interior box for separation.  The only reason to
// split such a box is the size function.

void QMG::MG::Separation_Functor::completely_interior(NonConstActiveBox& b,
                                                      ActiveBoxVec& smallboxvec)  {

  int di = b.embedded_dim();
  
#ifdef DEBUGGING
  if (!b.completely_interior()) {
    throw_error("Wrong separation function called 1.");
    return;
  }
#endif
  if (logstr_.verbosity() >= 2) 
    logstr_.str() << "separation comp int on box " << b;  
  bool size_verified = b.size_function_verified();
  int split_reason = 0;
  Brep::Face_Spec apex(di, b.owner());
  Size_Return ret;
  if (!size_verified) {
    ret = ok_for_size_function_completely_interior_(b, apex);
    if (ret.ok)
      b.mark_size_verified();
    else
      split_reason = 2;
  }
  if (split_reason) {
    if (logstr_.verbosity() >= 2)
      logstr_.str() << " split reason = " << split_reason << '\n';
    if (b.iwidth() >= MAXLEV) {
      ostringstream os;
      os << "Quadtree has been split too finely because of the size"
       << " control function on region " << brep_.face_name(ret.fspec)
        << ".\nCould not match requested size of " << ret.val;
      throw_error(os.str());
    }
    ActiveBoxVec::split_completely_interior_box_and_push_children(b, smallboxvec,
    logstr_);
  }
  else {
    if (logstr_.verbosity() >= 2)
      logstr_.str() << " put in orbit of " << apex.fdim() << ":" << apex.faceind() << '\n';
    orbits_.retrieve_orbit(apex, b.iwidth()).push_box_completely_interior(b, logstr_);
  }
}


// Process a box during the separation stage.  Box not completely interior.

void QMG::MG::Separation_Functor::
not_interior(NonConstActiveBox& b,
           int phase,
           ActiveBoxVec& tempstack,
           ActiveBoxVec& smallboxstack,
           ActiveBoxVec& idlestack,
           ActiveBoxVec& smallboxidlestack)  {

#ifdef DEBUGGING
  if (b.completely_interior())
    throw_error("Wrong separation function called 2.");
#endif
  
  if (logstr_.verbosity() >= 2) { 
    logstr_.str() << "separation on box ";
    if (logstr_.verbosity() >= 4)
      b.full_dump(logstr_.str());
    else
      logstr_.str() << b;
  }
 
  // First test: is it crowded?
  
  Triple<Brep::Face_Spec, Brep::Face_Spec, bool> returnval = is_crowded_(b, phase);
  
  bool crowded = returnval.third;
  Brep::Face_Spec apex = returnval.first;
  Brep::Face_Spec apex2 = returnval.second;
  
  bool size_verified = b.size_function_verified();
  bool curvature_verified = b.curvature_verified();
  
  Size_Return size_ret, curve_ret;
  
  if (apex.fdim() > phase) {
    if (logstr_.verbosity() >= 2)
      logstr_.str() << " no apex (idle)\n";
    idlestack.push_box_shrink_ex(b, crt_wkspa_);
  }
  else {
    int split_reason = 0;
    if (crowded) {
      split_reason = 1;
    }

    // second test: is it too big for the size function?
    
    if (!split_reason && !size_verified) {
      size_ret = ok_for_size_function_(b, apex);
      if (size_ret.ok)
        b.mark_size_verified();
      else
        split_reason = 2;
    }
    // third test: is it too curvy for the curve control function?
      
    if (!split_reason && !curvature_verified) {
      curve_ret = ok_for_curvature_(b);
      if (curve_ret.ok)
        b.mark_curvature_verified(); 
      else
        split_reason = 3;
    }
    if (split_reason) {
      if (logstr_.verbosity() >= 2) {
        logstr_.str() << " split reason " << split_reason << '\n';
        if (split_reason == 1) {
          Brep::Face_Spec apex2 = returnval.second;
          logstr_.str() << " crowding caused by " << apex.fdim() << ":" << apex.faceind()
            << " and " << apex2.fdim() << ":" << apex2.faceind() << '\n';
        }
        else if (split_reason == 2) {
          logstr_.str() << " size function failed on facespec " << size_ret.fspec
            << '\n';
        }
        else if (split_reason == 3) {
          logstr_.str() << " curvature function failed on facespec " << curve_ret.fspec
            << "\n";
        }
      }
      if (b.iwidth() >= MAXLEV || 1.0 / static_cast<double>(1 << b.iwidth()) < 100 * tol_) {
        ostringstream os;
        os << "Quadtree has been split too finely, possibly indicating an error\n"
          << "in the input brep.  The cause of the last split is:\n";
        if (split_reason == 1) {
          os << "Unable to separate brep faces " << brep_.face_name(apex) 
            << " from " << brep_.face_name(apex2);
        }
        else if (split_reason == 2) {
          os << "Unable to match desired size function value of "
            << size_ret.val << " on topological entity " 
            << brep_.face_name(size_ret.fspec);
        }
        else if (split_reason == 3) {
          os << "Unable to match curvature value of " << curve_ret.val
            << " on topological entity " << brep_.face_name(curve_ret.fspec);
          os << " Patches involved are: ";
          for (int relnum = 0; relnum < b.num_patch(); ++relnum) {
            PatchTable::Index patchind = b.patch(relnum).patchind();
            if (patchtable_.owner(patchind) != curve_ret.fspec) continue;
            if (patchtable_.gdim(patchind) != curve_ret.fspec.fdim()) continue;
            os << patchtable_.orig_index(patchind) << " ";
          }
        }
        throw_error(os.str());
      }

      ActiveBoxVec::split_box_and_push_children(phase, tempstack, smallboxstack, 
        smallboxidlestack, b, inc_table_, crt_wkspa_);
    }
    else {
      if (logstr_.verbosity() >= 2)
        logstr_.str() << " moved to orbit of " << apex.fdim() << ":" << apex.faceind() << '\n';
      orbits_.retrieve_orbit(apex, b.iwidth()).push_box_shrink_ex(b, crt_wkspa_);
    }
  }
}