checktri.cpp

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

// ------------------------------------------------------------------
// checktri.cpp
//
// This file contains the routine for checking that a triangulation
// of a brep is correct.
// ------------------------------------------------------------------
// 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.
// ------------------------------------------------------------------

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

#include "qbrep.h"
#include "qsimpcomp.h"
#include "qmatvec.h"


// Function exported by this file.

namespace QMG {
  Real checktri(const Brep& brep,
    const SimpComplex& simpcomp,
    ostream& os,
    int checko,
    double tol);
}



namespace {
  using namespace QMG;

  const int CHECKTRI_MAXDIM = 3;
  
  // ------------------------------------------------------------------
  // struct Facet_Rec
  // A record of a facet of a simplex.  Used as a map data.  Keep track
  // of how many times a facet occurs.

  struct Facet_Rec {
    int pos_occurrence_count;
    int neg_occurrence_count;
    int occurs_on_bdry;
    int seqno_occ[2];
    Facet_Rec() : pos_occurrence_count(0),
      neg_occurrence_count(0),
      occurs_on_bdry(0) { }
  };
  
  // ------------------------------------------------------------------
  // struct Facet_Key
  // A record of a facet of a simplex.  Used as a map key. The key is
  // the list of vertices of the facet, sorted into increasing order.


  struct Facet_Key {
    static int facetdim;
    SimpComplex::VertexGlobalIndex vert[CHECKTRI_MAXDIM];
    bool operator<(const Facet_Key& other) const;
    bool operator==(const Facet_Key& other) const;
  };

  bool Facet_Key::operator<(const Facet_Key& other) const {
    for (int j = 0; j < facetdim + 1; ++j) {
      if (vert[j] < other.vert[j])
        return true;
      if (vert[j] > other.vert[j])
        return false;
    }
    return false;
  }

  int Facet_Key::facetdim;

  bool Facet_Key::operator==(const Facet_Key& other) const {
    for (int j = 0; j < facetdim + 1; ++j) {
      if (vert[j] != other.vert[j])
        return false;
    }
    return true;
  }

  class Dumpsimp {
  private:
    const SimpComplex& s;
    const Brep& g;
    Brep::Face_Spec fspec;
    int seqno;
  public:
    Dumpsimp(const SimpComplex& s1, 
      const Brep& g1,       
      const Brep::Face_Spec& f1,
      int seq1) :
    s(s1), g(g1), fspec(f1), seqno(seq1) {}
    friend QMG::ostream& operator<<(QMG::ostream& os, const Dumpsimp& ds);
  };

  QMG::ostream& operator<<(QMG::ostream& os, const Dumpsimp& ds) {
    os << "simplex #" << ds.seqno << " of topological entity "
      << ds.g.face_name(ds.fspec) << " (";
    ::operator<<(os, ds.fspec); 
    os << ") which has vertices ";
    for (int j = 0; j < ds.fspec.fdim() + 1; ++j)
      os << ds.s.node_of_meshface_on_face(ds.fspec, ds.seqno, j) << " ";
    os << '\n';
    return os;
  }
}

  
// ------------------------------------------------------------------
// checktri
// Argument g is a brep and argument s is a putative triangulation
// for g.  This routine checks whether the simplicial complex is
// correct.  Problems detected are printed on outstream.  
// The return value is the worst aspect ratio among all simplices.
// It is negative if there is an error.
// Set checko=1 if every elt should have forward orientation,
// 2 if every elt should have backward orientation, 0 not to check orientation.
// Note: orientation not checked for low-dimensional complexes.
// ------------------------------------------------------------------

QMG::Real 
QMG::checktri(const Brep& g,
              const SimpComplex& s,
              ostream& outstream,
              int checko,
              double tol) {
  
  int di = g.embedded_dim();
  if (di != s.embedded_dim()) {
    throw_error("Brep and simpcomplex have different embedded dims");
  }
  
  // initialize the random number generator.

  random_real(1);
  Point::InitStaticData(di);

  string brep_id1 = g.lookup_brep_propval(Brep::global_id_propname());
  string brep_id2 = s.lookup_prop_val(Brep::global_id_propname());
  bool match = compare_nocase(brep_id1, brep_id2) == 0;
  if (!match)
    throw_error("Brep and simpcomplex have different global id's");
  int sdim = s.gdim();

  if (sdim > CHECKTRI_MAXDIM) {
    throw_error("Dimension out of range for checktri");
  }
  if (sdim > g.gdim()) {
    throw_error("Intrinsic dimension of simp complex too high for brep");
  }
  double scfac = g.scale_factor();
  double scaled_tol = scfac * tol;
  
  Real minaltglobal = BIG_REAL;
  Real maxsideglobal = 0.0;
  Real maxaspglobal = 0.0;
  int worsttriseq = -999;
  Brep::Face_Spec worsttrifspec;
  bool imperfectflag = false;
  
  Matrix simplexmatrix(sdim, di);
  Matrix simplexvertex(sdim + 1, di);

  // hash table on facets; counts occurrences of facets on the
  // positive and negative sides.

  typedef map<Facet_Key, Facet_Rec> Facet_Map_Type;
  Facet_Map_Type facet_occurrence_map;
  map<SimpComplex::VertexGlobalIndex, int> checkoff;
  map<SimpComplex::VertexGlobalIndex, int> onface;
  int numvtx = s.num_nodes();
  int numelt1 = 0;
  vector<int> sortord(sdim + 1);
  Real hhtrans[4];

  for (int fdim = 0; fdim <= sdim; ++fdim) {
    if (g.level_size(fdim) != s.brep_levelsize(fdim))
      throw_error("Mismatch between brep and mesh (different level sizes)");
    Facet_Key::facetdim = fdim - 1;
    sortord.resize(fdim + 1);
    for (Brep::FaceIndex faceind = 0; faceind < g.level_size(fdim); ++faceind) {
      Brep::Face_Spec fspec(fdim, faceind);
      if (fdim < di) {
        int nnf = s.num_node_on_face(fspec);
        if (nnf == 0) {
          outstream << "?? Topological face " << g.face_name(fspec) << " (" << fspec 
            << ") does not have any mesh nodes\n";
          imperfectflag = true;
        }
        onface.clear();
      
        for (int seqno = 0; seqno < nnf; ++seqno) {
          SimpComplex::VertexGlobalIndex vnum = s.node_on_face(fspec,seqno);
          onface[vnum] = 1;
          SimpComplex::VertexOrdinalIndex vnumo = s.global_to_ordinal_nothrow(vnum);
          if (vnumo < 0) {
            outstream << "?? Global node #" << vnum << " listed as lying on face "
              << g.face_name(fspec) << "(" << fspec << ") but no node\n"
              << "  with that number occurs in the coordinate list.\n";
            imperfectflag = true;
            continue;
          }
          Point real1, param1;
          for (int j = 0; j < di; ++j)
            real1[j] = s.real_coord_o(vnumo, j);
          for (int j2 = 0; j2 < fdim; ++j2)
            param1[j2] = s.param_coord_on_face(fspec, seqno,j2);
          Brep::PatchIndex patchind = s.patchind_on_face(fspec, seqno);

          Point real2 = g.patchmath(fspec, patchind).get_real_point(param1);
          Real l2dist = Point::l2dist(real1,real2);
          if (l2dist > scfac * tol) {
            outstream << "?? Vertex " << vnum << " real coordinates too far away "
              << " from actual point on " << g.face_name(fspec) << " (" 
              << fspec << ")\n"
              << "  ( dist = " << l2dist << ")\n";
            imperfectflag = true;
          }
        }
      }
      if (fdim == 0) continue;

      int nmf = s.num_meshface_on_face(fspec);
      facet_occurrence_map.clear();
     
      {
        Facet_Key key;  
        for (int seqno = 0; seqno < nmf; ++seqno) {
          if (fdim < di) {
            for (int jj = 0; jj < fdim + 1; ++jj) {
              SimpComplex::VertexGlobalIndex vtx = 
                s.node_of_meshface_on_face(fspec, seqno,jj);
              if (!onface[vtx]) {
                outstream << "?? Vertex " << vtx << " occurs in " 
                  << Dumpsimp(s,g,fspec,seqno) << " but is\n"
                  << " not listed as a vertex on that topological entity.\n";
                imperfectflag = true;
              }
            }
          }
          bool notfound = false;
          for (int jj = 0; jj < fdim + 1; ++jj) {
            SimpComplex::VertexGlobalIndex vtx = 
              s.node_of_meshface_on_face(fspec, seqno, jj);
            SimpComplex::VertexOrdinalIndex vtxo = 
              s.global_to_ordinal_nothrow(vtx);
            if (vtxo < 0) {
              outstream << "?? Vertex " << vtx << " occurs in " 
                << Dumpsimp(s,g,fspec,seqno) << " but has no coordinates in table.\n";
              imperfectflag = true;
              notfound = true;
            }
            else {
              for (int j = 0; j < di; ++j)
                simplexvertex(jj,j) = s.real_coord_o(vtxo,j);
            }