qsimpcomp.cpp

算断裂的

    throw_error("Nonintegral data stored in vnum slot");
  }
  return vnum;
}


QMG::SimpComplex_From_FrontEnd::
SimpComplex_From_FrontEnd(const mxArray* mesh) {
  mesh_ = mesh;
  if (strcmp(mxGetClassName(mesh), "zba") != 0)
    throw_error("Attempt to use non-zba array as simplicial complex");
  c_entries_ = mxGetFieldByNumber(mesh,0,0);
  if (c_entries_ == 0)
    throw_error("Input array has no field 0 (not a simplicial complex)");
  if (!mxIsCell(c_entries_) || mxGetM(c_entries_) < 6 || 
      mxGetN(c_entries_) != 1)
    throw_error("Input array should be columnvec cell zba with at least 6 entries for simplicial complex");
  
  const mxArray* typecode = mxGetCell(c_entries_, 0);
  if (typecode == 0)
    throw_error("Input array has no cell 0 (not a simplicial complex)");
  int lth1 = mxGetM(typecode) * mxGetN(typecode);
  if (lth1 != strlen(io_header_code()))
    throw_error("Typecode has wrong size -- object not a simplicial complex");
  char buf[100];
  if (mxGetString(typecode, buf, 100) ||
      compare_nocase(string(buf), string(io_header_code())) != 0)
    throw_error("Incorrect typecode -- object not a simplicial complex");
  const mxArray* gdima = mxGetCell(c_entries_, 1);
  if (gdima == 0 || !mxIsDouble(gdima) || mxIsSparse(gdima) || mxGetPi(gdima)
      || mxGetM(gdima) * mxGetN(gdima) != 1)
    throw_error("Missing or ill-formatted intrinsic dim field");
  Real gdimr = mxGetScalar(gdima);
  c_gdim_ = static_cast<int>(gdimr);
  if (static_cast<Real>(c_gdim_) != gdimr)
    throw_error("Nonintegral data for intrinsic dim");
  const mxArray* edima = mxGetCell(c_entries_, 2);
  if (edima == 0 || !mxIsDouble(edima) || mxIsSparse(gdima) || mxGetPi(gdima)
      || mxGetM(edima) * mxGetN(edima) != 1)
    throw_error("Missing or ill formatted embedded dim field");
  Real edimr = mxGetScalar(edima);
  c_embedded_dim_ = static_cast<int>(edimr);
  if (static_cast<Real>(c_embedded_dim_) != edimr)
    throw_error("Nonintegral data for embedded dim");
  if (c_embedded_dim_ != 2 && c_embedded_dim_ != 3)
    throw_error("Embedded dim must be 2 or 3");
  if (c_gdim_ < 1 || c_gdim_ > c_embedded_dim_)
    throw_error("Intrinsic dimension out of range");
  if (mxGetM(c_entries_) != 6 + c_gdim_) {
    mexPrintf(" c_gdim_ = %d mxGetM(c_entries_) = %d\n",
              c_gdim_, mxGetM(c_entries_));
    throw_error("Mesh cell array has wrong number of entries");
  }
  const mxArray* vertices = mxGetCell(c_entries_, 4);
  if (vertices == 0 || !mxIsDouble(vertices) || mxIsSparse(vertices)
      || mxGetPi(vertices))
    throw_error("Missing or illformatted vertices");
  c_nodearray_ = mxGetPr(vertices);
  if (mxGetN(vertices) > 0 && mxGetM(vertices) != 1 + c_embedded_dim_)
    throw_error("Node array has incorrect number of rows");
  if (mxGetN(vertices) > 0 && c_nodearray_ == 0)
    throw_error("Unable to retrieve node array");
  const mxArray* pvlist = mxGetCell(c_entries_, 3);
  if (pvlist == 0 || !mxIsCell(pvlist) || 
      (mxGetN(pvlist) > 0 && mxGetM(pvlist) != 2))
    throw_error("Property-value list must be a 2-by-k cell array");
  for (int facedim = 0; facedim <= c_gdim_; ++facedim) {
    const mxArray* levarray = mxGetCell(c_entries_, 5 + facedim);
    if (levarray == 0 || !mxIsCell(levarray))
      throw_error("Level array of mesh must be cell array");
    int levsize = mxGetN(levarray);
    if (levsize > 0 && mxGetM(levarray) != 2)
      throw_error("Level array of mesh must be 2-by-k cell array");
    for (Brep::FaceIndex faceind = 0; faceind < levsize; ++faceind) {
      const mxArray* vtxlist = mxGetCell(levarray, 2 * faceind);
      if (vtxlist == 0 || !mxIsDouble(vtxlist) || mxIsSparse(vtxlist))
        throw_error("Vertex list of brep face must be a matrix");
      int numv = mxGetN(vtxlist);
      if (numv > 0 && mxGetM(vtxlist) != 
          1 + ((facedim == 0)? 0 : 1) + facedim)
        throw_error("Vertex list of brep face has wrong number of rows");
      const mxArray* smplist = mxGetCell(levarray, 2 * faceind + 1);
      if (smplist == 0 || !mxIsDouble(smplist) || mxIsSparse(smplist))
        throw_error("Simplex list of brep face must be a matrix");
      int nums = mxGetN(smplist);
      if (nums > 0 && mxGetM(smplist) != facedim + 1)
        throw_error("Vertex list of brep face has wrong number of rows");
    }
  }
  global_to_ordinal_ = 0;
}

void
QMG::SimpComplex_From_FrontEnd::make_global_to_ordinal_() const {
  if (global_to_ordinal_) 
    throw_error("Make global to ordinal called twice");
  global_to_ordinal_ = new GlobalToOrdinalMapType;
  int nn = mxGetN(mxGetCell(c_entries_, 4));
  for (VertexOrdinalIndex vnumo = 0;
       vnumo < nn;
       ++vnumo) {
    double vnumr = c_nodearray_[vnumo * (c_embedded_dim_ + 1)];
    VertexGlobalIndex vnum = static_cast<int>(vnumr);
    if (vnumr != static_cast<Real>(vnum))
      throw_error("Nonintegral data as vertex global index");
    if (vnum < 0)
      throw_error("global vertex id out of range");
    pair<GlobalToOrdinalMapType::iterator, bool> rval =
      global_to_ordinal_ -> 
      insert(pair<VertexGlobalIndex, VertexOrdinalIndex>(vnum,vnumo));
    if (!rval.second)
      throw_error("global vertex id occurs twice in input");
  }
}

QMG::SimpComplex_From_FrontEnd::~SimpComplex_From_FrontEnd() {
  delete global_to_ordinal_;
}


QMG::SimpComplex_From_FrontEnd::
SimpComplex_From_FrontEnd(SimpComplex_From_FrontEnd& other) :
  SimpComplex(other) {
  other.global_to_ordinal_ = 0;
}


QMG::SimpComplex_From_FrontEnd::
SimpComplex_From_FrontEnd(const SimpComplex_From_FrontEnd_Returnval& sucr) :
SimpComplex(sucr.sref_) {
  sucr.sref_.global_to_ordinal_ = 0;
}


bool 
QMG::SimpComplex_Under_Construction::Face_Bdry_Inc_::
operator<(const Face_Bdry_Inc_& other) const {
  if (meshv[0] < other.meshv[0]) return true;
  if (meshv[0] > other.meshv[0]) return false;
  if (meshv[1] < other.meshv[1]) return true;
  if (meshv[1] > other.meshv[1]) return false;
  return meshv[2] < other.meshv[2];
}

bool
QMG::SimpComplex_Under_Construction::Vertex_Bdry_Inc_::
operator<(const Vertex_Bdry_Inc_& other) const {
  if (vnum < other.vnum) return true;
  if (vnum > other.vnum) return false;
  return fspec < other.fspec;
}

void 
QMG::SimpComplex_Under_Construction::clobber_() {
  mesh_ = 0;
  global_to_ordinal_ = 0;
  vertinclookup_ = 0;  
  faceinclookup_ = 0;
  max_elems_per_face_ = 0;
  max_nodes_per_face_ = 0;
}

  
QMG::SimpComplex::VertexOrdinalIndex 
QMG::SimpComplex_Under_Construction::add_vertex(const vector<Real>& realcoord, 
                                                VertexGlobalIndex vnum) {
#ifdef RANGECHECK
  if (realcoord.size() != c_embedded_dim_)
    throw_error("Real coord wrong size in addvertex");
#endif

  mxArray* node_array = mxGetCell(c_entries_, 4); 
  int numnode = mxGetN(node_array);
  int p = numnode * (c_embedded_dim_ + 1);
  if (numnode == max_nodes_) {
    int newsize = numnode * 2 + 2;
    double* newpr = 
      static_cast<double*>(mxMalloc(newsize * (c_embedded_dim_ + 1) * sizeof(double)));
    for (int j = 0; j < p; ++j)
      newpr[j] = c_nodearray_[j];
    mxFree(const_cast<double*>(c_nodearray_));
    c_nodearray_ = newpr;
    mxSetPr(node_array, newpr);
    max_nodes_ = newsize;
  }
  double* nodearray1 = const_cast<double*>(c_nodearray_);
  nodearray1[p] = static_cast<double>(vnum);
  mxSetN(node_array, numnode + 1);
  for (int j1 = 0; j1 < c_embedded_dim_; ++j1)
    nodearray1[++p] = realcoord[j1];
  pair<GlobalToOrdinalMapType::iterator, bool> rval =
    global_to_ordinal_ -> 
    insert(pair<VertexGlobalIndex, VertexOrdinalIndex>(vnum,numnode));
  if (!rval.second)
    throw_error("global vertex id inserted twice");
  return numnode;
}


#ifndef BUG_IN_USING
using QMG::pair;
#endif

pair<QMG::SimpComplex::VertexOrdinalIndex, QMG::SimpComplex::VertexGlobalIndex> 
QMG::SimpComplex_Under_Construction::
add_vertex(const vector<Real>& realcoord) {
  VertexGlobalIndex thisid = next_vnum_;
  VertexOrdinalIndex vnum_o = add_vertex(realcoord, thisid);
  ++next_vnum_;
  return pair<VertexOrdinalIndex,VertexGlobalIndex>(vnum_o, thisid);
}


void 
QMG::SimpComplex_Under_Construction::add_fspec(const Brep::Face_Spec& fspec) {
  int fdim = fspec.fdim();
  if (fdim < 0 || fdim > c_gdim_)
    throw_error("Face dim out of range");
  const mxArray* lev1 = mxGetCell(c_entries_, 5 + fdim);
  mxArray* nclev1 = const_cast<mxArray*>(lev1);
  Brep::FaceIndex faceind = fspec.faceind();
  if (mxGetN(lev1) != faceind)
    throw_error("Faces added out of order");
  if (faceind == max_lev_size_[fdim]) {
    int newsize = faceind * 2 + 2;
    mxArray** oldcells = static_cast<mxArray**>(mxGetData(lev1));
    mxArray** newcells = 
      static_cast<mxArray**>(mxMalloc(2 * newsize * sizeof(mxArray*)));
    for (int j = 0; j < 2 * faceind; ++j)
      newcells[j] = oldcells[j];
    mxFree(oldcells);
    mxSetData(nclev1, static_cast<void*>(newcells));
    max_lev_size_[fdim] = newsize;
  }
  int numrows = 1 + ( (fdim == 0)? 0 : 1) + fdim;
  mxArray* newv = mxCreateDoubleMatrix(numrows, 1, mxREAL);
  mxSetN(newv, 0);
  (*max_nodes_per_face_)[fspec] = 1;
  mxSetN(nclev1, faceind + 1);
  mxSetCell(nclev1, faceind * 2, newv);

  mxArray* newe = mxCreateDoubleMatrix(fdim + 1, 1,  mxREAL);
  mxSetN(newe, 0);
  (*max_elems_per_face_)[fspec] = 1;
  mxSetCell(nclev1, faceind * 2 + 1, newe);
}


void 
QMG::SimpComplex_Under_Construction::
add_vertex_bdryinc(VertexGlobalIndex vnum,
                   const Brep::Face_Spec& fspec,
                   Brep::PatchIndex patchind,
                   const vector<Real>& param) {
  int fdim = fspec.fdim();
  Brep::FaceIndex faceind = fspec.faceind();
  if (global_to_ordinal_ -> find(vnum) == global_to_ordinal_ -> end())
    throw_error("Global vertex id not found");
  const mxArray* lev1 = mxGetCell(c_entries_, 5 + fdim);
#ifdef RANGECHECK
  if (fdim < 0 || fdim >= c_embedded_dim_ ||
      fdim > c_gdim_)
    throw_error("Face dim out of range");
  if (faceind < 0 || faceind >= mxGetN(lev1))
    throw_error("Face index out of range");
#endif
  if (param.size() != fdim) 
    throw_error("Wrong length for param coords");
  if (patchind < 0)
    throw_error("Patch index out of range");
  Vertex_Bdry_Inc_ inc;
  Vertex_Pos_ pos;
  inc.vnum = vnum;
  inc.fspec = fspec;
  for (int k1 = 0; k1 < fdim; ++k1)
    pos.param[k1] = param[k1];
  pos.patchind = patchind;
  pair<Vert_Bdry_Lookup_::iterator,bool> rval = 
    vertinclookup_ -> insert(pair<Vertex_Bdry_Inc_, Vertex_Pos_>(inc,pos));
  if (!rval.second) {
    return;
  }
  const mxArray* vtab1 = mxGetCell(lev1, 2 * faceind);
  mxArray* ncvtab1 = const_cast<mxArray*>(vtab1);
  int nn = mxGetN(vtab1);
  int numrow = 1 + ((fdim == 0)? 0 : 1) + fdim;
  int p = numrow * nn;
  int oldmax = (*max_nodes_per_face_)[fspec];
  if (nn == oldmax) {
    int newsize = oldmax * 2 + 2;
    double* newpr = 
      static_cast<double*>(mxMalloc(newsize * numrow * sizeof(double)));
    double* oldpr = mxGetPr(vtab1);
    for (int i = 0; i < p; ++i)
      newpr[i] = oldpr[i];
    mxFree(oldpr);
    mxSetPr(ncvtab1, newpr);
    (*max_nodes_per_face_)[fspec] = newsize;
  }
  mxSetN(ncvtab1, nn + 1);
  double* pr = mxGetPr(vtab1);
  pr[p] = static_cast<double>(vnum);
  if (fdim > 0)
    pr[++p] = static_cast<double>(patchind);
  for (int j = 0; j < fdim; ++j)
    pr[++p] = param[j];
}



QMG::Brep::PatchIndex 
QMG::SimpComplex_Under_Construction::
retrieve_vertex_bdryinc_patchind(VertexGlobalIndex vnum,
                                 const Brep::Face_Spec& fspec) const {
  Vertex_Bdry_Inc_ vp;
  vp.vnum = vnum;
  vp.fspec = fspec;
  Vert_Bdry_Lookup_::const_iterator it = 
    vertinclookup_ -> find(vp);
  if (it == vertinclookup_ -> end())
    throw_error("Vertex incidence not found");
  return it -> second.patchind;
}


QMG::Point
QMG::SimpComplex_Under_Construction::
retrieve_vertex_bdryinc_param(VertexGlobalIndex vnum,
                              const Brep::Face_Spec& fspec) const {
  Vertex_Bdry_Inc_ vp;
  vp.vnum = vnum;
  vp.fspec = fspec;
  Vert_Bdry_Lookup_::const_iterator it = 
    vertinclookup_ -> find(vp);
  if (it == vertinclookup_ -> end())
    throw_error("Vertex incidence not found");