align.cpp

算断裂的

  typedef map<DeduplicationMapKey, ActiveBoxVec::Index> DeduplicationMap;
  

}

void QMG::MG::alignInitStaticData(int di) {
  PriorityQueueKey_lastphase::di = di;
}

#ifdef NSF_DEMO
namespace QMG {
  namespace Nsf_Demo {
    using namespace QMG;
    using namespace QMG::MG;
    extern void show_alignment(const ActiveBoxVec& vec,
      const vector<ActiveBoxVec::Index>& boxindvec,
      const Point& alignpoint);
  }
}
#endif


// ------------------------------------------------------------------
// align_orbit
// This routine aligns an orbit of a brep.  First it puts all the
// closefaces in a priority queue.  Then boxes whose closefaces
// are completely covered are protected, their subfaces launched,
// and a new vertex appears in the simplicial complex.




void QMG::MG::align_orbit(ActiveBoxVec& alignvec,
                          int closefacedim,
                          Real maxwarpdist,
                          ActiveBoxVec& idlevec,
                          ActiveBoxVec& smallboxvec,
                          ActiveBoxVec& small_idle_boxvec,
                          const Brep::Face_Spec& alignfspec,
                          BoxFaceDag& dag,
                          SimpComplex_Under_Construction& sc,
                          vector<Brep::Face_Spec>& vtx_owner,
                          const vector<Real>& asp_degrade,
                          Logstream& logstr,
                          Meshgen_gui& gui,
                          const PatchTable& patchtable, 
                          IncidenceTable& inc_table,
                          ActiveBoxVec::Create_Subbox_Workspace& wkspa) {

  debug_ostream = &logstr.str();

  PriorityQueueKey::base_idx = 0;
  

  Real bxwidth = BoxAddress::multiplier() / ( (Real)(1 << alignvec.iwidth()) );
  Real maxwarpdist_act = maxwarpdist * bxwidth;
  
  IncidencePriorityQueue ipq(closefacedim, maxwarpdist_act, alignfspec,
    patchtable, inc_table, logstr);

  // Put box closefaces into the queue.

  {
    for (ActiveBoxVec::Loop_over_boxes_nonconst loop(alignvec);
    loop.notdone();
    ++loop) {
      NonConstActiveBox b = loop.getbox_nonconst();
#ifdef DEBUGGING
      if (b.completely_interior()) {
        throw_error("Completely interior box in wrong orbit?");
        return;
      }
#endif
      ipq.insert_box_closefaces(b);
      b.assign_close_face(false);
      gui.check_update();
    }
  }
  
  vector<PriorityQueueKey> boxes_to_align(0);
#ifdef NSF_DEMO
  vector<ActiveBoxVec::Index> boxind_vec;
#endif

  bool too_far_detected = false;
  int di = patchtable.embedded_dim();
  int phase = alignfspec.fdim();
  vector<SimpComplex::VertexOrdinalIndex> volist(di + 1);
  vector<Real> warp_point_re(di);
  vector<Real> warp_point_pa(alignfspec.fdim());
  Point warp_point_pa0;

  // Loop over the priority queue, either until the distance gets 
  // so big that an unacceptable warp is found, or the queue is empty.

  while (!too_far_detected && ipq.size() > 0) {
    
    gui.check_update();
    
    
    // Copy all the highest priority box faces with the same incidence index
    // to stack boxes_to_align.

    boxes_to_align.resize(0);
    Mask flatdim1 = ipq.top().flatdim;
    IncidenceTable::Index incind1 = ipq.top().incind;
    do {
      if (!(alignvec.iskilled(ipq.top().boxind)))
        boxes_to_align.push_back(ipq.top());
      ipq.pop();
    } while (ipq.size() > 0 && ipq.top().flatdim == flatdim1 &&
      ipq.top().incind == incind1);
    int basize = boxes_to_align.size();
    if (basize == 0)
      continue;

    if (logstr.verbosity() >= 2) { ////
      logstr.str() << "Attempting to align boxes\n";
      for (int jj = 0; jj < basize; ++jj) {
        ActiveBoxVec::Index boxind = boxes_to_align[jj].boxind;
        const ActiveBox& b1 = alignvec.getbox(boxind);
        logstr.str() << "    " << b1 << " ";
        boxes_to_align[jj].rindex.output(logstr.str(), b1.dim());
      }
      logstr.str() << " to incidence incd#" << inc_table.seqno(incind1) << "#\n";
    }

    
    // Test all of these boxes for OK aspect ratio.
    
    Real dist1 = boxes_to_align[0].dist;
    if (dist1 == BIG_REAL) {
      if (logstr.verbosity() >= 2) { ////
        logstr.str() << " distance = " << dist1 << " too far\n";
      }
      too_far_detected = true;
      break;
    }


    // Check if we are trying to align a too low-dimensional box face.

    {
      for (int jj = 0; jj < basize; ++jj) {
        ActiveBoxVec::Index boxind = boxes_to_align[jj].boxind;
        const ActiveBox& b1 = alignvec.getbox(boxind);
        if (b1.dim() + phase < di) {
          too_far_detected = true;
#ifdef DEBUGGING
          if (dist1 <= maxwarpdist_act * .05 && logstr.verbosity() >= 1) {
            logstr.str() << "Warning: warping anomaly may lead to poor aspect ratio.\n"
              << "Suggest you re-run QMG with smaller curvature setting.";
          }
#endif
          break;
        }
      }
    }

    if (too_far_detected)
      break;

    Point warp_point = inc_table.real_coord(incind1);
    if (dist1 > 0.0) {
      for (int jj = 0; jj < basize; ++jj) {
        ActiveBox b = alignvec.getbox(boxes_to_align[jj].boxind);
        too_far_detected = !warp_ok(dist1, warp_point, b,
          boxes_to_align[jj].rindex, dag, b.link_parent(), sc, 
          asp_degrade[b.dim()], logstr, wkspa.scfac, wkspa.tol, volist);
        if (too_far_detected) {
          if (logstr.verbosity() >= 2) { ////
            logstr.str() << " distance = " << dist1 << " too far detected on " << jj << "th box\n";
          }
          break;
        }
      }
    }
    
    if (too_far_detected)
      break;
      
    // Check if we can align these boxes by adding up their weights.
    
    Weight_t sumweight = 0.0;
    for (int kk = 0; kk < basize; ++kk) {
      sumweight += alignvec.getbox(boxes_to_align[kk].boxind).weight(boxes_to_align[kk].rindex);
    }
    
    if (logstr.verbosity() >= 2) ////
      logstr.str() << "weightsum = " << sumweight << '\n';
    if (sumweight > 1.00001)
      throw_error("Excess weight sum");

    // If the weight sum is 1, protect all the boxes, launch subfaces,
    // and put a new vertex in the simplicial complex.

    if (sumweight > 0.99999) {
      PatchTable::Index patchind = inc_table.patchind(incind1);
      warp_point_pa0 = inc_table.param_coord(incind1);
      for (int ll = 0; ll < di; ++ll)
        warp_point_re[ll] = warp_point[ll];
      for (int mm = 0; mm < alignfspec.fdim(); ++mm)
        warp_point_pa[mm] = warp_point_pa0[mm];
      pair<SimpComplex::VertexOrdinalIndex, SimpComplex::VertexGlobalIndex> rval = 
        sc.add_vertex(warp_point_re);
      SimpComplex::VertexOrdinalIndex vertnum_o = rval.first;
      SimpComplex::VertexGlobalIndex vertnum = rval.second;
      sc.add_vertex_bdryinc(vertnum, alignfspec,
        patchtable.orig_index(patchind), warp_point_pa);
      if (vtx_owner.size() <= vertnum_o) {
        vtx_owner.resize(vertnum_o + 1, alignfspec);
      }
      else 
        vtx_owner[vertnum_o] = alignfspec;
      if (logstr.verbosity() >= 2) { ////
        logstr.str() << "Alignment successful; new vertex vtx#" << vertnum 
          << "#   ord#" << vertnum_o << " at coords " << warp_point << " on patch " << patchind << '\n';
      }
#ifdef NSF_DEMO
      boxind_vec.resize(basize);
      for (int jj2 = 0; jj2 < basize; ++jj2) {
        boxind_vec[jj2] = boxes_to_align[jj2].boxind;
      }
      QMG::Nsf_Demo::show_alignment(alignvec, boxind_vec, warp_point);
#endif
      for (int jj = 0; jj < basize; ++jj) {
        ActiveBoxVec::Index boxind = boxes_to_align[jj].boxind;
        ActiveBox b = alignvec.getbox(boxind);
        Real maxwarpdist_act = maxwarpdist * b.real_width();
        Base3 rindex = boxes_to_align[jj].rindex;
        protect_box_and_launch_subfaces(b, vertnum_o, rindex,
          alignfspec, alignvec, idlevec, dag, ipq, wkspa);
        alignvec.killbox(boxind, logstr);
      }

    } 
    else {
      // If the weight sum is < 1, then mark the boxes to indicate that they
      // had a close face.  This forces them to be split.

      for (int jj = 0; jj < basize; ++jj) {
        alignvec.getbox_nonconst(boxes_to_align[jj].boxind).assign_close_face(true);
        // check_containment(alignvec.getbox(boxes_to_align[jj].boxind), idlevec, logstr.str());
      }
    }
  }
  
  
  // All the boxes that have a close subface and weren't aligned must be split.
  // Must split if this is the last pass.

  bool last_pass = closefacedim == di - phase;
  {
    for (ActiveBoxVec::Loop_over_boxes_nonconst loop(alignvec);
    loop.notdone();
    ++loop) {
      gui.check_update();
      ActiveBox b = loop.getbox();
      if (b.has_close_face() || last_pass) {
        if (b.iwidth() >= MAXLEV) 
          throw_error("Oversplitting for alignment (1)");
        ActiveBoxVec::split_box_and_push_children(phase, alignvec, 
          smallboxvec, small_idle_boxvec, b,  inc_table,  wkspa);
        alignvec.killbox(loop.index(), logstr);
      }
    }
  }
}



// ------------------------------------------------------------------
// align_orbit_lastphase
// This aligns boxes during the last phase.  It's quite similar
// to align, except there are some simplifications because
// there is no warping.



void QMG::MG::align_orbit_lastphase(ActiveBoxVec& alignvec,
                                    ActiveBoxVec& smallboxvec,
                                    const Brep::Face_Spec& alignfspec,
                                    BoxFaceDag& dag,
                                    SimpComplex_Under_Construction& sc,
                                    vector<Brep::Face_Spec>& vtx_owner,
                                    Logstream& logstr,
                                    Meshgen_gui& gui) {
  
  
  PriorityQueueKey_lastphase::base_idx = 0;
  
  IncidencePriorityQueue_LastPhase ipq(logstr);
  int di = alignfspec.fdim();

  // Put all the vertices of the boxes into the priority queue.

  {
    for (ActiveBoxVec::Loop_over_boxes loop(alignvec);
    loop.notdone();
    ++loop) {
      gui.check_update();
      ActiveBox b = loop.getbox();
#ifdef DEBUGGING
      if (!b.completely_interior()) {
        throw_error("Not completely interior box in wrong orbit");
        return;
      }
#endif
      ipq.insert_box_closefaces(b);
    }
  }
  
  vector<PriorityQueueKey_lastphase> boxes_to_align(0);
#ifdef NSF_DEMO
  vector<ActiveBoxVec::Index> boxind_vec;
#endif

  BoxAddress bxa;
  vector<Real> new_point_re(di);
  vector<Real> new_point_pa(0);

  while (ipq.size() > 0) {
    
    // Copy all the highest priority box faces with the same incidence index
    // to stack boxes_to_align.

    gui.check_update();
    
    boxes_to_align.resize(0);
    bxa.lowerleft() = ipq.top().position;
    do {
      if (!alignvec.iskilled(ipq.top().boxind))
        boxes_to_align.push_back(ipq.top());
      ipq.pop();
    } while (ipq.size() > 0 && ipq.top().position == bxa.lowerleft());
    int basize = boxes_to_align.size();
    if (basize == 0)
      continue;

    if (logstr.verbosity() >= 2) { ////
      logstr.str() << "Attempting to (last phase) align boxes\n";
      for (int jj = 0; jj < basize; ++jj) {
        ActiveBoxVec::Index boxind = boxes_to_align[jj].boxind;
        const ActiveBox& b1 = alignvec.getbox(boxind);
        logstr.str() << "    " << b1 << " ";
        boxes_to_align[jj].rindex.output(logstr.str(), b1.dim());
        logstr.str() << '\n';
      }
    }
    
    // Check if we can align these boxes by adding up their weights.
    
    Weight_t sumweight = 0.0;
    for (int jj = 0; jj < basize; ++jj) {
      sumweight += alignvec.getbox(boxes_to_align[jj].boxind).weight(boxes_to_align[jj].rindex);
    }
        
    if (logstr.verbosity() >= 2) ////
      logstr.str() << "weightsum = " << sumweight << '\n';
    if (sumweight > 1.00001)
      throw_error("Excess weight sum last phase");
    if (sumweight > 0.99999) {
      Point pt = bxa.real_lowerleft();
      for (int ll = 0; ll < di; ++ll)
        new_point_re[ll] = pt[ll];

      pair<SimpComplex::VertexOrdinalIndex, SimpComplex::VertexGlobalIndex> rval = 
        sc.add_vertex(new_point_re);

      SimpComplex::VertexOrdinalIndex vertnum_o = rval.first;
      SimpComplex::VertexGlobalIndex vertnum = rval.second;
      if (vtx_owner.size() <= vertnum_o) {
        vtx_owner.resize(vertnum_o + 1, alignfspec);
      }
      else 
        vtx_owner[vertnum_o] = alignfspec;
      if (logstr.verbosity() >= 2) { ////
        logstr.str() << "Alignment successful; new vertex vtx#" << vertnum 
          << "# ord #" << vertnum_o << " at coords " << pt << '\n';
      }
#ifdef NSF_DEMO
      boxind_vec.resize(basize);
      for (int jj2 = 0; jj2 < basize; ++jj2) {
        boxind_vec[jj2] = boxes_to_align[jj2].boxind;
      }
      QMG::Nsf_Demo::show_alignment(alignvec, boxind_vec, pt);
#endif
      for (int jj = 0; jj < basize; ++jj) {
        int boxind = boxes_to_align[jj].boxind;
        ActiveBox b = alignvec.getbox(boxind);
        protect_box_and_launch_subfaces_lastphase(b, vertnum_o,
          boxes_to_align[jj].rindex, alignfspec, alignvec, dag, ipq, logstr);
        alignvec.killbox(boxind, logstr);
      }
    } 
  }
  
  // All the boxes that have a close subface and weren't aligned must be split.
  
  {
    for (ActiveBoxVec::Loop_over_boxes loop(alignvec);
    loop.notdone();
    ++loop) {
      gui.check_update();
      ActiveBox b = loop.getbox();
      if (b.iwidth() >= MAXLEV)
        throw_error("Oversplitting for alignment (2)");
      ActiveBoxVec::split_completely_interior_box_and_push_children(b,
        smallboxvec, logstr);
      alignvec.killbox(loop.index(), logstr);
    }
  }
}