meshgen.cpp

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

          b.full_dump2(*debug_ostream, wkspa.patchtable, wkspa.inc_table);
        else
          *debug_ostream << b;
        *debug_ostream << " face ";
        c_rindex.output(*debug_ostream, b.dim());
        *debug_ostream << " with weight " << b.weight(c_rindex);
        ActiveBoxVec tmpvec(.01,iw);
        ActiveBoxVec::create_subbox(ActiveBox_As_Parent(b), bxa, iw,
          di, 0, 0, 0, tmpvec, tmpvec, wkspa);
        if (tmpvec.numboxes() > 0) {
          *debug_ostream << " actually inside\n";
          sum += b.weight(c_rindex);
        }
        else {
          *debug_ostream << " not inside\n";
        }
      }
    }
    return sum;
  }
  
  Weight_t sum_weight_bvl(Boxvecs_by_level& sep,
                          const IntCoord& coord,
                          int iw,
                          ActiveBoxVec::Create_Subbox_Workspace& wkspa) {
    Weight_t sum = 0;
    for (int level = 0; level < sep.size(); ++level) {
      ActiveBoxVec* vecp = sep[level];
      if (vecp)
        sum += sum_weight_stack(*vecp, coord, iw, wkspa);
    }
    return sum;
  }

  Weight_t sum_weight_orbit(Brep_Orbits& all_orbits,
                            const Brep::Face_Spec& fspec,
                            const IntCoord& coord,
                            int iw,
                            ActiveBoxVec::Create_Subbox_Workspace& wkspa){
    int level = 0;
    ActiveBoxVec* vecp = all_orbits.find_next_nonempty_level(fspec, level);
    Weight_t sum = 0;
    while (vecp) {
      level = vecp -> iwidth();
      sum += sum_weight_stack(*vecp, coord, iw, wkspa);
      vecp = all_orbits.find_next_nonempty_level(fspec, level+1);
    }
    return sum;
  } 

  void detect_weight_problem_at_point(Boxvecs_by_level& sep0, 
                                      Boxvecs_by_level& sep1, 
                                      Brep_Orbits& all_orbits,
                                      const Brep& g,
                                      const IntCoord& coord,
                                      int testiw,
                                      ActiveBoxVec::Create_Subbox_Workspace& wkspa,
                                      bool throw_if_error) {
    

    *debug_ostream << "Coord = ";
    for (int i1 = 0; i1 < g.embedded_dim(); ++i1) {
      *debug_ostream << std::hex << coord[i1] << " " << std::dec;
    }
    *debug_ostream << "\nreal coord = ";
    for (int i2 = 0; i2 < g.embedded_dim(); ++i2) {
      *debug_ostream << std::setprecision(17) 
        << coord[i2] * BoxAddress::multiplier() / (1 << MAXLEV)
        + BoxAddress::global_coordinate_origin(i2) << " ";
    }
    *debug_ostream << '\n';
  


    Weight_t sum = 0;
    *debug_ostream << " weight test sep 0\n";
    sum += sum_weight_bvl(sep0, coord, testiw, wkspa);
    *debug_ostream << " weight test sep 1\n";
    sum += sum_weight_bvl(sep1, coord, testiw, wkspa);
    for (Brep::Face_Spec_Loop_Over_Faces fspec(g); fspec.notdone(); ++fspec) {
      *debug_ostream << " weight test orbit " << fspec << "\n";
      sum += sum_weight_orbit(all_orbits, fspec, coord, testiw, wkspa);
    }
    *debug_ostream << " weight sum in detect = " << sum << '\n';
    if (sum  < .99999 || sum > 1.00001) {
      *debug_ostream << " illegal sum = " << sum << std::endl;
      if (throw_if_error)
        throw_error("Weight error detected");
    }
  }


  // ------------------------------------------------------------------
  // This class exists to cause an error message to be written to the log,
  // and then for the log to flushed whenever an error is thrown.  

  class OnErrorFlushLog : public QMG::Error_Task {
  private:
    Logstream& log_;
    // Override base case function
    virtual void on_error_do(const string& errmsg) { 
      log_.str() << errmsg << '\n';
      QMG::Error_Message::put_error_messages_into_stream(log_.str());
      log_.str() << std::flush;
    }
  public:
    OnErrorFlushLog(Logstream& logs) : log_(logs) { }
    ~OnErrorFlushLog() { }
  };

  // ------------------------------------------------------------------
  // This class creates a button labeled "dismiss" when an error
  // occurs during mesh generation.

  class OnErrorCreateDismissButton : public QMG::Error_Task {
  private:
    Meshgen_gui& gui_;
    virtual void on_error_do(const string& errmsg) {
      gui_.indicate_error_has_occurred(errmsg);
    }
  public:
    OnErrorCreateDismissButton(Meshgen_gui& gui) : gui_(gui) { }
    ~OnErrorCreateDismissButton() { }
  };
}

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

#ifdef NSF_DEMO
namespace QMG {
  namespace Nsf_Demo {
    using namespace QMG;
    using namespace QMG::MG;
    extern void indicate_phase_start(const string& phasename);
  }
}
#endif
          

// ------------------------------------------------------------------
// meshgen is the main mesh generation routine.


QMG::SimpComplex_Under_Construction
QMG::MG::meshgen(const Brep& g,
                 const string& default_size_control,
                 const string& default_curv_control,
                 const string& user_data,
                 const vector<Real>& expa,
                 const vector<Real>& asp_degrade,
                 ostream& logfile,
                 int verbosity1,
                 QMG::FrontEnd::Interpreter& interp,
                 bool gui_active,
                 double tol,
                 int ori) {




  using std::endl;

  // initialize the log

  Logstream logstr(logfile, verbosity1);

  // Initialize the object to cause errors to go to the log.

  OnErrorFlushLog unnamed1_(logstr);
  debug_ostream = &logstr.str();

  // initialize the random number generator.

  random_real(1);


  int di = g.embedded_dim();

  // Check the ranges of parameters.

  if (g.gdim() != di) 
    throw_error("Intrinsic and embedded dims of input brep must be equal.");

  if (di < MINDIM || di > MAXDIM) 
    throw_error("Dimension of brep is out of range.");

  if (g.level_size(di) == 0) 
    throw_error("Brep has no chambers.");

  if (expa.size() != di + 1)
    throw_error("Wrong size for expa");
  if (asp_degrade.size() != di + 1)
    throw_error("Wrong size for asp_degrade");
  {
    for (int i = 0; i < di + 1; ++i) {
      if (expa[i] < 0 || expa[i] > 1) 
        throw_error("expansion factor is out of range (must be in [0,1])");
      if (asp_degrade[i] < 0 || asp_degrade[i] > 1) 
        throw_error("asp_degrade factor is out of range (must be in [0,1])");
    }
  }

  if (ori != 1 && ori != 2)
    throw_error("invalid value for 'orientiation' argument");


  // Print the calling arguments in the log.

  if (logstr.verbosity() >= 1) {
    logstr.str() << "QMG 2.0 MESH GENERATION LOG (verbosity = " 
      << logstr.verbosity () << ") started on ";
    logstr.logdate();
    logstr.str() << endl << "Input brep:\n" << g << endl;
    logstr.str() << "Default size function = '" << default_size_control << "'" << endl;
    logstr.str() << "Default curvature function = '" << default_curv_control << "'" << endl;
    logstr.str() << "User data = " << user_data << endl;
    logstr.str() << "expa = ";
    {
      for (int j = 0; j < asp_degrade.size(); ++j)
        logstr.str() << expa[j] << " ";
    }
    logstr.str() << "asp degrade = ";
    {
      for (int j = 0; j < asp_degrade.size(); ++j)
        logstr.str() << asp_degrade[j] << " ";
    }
    logstr.str() << endl;
    logstr.str() << "verbosity = " << verbosity1 << endl;
    logstr.str() << "gui_active = " << gui_active << endl;
    logstr.str() << "tol = " << tol << endl;
  }



  //Initialize static data.

  small_initialize_static_data(di);
  boxstack_init_static_data(di); 
  alignInitStaticData(di);
  Point::InitStaticData(di);
#ifdef DEBUGGING
  Local_to_PatchMath_CPP::dump_everything = false;
  Local_to_PatchMath_CPP::debug_str = &logstr.str();
#endif


 /*
 
IntCoord coord;
coord[0] = 0x4000000U;
coord[1] = 0xc000000U;
int testiw = 4;
*/

  // Create the mesh generation gui panel.

  Meshgen_gui gui(gui_active, interp, g, logstr);

  // Set up a class to insert a button labeled "Dismiss"
  // if an error occurs.

  OnErrorCreateDismissButton unnamed2_(gui);


  gui.operation(string("Preprocessing Brep"));
  double scfac = g.scale_factor();

  vector<int> level_size0(di + 1);
  for (int fdim = 0; fdim <= di; ++fdim)
    level_size0[fdim] = g.level_size(fdim);

  SimpComplex_Under_Construction sc(di, di, level_size0);

  
  string brep_id = g.lookup_brep_propval(Brep::global_id_propname());

  if (brep_id.length() == 0) {
    if (logstr.verbosity() >= 1) {
      logstr.str() << "global id propval not found\n";
    }
  }
  else {
    if (logstr.verbosity() >= 1) {
      logstr.str() << "global id propval value = " << brep_id << '\n';
    }
    sc.add_prop_val(Brep::global_id_propname(), brep_id);
  }

  vector<Brep::Face_Spec> vtx_owner;

  BoxFaceDag dag(di);

  gui.associate_sc_and_dag(sc, dag);

 
  gui.operation("Quadtree generation");

  if (logstr.verbosity() >= 1) 
    logstr.str() << "\n\n-------------- Quadtree generation --------------------" << endl;
  
  // Open brace to make separation stacks and orbits vanish after this code block.

  {
  
    // Make the patch table.  Initialize the incidence table.
    PatchTable patchtable(g, logstr, scfac, tol);    
    IncidenceTable inc_table;


    // Gather size control info.
    
    SizeControl size_control(g, 
      patchtable,
      "sizecontrol",
      "gm_sizecontrol",
      default_size_control,
      user_data,
      interp,
      logstr);
    
    SizeControl curvature_control(g,
      patchtable,