gmchecknormals.cpp

来自「算断裂的」· C++ 代码 · 共 356 行

// ------------------------------------------------------------------
// gmchecknormals.cpp
//
// This file contains a function for checking how smooth the normals
// are across patch boundaries (departure from G^1).
// ------------------------------------------------------------------
// 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.h"
#include "qpatchtab.h"
#include "qfrontend.h"

#define GM_ROUTINE_NAME gmchecknormals
#define GM_ROUTINE_NAME_Q "gmchecknormals"


namespace {
  using namespace QMG;

  // Duplication of class defined in double.cpp.  Fix later.

  
  struct PatchInBrepAddress {
    Brep::Face_Spec fspec;
    Brep::PatchIndex patchind;
    PatchInBrepAddress(const Brep::Face_Spec& f, Brep::PatchIndex p) : 
    fspec(f), patchind(p) { }
    PatchInBrepAddress(const PatchInBrepAddress& o) : fspec(o.fspec), patchind(o.patchind) { }
    PatchInBrepAddress() { }
    bool operator<(const PatchInBrepAddress& ip) const;
    bool operator==(const PatchInBrepAddress& ip) const;
  };

  
  
  bool PatchInBrepAddress::operator<(const PatchInBrepAddress& ip) const {
    if (fspec < ip.fspec) return true;
    if (fspec > ip.fspec) return false;
    return patchind < ip.patchind;
  }
  
  bool PatchInBrepAddress::operator==(const PatchInBrepAddress& ip) const {
    return fspec == ip.fspec && patchind == ip.patchind;
  }



  // 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};
  







  struct CpnumFaceSpec {
    Brep::ControlPointIndex cpnum;
    Brep::FaceIndex faceind;
    CpnumFaceSpec() { }
    CpnumFaceSpec(Brep::ControlPointIndex cp1, Brep::FaceIndex f) :
     cpnum(cp1), faceind(f) { }
    CpnumFaceSpec(const CpnumFaceSpec& o) :
     cpnum(o.cpnum), faceind(o.faceind) { }
    bool operator<(const CpnumFaceSpec& o) const;
  };

  bool CpnumFaceSpec::operator<(const CpnumFaceSpec& o) const {
    if (cpnum < o.cpnum) return true;
    if (cpnum > o.cpnum) return false;
    return faceind < o.faceind;
  }


    
  struct FacePatchInd {
    Brep::PatchIndex patchind;
    int whichvert;
    FacePatchInd() { }
    FacePatchInd(Brep::PatchIndex p, int w) :
     patchind(p), whichvert(w) { }
    FacePatchInd(const FacePatchInd& o) :
     patchind(o.patchind), whichvert(o.whichvert) { }
  };

  typedef multimap<CpnumFaceSpec, FacePatchInd> Maptype;

  struct Check_normals_return_val {
    Real max_deviation_on_face;
    Brep::Face_Spec face_with_max_deviation;
    Brep::ControlPointIndex node_with_max_deviation_on_face;

    /*
    Real max_deviation_overall;
    Brep::ControlPointIndex node_with_max_deviation_overall;
    */
  };



  Check_normals_return_val check_normals(const Brep& b) {

    Maptype cpmap;
    map<PatchInBrepAddress, bool> patch_is_flipped;
    int di = b.embedded_dim();
    Point::InitStaticData(di);

    int gdim = b.gdim();
    if (gdim < di - 1)
      throw_error("Brep must be full dimensional");
    {  

      ostringstream nullstr;
      QMG::MG::Logstream null_log(nullstr,0);
      // Make a map to determine which patch is flipped.
      QMG::MG::PatchTable patchtable(b, null_log, -1, -1);
      for (QMG::MG::PatchTable::Index patchind = 0; 
      patchind < patchtable.numpatch();
      ++patchind) {
        if (patchtable.gdim(patchind) != di - 1) continue;
        patch_is_flipped[PatchInBrepAddress(patchtable.owner(patchind), 
          patchtable.orig_index(patchind))] = patchtable.is_forward(patchind);
      }
    }

    for (Brep::Face_Spec_Loop_Over_Faces_Of_Dim fspec(b, di - 1);
    fspec.notdone();
    ++fspec) {  
      Brep::FaceIndex faceind = fspec.faceind();
      for (Brep::Loop_over_patches_of_face patchloop(b, fspec);
      patchloop.notdone();
      ++patchloop) {
        Brep::PatchIndex patchind = patchloop.index();
        if (di == 2) {
          int deg = patchloop.degree1();
          cpmap.insert(pair<CpnumFaceSpec, FacePatchInd>
            (CpnumFaceSpec(patchloop.control_point(0), faceind),
            FacePatchInd(patchind, 0)));
          cpmap.insert(pair<CpnumFaceSpec, FacePatchInd>
            (CpnumFaceSpec(patchloop.control_point(deg), faceind),
            FacePatchInd(patchind, 1)));
        }
        else { // di == 3
          if (patchloop.ptype() == BEZIER_TRIANGLE) {
            int deg = patchloop.degree1();
            cpmap.insert(pair<CpnumFaceSpec, FacePatchInd>
              (CpnumFaceSpec(patchloop.control_point(0), faceind),
              FacePatchInd(patchind, 0)));
            cpmap.insert(pair<CpnumFaceSpec, FacePatchInd>
              (CpnumFaceSpec(patchloop.control_point(deg * (deg + 1) / 2), faceind), 
              FacePatchInd(patchind, 1)));
            cpmap.insert(pair<CpnumFaceSpec, FacePatchInd>
              (CpnumFaceSpec(patchloop.control_point((deg + 1) * (deg + 2) / 2 - 1), faceind), 
              FacePatchInd(patchind, 2)));
          }
          else { // patchtype == BEZIER_QUAD
            int deg1 = patchloop.degree1();
            int deg2 = patchloop.degree2();
            cpmap.insert(pair<CpnumFaceSpec, FacePatchInd>
              (CpnumFaceSpec(patchloop.control_point(0), faceind),
              FacePatchInd(patchind, 0)));
            cpmap.insert(pair<CpnumFaceSpec, FacePatchInd>
              (CpnumFaceSpec(patchloop.control_point(deg1), faceind),  
              FacePatchInd(patchind, 1)));
            cpmap.insert(pair<CpnumFaceSpec, FacePatchInd>
              (CpnumFaceSpec(patchloop.control_point(deg2 * (deg1 + 1)), faceind),
              FacePatchInd(patchind, 2)));
            cpmap.insert(pair<CpnumFaceSpec, FacePatchInd>
              (CpnumFaceSpec(patchloop.control_point((deg2 + 1) * (deg1 + 1) - 1), faceind),
              FacePatchInd(patchind, 3)));
          }
        }
      }
    }

    // Now find maximum deviations.

    Check_normals_return_val rval;
    rval.max_deviation_on_face = -1.0;

    //    rval.max_deviation_overall = -1.0;

    Real maxip[NUM_CURVATURE_TEST_DIRECTION];
    Real minip[NUM_CURVATURE_TEST_DIRECTION];
    int num_test_directions = (di == 2)?    
      NUM_TEST_DIRECTIONS_2D : NUM_CURVATURE_TEST_DIRECTION;

    vector<Point> test_directions(num_test_directions);
    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();
    }
    
    Point paramcoord, normal_or_tan, direc;
    direc[0] = 1.0;


    // Pass 2 looks for maxdeviation per node/fspec pair.
    
    
    
    // for (int passcount = 0;
    // passcount < 2;
    // ++passcount) {
    
    Brep::ControlPointIndex lastv = -1;
    Brep::FaceIndex lastface = -1;
    Maptype::const_iterator start_search = cpmap.begin();
    while (start_search != cpmap.end()) {
      Brep::ControlPointIndex thisv = start_search -> first.cpnum;
      Brep::FaceIndex thisf = start_search -> first.faceind;
      Maptype::const_iterator end_search =
        
        // (passcount == 0) ?
        // cpmap.lower_bound(CpnumFaceSpec(thisv + 1, 0)) :
        
        cpmap.upper_bound(CpnumFaceSpec(thisv, thisf));        
      for (int j = 0; j < num_test_directions; ++j) {
        maxip[j] = -100;
        minip[j] = 100;
      }
      for (Maptype::const_iterator it = start_search; 
      it != end_search;
      ++it) {
        Brep::FaceIndex thisf2 = it -> first.faceind;
        Brep::Face_Spec owner(di - 1, thisf2);
        Brep::PatchIndex patchind = it -> second.patchind;
        int whichvert = it -> second.whichvert;
        if (di == 2) {
          paramcoord[0] = static_cast<double>(whichvert);
          normal_or_tan = b.patchmath(owner, patchind).direc_deriv(paramcoord, direc);
          if (normal_or_tan.l2norm() == 0)
            throw_error("Tangential derivative vanishes at edge endpoint -- illegal Bezier curve");
          normal_or_tan.normalize();
        }
        else { //di == 3
          if (b.patch_type(owner, patchind) == BEZIER_TRIANGLE) {
            if (whichvert == 0) {
              paramcoord[0] = 0.0;
              paramcoord[1] = 1.0;
            } 
            else if (whichvert == 1) {
              paramcoord[0] = 0.0;
              paramcoord[1] = 0.0;
            }
            else { //whichvert == 2
              paramcoord[0] = 1.0;
              paramcoord[1] = 0.0;
            }
          }
          else { // patchtype == BEZIER_QUAD
            if (whichvert == 0) {
              paramcoord[0] = 0.0;
              paramcoord[1] = 0.0;
            }
            else if (whichvert == 1) {
              paramcoord[0] = 1.0;
              paramcoord[1] = 0.0;
            }
            else if (whichvert == 2) {
              paramcoord[0] = 0.0;
              paramcoord[1] = 1.0;
            }
            else { //whichvert == 3
              paramcoord[0] = 1.0;
              paramcoord[1] = 1.0;
            }
          }
          normal_or_tan = b.patchmath(owner,patchind).normal(paramcoord);
        }
        
        if (patch_is_flipped[PatchInBrepAddress(owner, patchind)]) {
          normal_or_tan[0] = -normal_or_tan[0];
          normal_or_tan[1] = -normal_or_tan[1];
          if (di == 3)
            normal_or_tan[2] = -normal_or_tan[2];
        }
        
        // Update maxip and minip with normal_or_tan.
        
        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;
        }
      }
      
      // Done with all relevant directions; now find the max gap in the
      // inner products.
      
      Real maxgap = 0.0;
      for (int j2 = 0; j2 < num_test_directions; ++j2) {
        Real gap = maxip[j2] - minip[j2];
        if (gap > maxgap) {
          maxgap = gap;
        }
      }
      
      // if (passcount == 0) {
      //  if (maxgap > rval.max_deviation_overall) {
      //    rval.max_deviation_overall = maxgap;
      //     rval.node_with_max_deviation_overall = thisv;
      //  }
      //  }
      //  else { //passcount == 1
      if (maxgap > rval.max_deviation_on_face) {
        rval.max_deviation_on_face = maxgap;
        rval.face_with_max_deviation = Brep::Face_Spec(di - 1, thisf);
        rval.node_with_max_deviation_on_face = thisv;
      }
      //  }
      start_search = end_search;
      }
      // }
      return rval;
  }

  void worker(const QMG::FrontEnd::ArgValType& argvalhandle,
    QMG::FrontEnd::ReturnValType& returnvalhandle,
    QMG::FrontEnd::Interpreter& interp) {

    argvalhandle.verify_nargin(1, 1, GM_ROUTINE_NAME_Q);
    returnvalhandle.verify_nargout(3, 3, GM_ROUTINE_NAME_Q);
    Brep_From_FrontEnd b = argvalhandle.get_brep(0);
    Check_normals_return_val rval = check_normals(b);
    returnvalhandle.put_double(0, rval.max_deviation_on_face);
    returnvalhandle.put_string(1, b.face_name(rval.face_with_max_deviation));
    returnvalhandle.put_int(2, rval.node_with_max_deviation_on_face);
    // returnvalhandle.put_double(3, rval.max_deviation_overall);
    // returnvalhandle.put_int(4, rval.node_with_max_deviation_overall);
  }
}

#include "gm_entrypoint.cpp"