qsimpcomp.cpp

算断裂的

// ------------------------------------------------------------------
// qsimpcomp.cpp
//
// This file contains member functions the class SimpComp and SimpComp_
// Under_Construction, which are for simplicial complexes (meshes).
// This version of qsimpcomp.cpp is for the Matlab front end.
// ------------------------------------------------------------------
// 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 "qsimpcomp.h"
#include "string.h"


namespace {
  void pause(const char* s) {
    mexPrintf("%s\nPress any key to continue\n", s);
    mexEvalString("pause");
  }
  void dump610(const mxArray* a) {
    mexPrintf("Dumping m{6}{1,0}\n");
    mxArray* a0 = mxGetCell(a,6);
    mxArray* a1 = mxGetCell(a0,1);
    int m = mxGetM(a1);
    int n = mxGetN(a1);
    mexPrintf("size = %d-by-%d\n", m,n);
    for (int i = 0; i < m; ++i) {
      for (int j = 0; j < n; ++j)
        mexPrintf("%10.5f ", mxGetPr(a1)[i + j * m]);
      mexPrintf("\n");
    }
  }
}



QMG::SimpComplex::VertexOrdinalIndex 
QMG::SimpComplex::process_ordinal_vnum_(VertexOrdinalIndex vnumo) const {
  if (vnumo < 0 || vnumo >= mxGetN(mxGetCell(c_entries_, 4))) {
    mexPrintf("%d %d\n", vnumo, mxGetN(mxGetCell(c_entries_, 4)));
    throw_error("ordinal vnum out of range");
  }
  return vnumo;
}

QMG::SimpComplex::SimpComplex(const SimpComplex& sc) :
  mesh_(sc.mesh_),
  c_embedded_dim_(sc.c_embedded_dim_),
  c_gdim_(sc.c_gdim_),
  c_entries_(sc.c_entries_),
  c_nodearray_(sc.c_nodearray_),
  global_to_ordinal_(sc.global_to_ordinal_)
{
 }


int 
QMG::SimpComplex::num_nodes() const {
  return mxGetN(mxGetCell(c_entries_, 4));
}


int
QMG::SimpComplex::num_elements() const {
  const mxArray* toplev = mxGetCell(c_entries_, c_gdim_ + 5);
  int numface = mxGetN(toplev);
  int nume = 0;
  for (Brep::FaceIndex faceind = 0; faceind < numface; ++faceind)
    nume += mxGetN(mxGetCell(toplev, 2 * faceind + 1));
  return nume;
}

int
QMG::SimpComplex::brep_levelsize(int fdim) const {
#ifdef RANGECHECK
  if (fdim < 0 || fdim > c_gdim_)
    throw_error("fdim out of range");
#endif
  return mxGetN(mxGetCell(c_entries_, fdim + 5));
}

QMG::string
QMG::SimpComplex::lookup_prop_val(const string& prop) const {
  const mxArray* pvlist = mxGetCell(c_entries_, 3);
  int npv = mxGetN(pvlist);
  for (int k = 0; k < npv; ++k) {
    const mxArray* propa = mxGetCell(pvlist, 2 * k);
    if (propa == 0 || !mxIsChar(propa) || mxGetM(propa) != 1
      || mxIsSparse(propa))
      throw_error("Illegal data type stored in property slot");
    int sz = mxGetN(propa);
    if (sz != prop.length()) continue;
    char* buf = new char[sz + 2];
    mxGetString(propa, buf, sz + 1);
    string props(buf);
    delete[] buf;
    if (compare_nocase(props, prop) == 0) {
      const mxArray* vala = mxGetCell(pvlist, 2 * k + 1);
      if (vala == 0 || !mxIsChar(vala) || mxGetM(vala) != 1
          || mxIsSparse(vala))
        throw_error("Illegal data type stored in value slot");
      sz = mxGetN(vala);
      buf = new char[sz + 2];
      mxGetString(vala, buf, sz + 1);
      string rval(buf);
      delete[] buf;
      return rval;
    }
  }
  return string("");
}

int
QMG::SimpComplex::num_prop_val() const {
  return mxGetN(mxGetCell(c_entries_, 3));
}

QMG::string 
QMG::SimpComplex::prop(int i) const {
  const mxArray* pvlist = mxGetCell(c_entries_, 3);
  int npv = mxGetN(pvlist);
#ifdef RANGECHECK
  if (i < 0 || i >= npv)
    throw_error("Index out of range in prop");
#endif
  const mxArray* propa = mxGetCell(pvlist, 2 * i);
  if (propa == 0 || !mxIsChar(propa) || mxGetM(propa) != 1
      || mxIsSparse(propa))
    throw_error("Illegal data type stored in property slot");
  int sz = mxGetN(propa);
  char* buf = new char[sz + 2];
  mxGetString(propa, buf, sz + 1);
  string prop_s(buf);
  delete[] buf;
  return prop_s;
}

QMG::string 
QMG::SimpComplex::val(int i) const {
  const mxArray* pvlist = mxGetCell(c_entries_, 3);
  int npv = mxGetN(pvlist);
#ifdef RANGECHECK
  if (i < 0 || i >= npv)
    throw_error("Index out of range in val");
#endif
  const mxArray* vala = mxGetCell(pvlist, 2 * i + 1);
  if (vala == 0 || !mxIsChar(vala) || mxGetM(vala) != 1
      || mxIsSparse(vala))
    throw_error("Illegal data type stored in value slot");
  int sz = mxGetN(vala);
  char* buf = new char[sz + 2];
  mxGetString(vala, buf, sz + 1);
  string val_s(buf);
  delete[] buf;
  return val_s;
}

QMG::SimpComplex::VertexGlobalIndex 
QMG::SimpComplex::ordinal_to_global(VertexOrdinalIndex vnumo) const {
  // mexPrintf("Reached ordinal to global vnumo = %d\n", vnumo);
  double vnumr = c_nodearray_[process_ordinal_vnum_(vnumo) *
                             (c_embedded_dim_ + 1)];
  VertexGlobalIndex vnum = static_cast<int>(vnumr);
  if (vnumr != static_cast<Real>(vnum))
    throw_error("Noninteger data stored as vertex global index");
  // mexPrintf("return\n");
  return vnum;
}

QMG::Real 
QMG::SimpComplex::real_coord_o(VertexOrdinalIndex vnumo, int jj) const {
#ifdef RANGECHECK
  if (jj < 0 || jj >= c_embedded_dim_) 
    throw_error("Range err in real_coord_o");
#endif
  // mexPrintf("reached real_coord_o vnumo = %d\n", vnumo);
  VertexOrdinalIndex vnumo1 = process_ordinal_vnum_(vnumo);
  // mexPrintf("returning\n");
  return c_nodearray_[vnumo1 *
                     (c_embedded_dim_ + 1) + jj + 1];
}

int 
QMG::SimpComplex::num_node_on_face(const Brep::Face_Spec& fspec) const {
  int fdim = fspec.fdim();
  int faceind = fspec.faceind();
#ifdef RANGECHECK
  if (fdim < 0 || fdim > c_gdim_)
    throw_error("fdim out of range");
#endif
  const mxArray* levarray = mxGetCell(c_entries_, 5+fdim);
#ifdef RANGECHECK
  if (faceind < 0 || faceind >= mxGetN(levarray))
    throw_error("faceind out of range");
#endif
  return mxGetN(mxGetCell(levarray, 2 * faceind));
}

QMG::SimpComplex::VertexGlobalIndex
QMG::SimpComplex::node_on_face(const Brep::Face_Spec& fspec,
                               int seqno) const {
  int fdim = fspec.fdim();
  int faceind = fspec.faceind();
#ifdef RANGECHECK
  if (fdim < 0 || fdim > c_gdim_)
    throw_error("fdim out of range");
#endif
  const mxArray* levarray = mxGetCell(c_entries_, 5+fdim);
#ifdef RANGECHECK
  if (faceind < 0 || faceind >= mxGetN(levarray))
    throw_error("faceind out of range");
#endif
  const mxArray* nodear = mxGetCell(levarray, 2 * faceind);
#ifdef RANGECHECK
  if (seqno < 0 || seqno >= mxGetN(nodear))
    throw_error("seqno out of range");
#endif
  double vnumr = mxGetPr(nodear)[seqno * mxGetM(nodear)];
  VertexGlobalIndex vnum = static_cast<int>(vnumr);
  if (static_cast<Real>(vnum) != vnumr)
    throw_error("Nonintegral data stored in node number slot");
  return vnum;
}

QMG::Brep::PatchIndex
QMG::SimpComplex::patchind_on_face(const Brep::Face_Spec& fspec,
                                   int seqno) const {
  int fdim = fspec.fdim();
  if (fdim == 0)
    return 0;
  int faceind = fspec.faceind();
#ifdef RANGECHECK
  if (fdim < 0 || fdim > c_gdim_)
    throw_error("fdim out of range");
#endif
  const mxArray* levarray = mxGetCell(c_entries_, 5+fdim);
#ifdef RANGECHECK
  if (faceind < 0 || faceind >= mxGetN(levarray))
    throw_error("faceind out of range");
#endif
  const mxArray* nodear = mxGetCell(levarray, 2 * faceind);
#ifdef RANGECHECK
  if (seqno < 0 || seqno >= mxGetN(nodear))
    throw_error("seqno out of range");
#endif
  double patchindr = mxGetPr(nodear)[seqno * mxGetM(nodear) + 1];
  Brep::PatchIndex patchind = static_cast<int>(patchindr);
  if (static_cast<Real>(patchind) != patchindr)
    throw_error("Nonintegral data stored in patch number slot");
  return patchind;
}

QMG::Real 
QMG::SimpComplex::param_coord_on_face(const Brep::Face_Spec& fspec, 
                                      int seqno, 
                                      int jj) const {
  int fdim = fspec.fdim();
  if (fdim == 0)
    return 0;
  int faceind = fspec.faceind();
#ifdef RANGECHECK
  if (fdim < 0 || fdim > c_gdim_)
    throw_error("fdim out of range");
  if (jj < 0 || jj >= fdim)
    throw_error("jj out of range");
#endif
  const mxArray* levarray = mxGetCell(c_entries_, 5+fdim);
#ifdef RANGECHECK
  if (faceind < 0 || faceind >= mxGetN(levarray))
    throw_error("faceind out of range");
#endif
  const mxArray* nodear = mxGetCell(levarray, 2 * faceind);
#ifdef RANGECHECK
  if (seqno < 0 || seqno >= mxGetN(nodear))
    throw_error("seqno out of range");
#endif
  return mxGetPr(nodear)[seqno * mxGetM(nodear) + jj + 2];
}


int 
QMG::SimpComplex::num_meshface_on_face(const Brep::Face_Spec& fspec) const {
  int fdim = fspec.fdim();
  int faceind = fspec.faceind();
#ifdef RANGECHECK
  if (fdim < 0 || fdim > c_gdim_)
    throw_error("fdim out of range");
#endif
  const mxArray* levarray = mxGetCell(c_entries_, 5+fdim);
#ifdef RANGECHECK
  if (faceind < 0 || faceind >= mxGetN(levarray))
    throw_error("faceind out of range");
#endif
  if (fdim == 0)
    return mxGetN(mxGetCell(levarray, 2 * faceind));
  return mxGetN(mxGetCell(levarray, 2 * faceind + 1));
}

QMG::SimpComplex::VertexGlobalIndex 
QMG::SimpComplex::node_of_meshface_on_face(const Brep::Face_Spec& fspec, 
                                           int seqno, 
                                           int jj) const {
  int fdim = fspec.fdim();
  int faceind = fspec.faceind();
#ifdef RANGECHECK
  if (fdim < 0 || fdim > c_gdim_)
    throw_error("fdim out of range");
  if (jj < 0 || jj > fdim)
    throw_error("jj  out of range");
#endif
  const mxArray* levarray = mxGetCell(c_entries_, 5+fdim);
#ifdef RANGECHECK
  if (faceind < 0 || faceind >= mxGetN(levarray))
    throw_error("faceind out of range");
#endif
  const mxArray* simparray = (fdim == 0)?
    mxGetCell(levarray, 2 * faceind) :
    mxGetCell(levarray, 2 * faceind + 1);
    
#ifdef RANGECHECK
  if (seqno < 0 || seqno >= mxGetN(simparray))
    throw_error("Seqno out of range");
#endif
  double vnumr = mxGetPr(simparray)[seqno * (fdim + 1) + jj];
  VertexGlobalIndex vnum = static_cast<int>(vnumr);
  if (static_cast<Real>(vnum) != vnumr) {
    mexPrintf("%d %e\n", vnum, vnumr);