meshgen.cpp

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

// ------------------------------------------------------------------
// meshgen.cpp
//
// This file contains the main routine for mesh generation.
// ------------------------------------------------------------------
// 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.
// ------------------------------------------------------------------

#include <iomanip>

#ifdef _MSC_VER
#if _MSC_VER == 1200
#include "qmatvec.h"
#endif
#endif

#include "qboxstack.h"
#include "qlog.h"
#include "qsizectl.h"
#include "qerr.h"
#include "qseparation.h"
#include "qmg_gui.h"
#include "qlifter.h"


// declare extern functions.


QMG::ostream* debug_ostream;

namespace QMG {
  namespace MG {

    using namespace QMG;


    // This is the function exported from this file.

    SimpComplex_Under_Construction 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);

    // See align.cpp for these.

    extern void alignInitStaticData(int di); 


    
    extern void 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);


    extern void 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);


    // see qboxstack.cpp
    extern void boxstack_init_static_data(int di); 

    // see qsmall.cpp
    extern void small_initialize_static_data(int di);
  }
}


// Variables and defs with file scope.


namespace {
  using namespace QMG;
  using namespace QMG::MG;


  // ------------------------------------------------------------------
  // Boxvecs_by_level
  // An array of ActiveBoxVec's, indexed by level.

  class Boxvecs_by_level {
  private:
    vector<ActiveBoxVec*> vecs_;
    Real expa_;
    void enlarge_(int lev) {
      if (lev >= vecs_.size()) {
        vecs_.resize(lev + 1, 0);
      }
    }
    Boxvecs_by_level(const Boxvecs_by_level&) { }
    void operator=(const Boxvecs_by_level&) { }
  public:
    ActiveBoxVec* operator[](int lev) {
      enlarge_(lev);
      return vecs_[lev];
    }
    ActiveBoxVec& force_creation(int lev) {
      enlarge_(lev);
      if (!vecs_[lev]) {
        vecs_[lev] = new ActiveBoxVec(expa_, lev);
      }
      return *(vecs_[lev]);
    }
    void reset_expa(Real expa) {
      for (int i = 0; i < vecs_.size(); ++i) {
        if (vecs_[i])
          vecs_[i] -> reset_expa(expa);
      }
      expa_ = expa;
    }
    void erase(int lev) {
      enlarge_(lev);
      if (vecs_[lev])
        delete vecs_[lev];
      vecs_[lev] = 0;
    }
    int size() const {return vecs_.size();}
    Boxvecs_by_level(Real expa) : expa_(expa) { }
    ~Boxvecs_by_level() {
      for (int lev = vecs_.size() - 1; lev >= 0; --lev) {
        if (vecs_[lev]) delete vecs_[lev];
      }
    }

    void dump_size(ostream& os) {
      for (int lev = 0; lev < vecs_.size(); ++lev) {
        os << "Level " << lev;
        if (vecs_[lev]) {
          os << " size " << vecs_[lev] -> numboxes() << '\n';
        }
        else {
          os << " not created\n";
        }
      }
    }
  };


  // ------------------------------------------------------------------
  // This local routine adds a simplex to the simplicial complex.
  // It computes the face owners and lifts parameters


  void add_simplex(SimpComplex_Under_Construction& sc,
    const vector<SimpComplex::VertexOrdinalIndex>& vertices,
    int surf_ori,
    int ori,
    const vector<Brep::Face_Spec>& vtx_owner,
    const Brep_Parameter_Lifter& lifter,
    // workspaces:
    vector<SimpComplex::VertexGlobalIndex>& vlist1,
    vector<Real>& param1,
    vector<Point>& params,
    vector<Brep::PatchIndex>& patchinds
    ) {
    int di = sc.embedded_dim();

    SimpComplex::VertexOrdinalIndex vnumf_o = vertices[di];
    Brep::Face_Spec lastowner = vtx_owner[vnumf_o];
    SimpComplex::VertexGlobalIndex vnumf = sc.ordinal_to_global(vnumf_o);
    if (lastowner.fdim() < di) {
      patchinds[di] = sc.retrieve_vertex_bdryinc_patchind(vnumf, lastowner);
      params[di] = sc.retrieve_vertex_bdryinc_param(vnumf,lastowner);
    }
    vlist1.resize(1);
    vlist1[0] = vnumf;
    for (int j1 = di - 1; j1 >= 0; --j1) {
      SimpComplex::VertexOrdinalIndex vnum_o = vertices[j1];
      SimpComplex::VertexGlobalIndex vnum = sc.ordinal_to_global(vnum_o);
      Brep::Face_Spec newowner = vtx_owner[vnum_o];
      int newfdim = newowner.fdim();
#ifdef DEBUGGING
      if (newfdim < lastowner.fdim() ||
          (newfdim == lastowner.fdim() && newowner != lastowner)) {
        throw_error("problem with vertex owner sequence");
      }
#endif
      if (newfdim < di) {
        patchinds[j1] = sc.retrieve_vertex_bdryinc_patchind(vnum, newowner);
        params[j1] = sc.retrieve_vertex_bdryinc_param(vnum,newowner);
        if (newowner != lastowner) {
          for (int j2 = j1 + 1; j2 < di + 1; ++j2) {
            SimpComplex::VertexOrdinalIndex vnum_o2 = vertices[j2];
            SimpComplex::VertexGlobalIndex vnum2 = sc.ordinal_to_global(vnum_o2); 
            pair<Point, Brep::PatchIndex> rval = lifter.lift_parameters(lastowner,
              patchinds[j2], params[j2], newowner);
            param1.resize(newfdim);
            for (int j3 = 0; j3 < newfdim; ++j3)
              param1[j3] = rval.first[j3];
            sc.add_vertex_bdryinc(vnum2, newowner, rval.second, param1);
            params[j2] = rval.first;
            patchinds[j2] = rval.second;
          }
        }
      }
      lastowner = newowner;
      vlist1.push_back(vnum);
      if (newfdim == di - j1) {
        bool swap1 = false;
        if (newfdim == di && ori) 
          swap1 = true;
        if (newfdim == di - 1) {
          if (surf_ori > 1) {
            throw_error("Missing surf ori");
          }
          if (di % 2 == 0)
            surf_ori = 1 - surf_ori;
          if ((surf_ori && ori) || (!surf_ori && !ori) )
            swap1 = true;
        }
        if (swap1) {
          SimpComplex::VertexGlobalIndex tmp = vlist1[0];
          vlist1[0] = vlist1[1];
          vlist1[1] = tmp;
        }
        sc.add_simplex_face(vlist1, newowner);
        if (swap1) {
          SimpComplex::VertexGlobalIndex tmp = vlist1[0];
          vlist1[0] = vlist1[1];
          vlist1[1] = tmp;
        }
      }
    }
  }

  
  // ------------------------------------------------------------------
  // The next few functions were used for debugging.  They compute
  // the total weight in all ActiveBoxVec's of a certain point.

  pair<bool, Base3>
    coord_is_inside_box(const ActiveBox& b,
                        const IntCoord& coord) {
    bool is_inside = true;
    Base3 c_rindex = 0;
    IntCoord lb = b.boxaddress().lowerleft();
    IntCoord ub = lb;
    int di = b.embedded_dim();
    int width = b.iwidth();

    for (int jj = 0; jj < di; ++jj) {
      if (!b.flatdim(jj))
        ub.increment(jj, width);
    }
    int reldim = 0;
    for (int j = 0; j < di; ++j) {
      if (coord[j] < lb[j]) {
        is_inside = false;
        break;
      }
      else if (coord[j] > ub[j]) {
        is_inside = false;
        break;
      }
      else if (!b.flatdim(j)) {
        int newdig;
        if (coord[j] == lb[j])
          newdig = 0;
        else if (coord[j] == ub[j])
          newdig = 1;
        else
          newdig = 2;
        c_rindex.set_digit(reldim++,newdig);
      }
    }
    return pair<bool, Base3>(is_inside, c_rindex);
  }
  
  
  Weight_t sum_weight_stack(const ActiveBoxVec& vec,
                            const IntCoord& coord, 
                            int iw, 
                            ActiveBoxVec::Create_Subbox_Workspace& wkspa) {
    Weight_t sum = 0;
    int di = vec.embedded_dim();
    int width = vec.iwidth();
    BoxAddress bxa;
    bxa.lowerleft()[0] = coord[0];
    bxa.lowerleft()[1] = coord[1];
    if (di == 3) {
      bxa.lowerleft()[2] = coord[2];
      bxa.flatdim() = 7;
    } 
    else {
      bxa.flatdim() = 3;
    }
    for (ActiveBoxVec::Loop_over_boxes loop(vec); loop.notdone(); ++loop) {
      ActiveBox b = loop.getbox();
      pair<bool, Base3> rval = 
        coord_is_inside_box(b, coord);
      bool is_inside = rval.first;
      if (is_inside) {
        Base3 c_rindex = rval.second;
        *debug_ostream << "***contained by ";
        //      b.full_dump2(*debug_ostream, patchtable, inc_table);
        if (width > 6)