qsimpcomp.cpp

算断裂的

  Point pt;
  pt[0] = it -> second.param[0];
  pt[1] = it -> second.param[1];
  return pt;
}


void
QMG::SimpComplex_Under_Construction::
add_simplex_face(const vector<VertexGlobalIndex>& verts,
                 const Brep::Face_Spec& fspec) { 
  int fdim = fspec.fdim();
  Brep::FaceIndex faceind = fspec.faceind();

#ifdef RANGECHECK
  if (fdim <= 0 || fdim > c_gdim_)
    throw_error("Face dim out of range");
#endif

  const mxArray* lev1 = mxGetCell(c_entries_, 5 + fdim);

#ifdef RANGECHECK
  if (faceind < 0 || faceind >= mxGetN(lev1))
    throw_error("Face index out of range");
  if (verts.size() != fdim + 1) 
    throw_error("Wrong length for vertex list");
#endif

  for (int jj = 0; jj < fdim + 1; ++jj) {
    if (global_to_ordinal_ -> find(verts[jj]) == global_to_ordinal_ -> end())
      throw_error("Global vertex id not found");
  }
  Face_Bdry_Inc_ inc;
  if (fdim < c_embedded_dim_) {
    {
      for (int j1 = 0; j1 < fdim + 1; ++j1)
        inc.meshv[j1] = verts[j1];
    }
    {
      for (int j1 = fdim + 1; j1 < 3; ++j1)
        inc.meshv[j1] =  -1;
    }
    {
      for (int j1 = 0; j1 < fdim; ++j1) {
        for (int j2 = j1 + 1; j2 < fdim + 1; ++j2) {
          if (inc.meshv[j2] < inc.meshv[j1]) {
            VertexGlobalIndex tmp = inc.meshv[j2];
            inc.meshv[j2] = inc.meshv[j1];
            inc.meshv[j1] = tmp;
          }
        }
      }
    }
    pair<Face_Bdry_Lookup_::iterator, bool> rval =
      faceinclookup_ -> insert(pair<Face_Bdry_Inc_,Brep::Face_Spec>(inc,fspec));
    if (!rval.second) {
      if (rval.first -> second != fspec) {
        throw_error("Attempt to add simplex face with mismatched fspecs");
      }
      return;
    }
  }
  const mxArray* etab1 = mxGetCell(lev1, 2 * faceind + 1);
  mxArray* ncetab1 = const_cast<mxArray*>(etab1);
  int ne = mxGetN(etab1);
  int numrow = fdim + 1;
  int p = numrow * ne;
  int oldmax = (*max_elems_per_face_)[fspec];
  if (ne == oldmax) {
    int newsize = oldmax * 2 + 2;
    double* newpr = 
      static_cast<double*>(mxMalloc(newsize * numrow * sizeof(double)));
    double* oldpr = mxGetPr(etab1);
    for (int i = 0; i < p; ++i)
      newpr[i] = oldpr[i];
    mxFree(oldpr);
    mxSetPr(ncetab1, newpr);
    (*max_elems_per_face_)[fspec] = newsize;
  }
  mxSetN(ncetab1, ne + 1);
  double* pr = mxGetPr(etab1);
  for (int j = 0; j < fdim + 1; ++j)
    pr[p++] = static_cast<double>(verts[j]);
  //dump610(c_entries_);
}


void
QMG::SimpComplex_Under_Construction::remove_unused_vertices() {
  map<VertexGlobalIndex, int> lookupnode;
  int levsize = mxGetN(mxGetCell(c_entries_, 5 + c_gdim_));
  {
    for (Brep::FaceIndex faceind = 0; faceind < levsize; ++faceind) {
      Brep::Face_Spec fspec(c_gdim_, faceind);
      int nf = this -> num_meshface_on_face(fspec);
      for (int seqno = 0; seqno < nf; ++seqno) {
        for (int l = 0; l < c_gdim_ + 1; ++l) {
          VertexGlobalIndex vnum = 
            this -> node_of_meshface_on_face(fspec, seqno, l);
          lookupnode[vnum] = 1;
        }
      }
    }
  }
  VertexOrdinalIndex dest_vnumo = 0;
  mxArray* node_array = mxGetCell(c_entries_, 4);
  int nn = mxGetN(node_array);
  double* nodepr = mxGetPr(node_array);
  int numrows = c_embedded_dim_ + 1;
  for (VertexOrdinalIndex source_vnumo = 0;
  source_vnumo < nn;
  ++source_vnumo) {
    VertexGlobalIndex vnum = 
      static_cast<int>(nodepr[source_vnumo * numrows]);
    if (lookupnode[vnum]) {
      for (int jj = 0; jj < numrows; ++jj)
        nodepr[dest_vnumo * numrows + jj] = 
          nodepr[source_vnumo * numrows + jj];
      ++dest_vnumo;
    }
  }
  mxSetN(node_array, dest_vnumo);


  for (int fdim = 0; fdim < c_gdim_; ++fdim) {
    const mxArray* lev1 = mxGetCell(c_entries_, fdim + 5);
    int levsize = mxGetN(lev1);
    for (Brep::FaceIndex faceind = 0; faceind < levsize; ++faceind) {
      Brep::Face_Spec fspec(fdim, faceind);
      mxArray* facenodes = mxGetCell(lev1, 2 * faceind);
      double* fnr = mxGetPr(facenodes);
      int nnf = mxGetN(facenodes);
      int numr1 = 1 + ((fdim == 0)? 0:1) + fdim;
      int destseqno = 0;
      for (int seqno = 0; seqno < nnf; ++seqno) {
        VertexGlobalIndex vnum = 
          static_cast<int>(fnr[seqno * numr1]);
        if (lookupnode[vnum]) {
          for (int j = 0; j < numr1; ++j)
            fnr[destseqno * numr1 + j] = fnr[seqno * numr1 + j];
          ++destseqno;
        }
      }
      mxSetN(facenodes, destseqno);
      if (fdim == 0) continue;
      mxArray* simpnodes = mxGetCell(lev1, 2 * faceind + 1);
      double* snr = mxGetPr(simpnodes);
      int nmf = mxGetN(simpnodes);
      numr1 = fdim + 1;
      int destseqno2 = 0;
      for (int seqno2 = 0; seqno2 < nmf; ++seqno2) {
        bool all_still_there = true;
        for (int jj = 0; jj < fdim + 1; ++jj) {
          VertexGlobalIndex vnum = static_cast<int>(snr[seqno2 * numr1 + jj]);
          if (!lookupnode[vnum]) {
            all_still_there = false;
            break;
          }
        }
        if (all_still_there) {
          for (int j2 = 0; j2 < fdim + 1; ++j2)
            snr[destseqno2 * numr1 + j2] = snr[seqno2 * numr1 + j2];
          ++destseqno2;
        }
      }
      mxSetN(simpnodes, destseqno2);
    }
  }
}

void 
QMG::SimpComplex_Under_Construction::add_prop_val(const string& prop, 
                                                  const string& val) {
  mxArray* pvlist = mxGetCell(c_entries_, 3);
  int npv = mxGetN(pvlist);
  if (npv == max_prop_val_) {
    int newsize = npv * 2 + 2;
    mxArray** oldcells = static_cast<mxArray**>(mxGetData(pvlist));
    mxArray** newcells = 
      static_cast<mxArray**>(mxMalloc(2 * newsize * sizeof(mxArray*)));
    for (int j = 0; j < 2 * npv; ++j)
      newcells[j] = oldcells[j];
    mxFree(oldcells);
    mxSetData(pvlist, static_cast<void*>(newcells));
    max_prop_val_ = newsize;
  }
  mxSetN(pvlist, npv + 1);
  mxSetCell(pvlist, 2 * npv, mxCreateString(prop.c_str()));
  mxSetCell(pvlist, 2 * npv + 1, mxCreateString(val.c_str()));
}



QMG::SimpComplex_Under_Construction::
SimpComplex_Under_Construction(int intrinsic_dim, 
			       int embedded_dim, 
                               const vector<int>& level_size) {
  if (embedded_dim != 2 && embedded_dim != 3)
    throw_error("Illegal embedded dimension");
  if (intrinsic_dim < 1 || intrinsic_dim > embedded_dim)
    throw_error("Illegal intrinsic dimension");
  int brep_gdim = level_size.size() - 1;
  if (brep_gdim < intrinsic_dim || brep_gdim > embedded_dim)
    throw_error("illegal brep intrinsic dimension");


  mxArray* plhsa[1];
  mxArray* prhsa[1];
  mexCallMATLAB(1,plhsa,0, prhsa, "make_empty_zba");
  mesh_ = plhsa[0];
  mxArray* oldentries = mxGetFieldByNumber(mesh_,0, 0);
  if (oldentries)
    mxFree(oldentries);

  c_gdim_ = intrinsic_dim;
  c_embedded_dim_ = embedded_dim;

  mxArray* entries = mxCreateCellMatrix(6+intrinsic_dim,1);
  c_entries_ = entries;
  mxSetFieldByNumber(const_cast<mxArray*>(mesh_),0,0,entries);


  mxArray* typecode = mxCreateString(io_header_code());
  mxSetCell(entries, 0, typecode);


  mxArray* gdima = mxCreateDoubleMatrix(1,1, mxREAL);
  *mxGetPr(gdima) = static_cast<double>(intrinsic_dim);
  mxSetCell(entries, 1, gdima);

  mxArray* edima = mxCreateDoubleMatrix(1,1,mxREAL);
  *mxGetPr(edima) = static_cast<double>(embedded_dim);
  mxSetCell(entries, 2, edima);

  mxArray* pvlist = mxCreateCellMatrix(2,1);
  mxSetN(pvlist, 0);
  max_prop_val_ = 1;
  mxSetCell(entries, 3, pvlist);

  mxArray* nodearray = mxCreateDoubleMatrix(1 + c_embedded_dim_, 1, mxREAL);
  mxSetN(nodearray, 0);
  max_nodes_ = 1;
  mxSetCell(entries, 4, nodearray);
  c_nodearray_ = mxGetPr(nodearray);


  max_lev_size_.resize(c_gdim_ + 1);
  max_nodes_per_face_ = new map<Brep::Face_Spec, int>;
  max_elems_per_face_ = new map<Brep::Face_Spec, int>;
  for (int fdim = 0; fdim <= c_gdim_; ++fdim) {
    int lsize = level_size[fdim];
    if (lsize < 0)
      throw_error("Level size out of bounds");
    int maxlevsize = (lsize == 0)? 1 : lsize;
    max_lev_size_[fdim] = maxlevsize;
    mxArray* lev = mxCreateCellMatrix(2, maxlevsize);
    mxSetN(lev, lsize);
    for (Brep::FaceIndex faceind = 0; faceind < lsize; ++faceind) {
      Brep::Face_Spec fspec(fdim, faceind);
      mxArray* nodes = 
        mxCreateDoubleMatrix(1 + ((fdim == 0)? 0 : 1) + fdim, 1, mxREAL);
      mxSetN(nodes, 0);
      mxSetCell(lev, 2 * faceind, nodes);
      (*max_nodes_per_face_)[fspec] = 1;
      mxArray* simps = 
        mxCreateDoubleMatrix(fdim + 1, 1, mxREAL);
      mxSetN(simps, 0);
      mxSetCell(lev, 2 * faceind + 1, simps);
      (*max_elems_per_face_)[fspec] = 1;
    }
    mxSetCell(entries, 5 + fdim, lev);
  }
  vertinclookup_ = new Vert_Bdry_Lookup_;
  faceinclookup_ = new Face_Bdry_Lookup_;
  global_to_ordinal_ = new GlobalToOrdinalMapType;
  next_vnum_ = 0;
}

mxArray*
QMG::SimpComplex_Under_Construction::release() {
  mxArray* tmp = const_cast<mxArray*>(mesh_);
  delete faceinclookup_;
  delete vertinclookup_;
  delete global_to_ordinal_;
  delete max_elems_per_face_;
  delete max_nodes_per_face_;
  clobber_();
  return tmp;
}


QMG::SimpComplex_Under_Construction::~SimpComplex_Under_Construction() {
  if (mesh_)
    mxDestroyArray(const_cast<mxArray*>(mesh_));
  delete faceinclookup_;
  delete vertinclookup_;
  delete global_to_ordinal_;
  delete max_elems_per_face_;
  delete max_nodes_per_face_;
}
  

QMG::SimpComplex_Under_Construction::
SimpComplex_Under_Construction(SimpComplex_Under_Construction& other) :
SimpComplex(other),
max_nodes_per_face_(other.max_nodes_per_face_),
max_elems_per_face_(other.max_elems_per_face_),
max_lev_size_(other.max_lev_size_),
max_nodes_(other.max_nodes_),
max_prop_val_(other.max_prop_val_),
next_vnum_(other.next_vnum_),
vertinclookup_(other.vertinclookup_),
faceinclookup_(other.faceinclookup_) {
  other.clobber_();
}


QMG::SimpComplex_Under_Construction::
SimpComplex_Under_Construction(const SimpComplex_Under_Construction_Returnval& sucr) :
SimpComplex(sucr.sref_),
max_nodes_per_face_(sucr.sref_.max_nodes_per_face_),
max_elems_per_face_(sucr.sref_.max_elems_per_face_),
max_lev_size_(sucr.sref_.max_lev_size_),
max_nodes_(sucr.sref_.max_nodes_),
max_prop_val_(sucr.sref_.max_prop_val_),
next_vnum_(sucr.sref_.next_vnum_),
vertinclookup_(sucr.sref_.vertinclookup_),
faceinclookup_(sucr.sref_.faceinclookup_) {
  sucr.sref_.clobber_();
}