gmcoarsetopo.cpp

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

// ------------------------------------------------------------------
// gmcoarsentopo.cpp
//
// This file contains a function for coarsening the topology of
// a brep.
// ------------------------------------------------------------------
// 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 "qbrep_constr.h"
#include "qpatchtab.h"
#include "qfrontend.h"
#include "qunionfind.h"
#include "qancestor.h"
#include "qsmall.h"

#define GM_ROUTINE_NAME gmcoarsetopo
#define GM_ROUTINE_NAME_Q "gmcoarsetopo"


namespace {
  using namespace QMG;


  // Duplicate of stuff from qseparation.cpp.  Fix later.

  const int NUM_TEST_DIRECTIONS_2D = 4;
  const int NUM_CURVATURE_TEST_DIRECTION = 13;

  const Real test_direction_init[] = {0,1,0, 1,0,0, 1,1,0, 1,-1,0,
                                      0,0,1,  0,1,1, 1,0,1, 1,1,1,
                                      0,-1,1, -1,0,1, 1,-1,1, -1,1,1,
                                      -1,-1,1};

 
  
  

  // Find the deviations along all edge points.  Similar to gmchecknormals.cpp
  
  Real maxip[NUM_CURVATURE_TEST_DIRECTION];
  Real minip[NUM_CURVATURE_TEST_DIRECTION];
  Point test_directions[NUM_CURVATURE_TEST_DIRECTION];
  int num_test_directions;

  void init_test_directions(int di) {
    num_test_directions = (di == 2)?    
      NUM_TEST_DIRECTIONS_2D : NUM_CURVATURE_TEST_DIRECTION;
    for (int j1 = 0; j1 < num_test_directions; ++j1) {
      for (int k = 0; k < di; ++k) 
        test_directions[j1][k] = test_direction_init[j1 * 3 + k];
      test_directions[j1].normalize();
    }
  }
  
  void init_maxmin_() {         
    for (int j = 0; j < num_test_directions; ++j) {
      maxip[j] = -100;
      minip[j] = 100;
    }
  }

              // Update maxip and minip with normal_or_tan.
              
  bool update_maxmin_(const Point& normal_or_tan, double allowable_curvature) {
    for (int j = 0; j < num_test_directions; ++j) {
      Real ip = Point::inner_product(test_directions[j], normal_or_tan);
      if (ip > maxip[j])
        maxip[j] = ip;
      if (ip < minip[j])
        minip[j] = ip;
      if (maxip[j] - minip[j] > allowable_curvature) {
        return false;
      }
    }
    return true;
  }

  bool twovec_ok_for_curve_bound(const Point& nrml1,
    const Point& nrml2,
    double allowable_curvature) {
    init_maxmin_();
    update_maxmin_(nrml1, allowable_curvature);
    return update_maxmin_(nrml2, allowable_curvature);
  }


  // Copied from qpatchtab.cpp.  Fix later.

  bool lookup_orientation_[33];
  void init_lookup_orientation() {
    lookup_orientation_[0] = true;
    lookup_orientation_[2] = true;
    lookup_orientation_[3] = false;
    lookup_orientation_[10] = false;
    lookup_orientation_[12] = false;
    lookup_orientation_[13] = true;
    lookup_orientation_[20] = false;
    lookup_orientation_[21] = true;
    lookup_orientation_[23] = true;
    lookup_orientation_[30] = true;
    lookup_orientation_[31] = false;
    lookup_orientation_[32] = false;
  }


  // New iterator: loop over descendents of a brep face, generally with
  // with repeats.

  class Fspec_loop_over_descendents : public Brep::Face_Spec {
  private:
    const Brep& b_;
    vector<Brep::Face_Spec> stack_;
    void set_this_();
    bool done_;
    // No copying, no assignment
    Fspec_loop_over_descendents(const Fspec_loop_over_descendents& o) :
    b_(o.b_) { }
    void operator=(const Fspec_loop_over_descendents&) { }
    Brep::Face_Spec& operator=(const Brep::Face_Spec& fspec) {return *this;}
  public:
    Fspec_loop_over_descendents(const Brep& b) :
        b_(b)  {set_this_();}
    void init(const Brep::Face_Spec& fspec);
    bool notdone();
    void operator++();
    ~Fspec_loop_over_descendents() {}
  };

  void Fspec_loop_over_descendents::set_this_() {
    if (stack_.size() == 0) {
      done_ = true;
      fdim_ = -1;
      faceind_ = -1;     
    } 
    else {
      fdim_ = stack_.back().fdim();
      faceind_ = stack_.back().faceind();
    }
  }


  void Fspec_loop_over_descendents::init(const Brep::Face_Spec& fspec) {
    stack_.resize(1,fspec);
    done_ = false;
    set_this_();
  }

  bool Fspec_loop_over_descendents::notdone() {
    set_this_();
    return !done_;
  }

  void Fspec_loop_over_descendents::operator++() {
    Brep::Face_Spec stacktop = stack_.back();
    stack_.pop_back();
    for (Brep::Face_Spec_Loop_Over_Face_Subfaces subfspec(b_, stacktop);
    subfspec.notdone();
    ++subfspec) {
      stack_.push_back(subfspec);
    }
  }
  

  Brep_Under_Construction coarsen_topo(const Brep& b, 
    double allowable_curvature) {


    // std::ofstream ctlog("ctlog");

    int di = b.embedded_dim();
    Point::InitStaticData(di);

    int gdim = b.gdim();
    if (gdim != di)
      throw_error("Brep must be full dimensional");
    ostringstream nullstr;

    using QMG::MG::Logstream;
    using QMG::MG::PatchTable;
    using QMG::MG::PatchTableIndexMap;
    using QMG::MG::PatchInBrepAddress;
    using QMG::MG::My_bool;

    Logstream null_log(nullstr,0);
    // Make a map to determine which parent-child patch relationships.
    PatchTable patchtable(b, null_log, -1, -1);
    PatchTableIndexMap<My_bool> already_seen(patchtable, My_bool(false));
    Brep::Ancestor_Lookup anc_lookup(b);
  
    init_test_directions(di);
    init_lookup_orientation();

    
    map<PatchInBrepAddress, PatchTable::Index> lookup_patchind; 
    {
      for (PatchTable::Index patchind = 0; 
      patchind < patchtable.numpatch(); 
      ++patchind) {
        Brep::Face_Spec owner = patchtable.owner(patchind);
        if (patchtable.gdim(patchind) != owner.fdim()) continue;
        Brep::PatchIndex patchind1 = patchtable.orig_index(patchind);
        lookup_patchind[PatchInBrepAddress(owner,patchind1)] = patchind;
      }
    }

    // Labeling to keep track of faces to delete.

    Brep::Face_Labeling<bool> ok_to_delete(b, true);

    // Union-find for merging faces.

    UnionFind uf[MAXDIM + 1];

    // Not ok to delete faces of dimension di or di - 1.

    {
      for (int fdim = di; fdim >= di - 1; --fdim) {
        uf[fdim].initialize(b.level_size(fdim));
        for (Brep::Face_Spec_Loop_Over_Faces_Of_Dim fspec(b, fdim);
        fspec.notdone();
        ++fspec) {
          ok_to_delete(fspec) = false;
        }
      }
    }

    vector<int> parent1_parentlabels;
    vector<int> parent2_parentlabels;

    // Loop over dimensions of faces that might be ok to delete.

    {
      for (int fdim = di - 2; fdim >= 0; --fdim) {
        uf[fdim].initialize(b.level_size(fdim));

        for (Brep::Face_Spec_Loop_Over_Faces_Of_Dim fspec(b, fdim);
        fspec.notdone();
        ++fspec) {

          // ctlog << " considering whether to delete " << fspec << " name " 
            // << b.face_name(fspec) << '\n';

          // Don't delete a face with prop-val pairs.

          int numpv = 0;
          {
            for (Brep::Loop_over_propvals_of_face pvloop1(b, fspec);
            pvloop1.notdone();
            ++pvloop1) {
              ++numpv;
            }
          }


          if (numpv > 0) { 
            ok_to_delete(fspec) = false;

            // ctlog << " not deleting; face has property-values\n";
            continue;
          }     
          
          // Count the number of parents

          int anc_count = 0;
          Brep::Face_Spec parent1, parent2;

          for (Brep::Ancestor_Lookup::Loop_over_ancestors ancloop(anc_lookup, fspec);
          ancloop.notdone();
          ++ancloop) {
            Brep::Face_Spec afspec = ancloop.ancestor_fspec();
            if (afspec.fdim() == fdim + 1 && !ok_to_delete(afspec)) {
              if (anc_count == 0)
                parent1 = afspec;
              if (ancloop.occurrence_count() > 1)
                ++anc_count;
              if (anc_count > 0)
                parent2 = afspec;
              ++anc_count;
            }
          }

          // If the number of parents is not 0 and not 2, can't delete,

          if (anc_count != 0 && anc_count != 2) {
            ok_to_delete(fspec) = false;
            // ctlog << " not deleting; # ancestors = " << anc_count << "\n";
            continue;
          }
          if (anc_count == 2) {

            // ctlog << " two parents " << parent1 << " (" << b.face_name(parent1)
             // << ") and " << parent2 << " (" << b.face_name(parent2) << ")\n";

            // Check that the two parents have the same p-v-pairs.
            
            for (Brep::Loop_over_propvals_of_face pvloop(b, parent1);
            pvloop.notdone();
            ++pvloop) {