qsimpcomp.cpp

来自「算断裂的」· C++ 代码 · 共 1,059 行 · 第 1/3 页

  max_node_surf_.resize(0);
  max_tri_surf_.resize(0);
  max_tet_region_.resize(0);
}



QMG::SimpComplex::VertexOrdinalIndex
QMG::SimpComplex_Under_Construction::add_vertex(const vector<Real>& realcoord, 
 VertexGlobalIndex thisid) {
  VertexOrdinalIndex nn = mesh_ -> num_nodes;
  if (nn >= max_nodes_) {
    int newsize = 2 * max_nodes_ + 1;
    Cptc_Mnode* new_nodes = 
      static_cast<Cptc_Mnode*>(malloc1(sizeof(Cptc_Mnode) * newsize));
    for (int jj = 0; jj < nn; ++jj)
      new_nodes[jj] = mesh_ -> nodes[jj];
    free(mesh_ -> nodes);
    mesh_ -> nodes = new_nodes;
    max_nodes_ = newsize;
  }
  mesh_ -> nodes[nn].id = thisid;
#ifdef RANGECHECK
  if (realcoord.size() != mesh_ -> embedded_dimension)
    throw_error("realcoord wrong size in add_vertex");
#endif
  for (int j = 0; j < mesh_ -> embedded_dimension; ++j)
    mesh_ -> nodes[nn].coord[j] = realcoord[j];
  ++(mesh_ -> num_nodes);
  pair<GlobalToOrdinalMapType::iterator, bool> rval =
    global_to_ordinal_ -> 
    insert(pair<VertexGlobalIndex, VertexOrdinalIndex>(thisid,nn));
  if (!rval.second)
    throw_error("global vertex id inserted twice");
  return nn;
}


#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();
  Brep::FaceIndex faceind = fspec.faceind();
  if (fdim < 0 || fdim > mesh_ -> intrinsic_dimension)
    throw_error("Face dimension out of range");
  if (faceind != this -> brep_levelsize(fdim))
    throw_error("Attempt to add faces out of order");
  if (fdim == 0) {
    if (faceind >= max_lev_size_[0]) {
      int newsize = faceind * 2 + 2;
      Cptc_MeshTopVertexInfo* newinfo = 
        static_cast<Cptc_MeshTopVertexInfo*>(malloc(sizeof(Cptc_MeshTopVertexInfo) * newsize));
      for (int j = 0; j < faceind; ++j)
        newinfo[j] = mesh_ -> vtxinfo[j];
      max_lev_size_[0] = newsize;
      free(mesh_ -> vtxinfo);
      mesh_ -> vtxinfo = newinfo;
    }
    ++(mesh_ -> num_top_vertices);
    mesh_ -> vtxinfo[faceind].vtx = -1;
  }
  else if (fdim == 1) {
    if (faceind >= max_lev_size_[1]) {
      int newsize = faceind * 2 + 2;
      Cptc_MeshTopEdgeInfo* newinfo = 
        static_cast<Cptc_MeshTopEdgeInfo*>(malloc(sizeof(Cptc_MeshTopEdgeInfo) * newsize));
      for (int j = 0; j < faceind; ++j)
        newinfo[j] = mesh_ -> edgeinfo[j];
      max_lev_size_[1] = newsize;
      free(mesh_ -> edgeinfo);
      mesh_ -> edgeinfo = newinfo;
    }
    max_node_edge_.resize(faceind + 1, 1);
    max_edge_edge_.resize(faceind + 1, 1);
    mesh_ -> edgeinfo[faceind].num_mnode = 0;
    mesh_ -> edgeinfo[faceind].num_medge = 0;
    mesh_ -> edgeinfo[faceind].incs =
        static_cast<Cptc_NodeCurveInc*>(malloc1(sizeof(Cptc_NodeCurveInc)));
    typedef VertexGlobalIndex Twoint[2];
    mesh_ -> edgeinfo[faceind].edgelist = 
      static_cast<Twoint*>(malloc1(sizeof(Twoint)));
    ++(mesh_ -> num_top_edges);
  }
  else if (fdim == 2) {
    if (faceind >= max_lev_size_[2]) {
      int newsize = faceind * 2 + 2;
      Cptc_MeshTopSurfaceInfo* newinfo = 
        static_cast<Cptc_MeshTopSurfaceInfo*>(malloc(sizeof(Cptc_MeshTopSurfaceInfo) * newsize));
      for (int j = 0; j < faceind; ++j)
        newinfo[j] = mesh_ -> surfinfo[j];
      max_lev_size_[2] = newsize;
      free(mesh_ -> surfinfo);
      mesh_ -> surfinfo = newinfo;
    }
    max_node_surf_.resize(faceind + 1, 1);
    max_tri_surf_.resize(faceind + 1, 1);
    mesh_ -> surfinfo[faceind].incs =
      static_cast<Cptc_NodePatchInc*>(malloc1(sizeof(Cptc_NodePatchInc)));
    mesh_ -> surfinfo[faceind].num_mnode = 0;
    typedef VertexGlobalIndex Threeint[3];
    mesh_ -> surfinfo[faceind].num_mtrifacet = 0;
    mesh_ -> surfinfo[faceind].trifacetlist = 
      static_cast<Threeint*>(malloc1(sizeof(Threeint)));
    mesh_ -> surfinfo[faceind].num_mquadfacet = 0;
    mesh_ -> surfinfo[faceind].quadfacetlist = 0;
    ++(mesh_ -> num_top_surfaces);
  }
  else if (fdim == 3) {
    if (faceind >= max_lev_size_[3]) {
      int newsize = faceind * 2 + 2;
      Cptc_MeshTopRegionInfo* newinfo = 
        static_cast<Cptc_MeshTopRegionInfo*>(malloc(sizeof(Cptc_MeshTopRegionInfo) * newsize));
      for (int j = 0; j < faceind; ++j)
        newinfo[j] = mesh_ -> regioninfo[j];
      max_lev_size_[3] = newsize;
      free(mesh_ -> regioninfo);
      mesh_ -> regioninfo = newinfo;
    }    
    max_tet_region_.resize(faceind + 1, 1);
    typedef VertexGlobalIndex Fourint[4];
    mesh_ -> regioninfo[faceind].tets = 
      static_cast<Fourint*>(malloc1(sizeof(Fourint)));
    mesh_ -> regioninfo[faceind].num_tet = 0;
    mesh_ -> regioninfo[faceind].num_prism = 0;
    mesh_ -> regioninfo[faceind].num_pyramid = 0;
    mesh_ -> regioninfo[faceind].num_hex = 0;
    mesh_ -> regioninfo[faceind].prisms = 0;
    mesh_ -> regioninfo[faceind].pyramids = 0;
    mesh_ -> regioninfo[faceind].hexs = 0;
    ++(mesh_ -> num_top_regions);
  }
}


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");
#ifdef RANGECHECK
  if (fdim < 0 || fdim >= mesh_ -> embedded_dimension || 
    fdim > mesh_ -> intrinsic_dimension)
    throw_error("Face dim out of range");
  if (faceind < 0)
    throw_error("Face index out of range");
  if ( (fdim == 0 && faceind >= mesh_ -> num_top_vertices) 
    || (fdim == 1 && faceind >= mesh_ -> num_top_edges)
    || (fdim == 2 && faceind >= mesh_ -> num_top_surfaces))
    throw_error("Face index out of range");
  if (param.size() != fdim) 
    throw_error("Wrong length for param coords");
#endif
  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) {
    // if (patchind != rval.first -> second.patchind)
    //   throw_error("Adding vertex twice; patchind mismatch");
    // for (int k2 = 0; k2 < fdim; ++k2) {
    //  Real p2 = rval.first -> second.param[k2];
    //  if (param[k2] != p2)
    //    throw_error("Adding vertex twice; param mismatch");
    // }
    return;
  }
  if (fdim == 0) {
    VertexGlobalIndex oldvnum = mesh_ -> vtxinfo[faceind].vtx;
    if (oldvnum < 0) 
      mesh_ -> vtxinfo[faceind].vtx = vnum;
    else if (oldvnum != vnum)
      throw_error("Distinct vertex numbers assigned to topological vertex");
  }
  else if (fdim == 1) {
    int nn = mesh_ -> edgeinfo[faceind].num_mnode;
    int mx = max_node_edge_[faceind];
    if (nn >= mx) {
      int newsize = 2 * mx + 1;
      Cptc_NodeCurveInc* newincs = 
        static_cast<Cptc_NodeCurveInc*>(malloc1(sizeof(Cptc_NodeCurveInc) * newsize));
      for (int jj = 0; jj < nn; ++jj)
        newincs[jj] = mesh_ -> edgeinfo[faceind].incs[jj];
      free(mesh_ -> edgeinfo[faceind].incs);
      mesh_ -> edgeinfo[faceind].incs = newincs;
      max_node_edge_[faceind] = newsize;
    }
    mesh_ -> edgeinfo[faceind].incs[nn].id = vnum;
    mesh_ -> edgeinfo[faceind].incs[nn].curve_index = patchind;
    mesh_ -> edgeinfo[faceind].incs[nn].param = param[0];
    ++(mesh_ -> edgeinfo[faceind].num_mnode);
  }
  else { // fdim == 2
    int nn = mesh_ -> surfinfo[faceind].num_mnode;
    int mx = max_node_surf_[faceind];
    if (nn >= mx) {
      int newsize = 2 * mx + 1;
      Cptc_NodePatchInc* newincs = 
        static_cast<Cptc_NodePatchInc*>(malloc1(sizeof(Cptc_NodePatchInc) * newsize));
      for (int jj = 0; jj < nn; ++jj)
        newincs[jj] = mesh_ -> surfinfo[faceind].incs[jj];
      free(mesh_ -> surfinfo[faceind].incs);
      mesh_ -> surfinfo[faceind].incs = newincs;
      max_node_surf_[faceind] = newsize;
    }
    mesh_ -> surfinfo[faceind].incs[nn].id = vnum;
    mesh_ -> surfinfo[faceind].incs[nn].patch_index = patchind;
    mesh_ -> surfinfo[faceind].incs[nn].param[0] = param[0];
    mesh_ -> surfinfo[faceind].incs[nn].param[1] = param[1];
    ++(mesh_ -> surfinfo[faceind].num_mnode);
  }
}

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");
  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 > mesh_ -> intrinsic_dimension)
    throw_error("Face dim out of range");
  if (faceind < 0)
    throw_error("Face index out of range");
  if ( (fdim == 3 && faceind >= mesh_ -> num_top_regions) 
    || (fdim == 1 && faceind >= mesh_ -> num_top_edges)
    || (fdim == 2 && faceind >= mesh_ -> num_top_surfaces))
    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");
  }
  if (fdim == 1) {
    Face_Bdry_Inc_ inc;
    if (verts[0] > verts[1]) {      
      inc.meshv[1] = verts[0];
      inc.meshv[0] = verts[1];
    }
    else {      
      inc.meshv[0] = verts[0];
      inc.meshv[1] = verts[1];
    }
    inc.meshv[2] = -1;
    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;
    }    
    int ne = mesh_ -> edgeinfo[faceind].num_medge;
    int mx = max_edge_edge_[faceind];
    if (ne >= mx) {
      int newsize = 2 * mx + 1;
      typedef VertexGlobalIndex Twoint[2];
      Twoint* newedges = 
        static_cast<Twoint*>(malloc1(2 * sizeof(Twoint) * newsize));
      for (int jj = 0; jj < ne; ++jj) {
        newedges[jj][0] = mesh_ -> edgeinfo[faceind].edgelist[jj][0];
        newedges[jj][1] = mesh_ -> edgeinfo[faceind].edgelist[jj][1];
      }
      free(mesh_ -> edgeinfo[faceind].edgelist);
      mesh_ -> edgeinfo[faceind].edgelist = newedges;
      max_edge_edge_[faceind] = newsize;
    }
    mesh_ -> edgeinfo[faceind].edgelist[ne][0] = verts[0];
    mesh_ -> edgeinfo[faceind].edgelist[ne][1] = verts[1];
    ++(mesh_ -> edgeinfo[faceind].num_medge);
    if (mesh_ -> intrinsic_dimension == 1)
      ++(mesh_ -> num_elems);
  }
  else if (fdim == 2) {
    Face_Bdry_Inc_ inc;
    for (int j = 0; j < 3; ++j)
      inc.meshv[j] = verts[j];
    if (inc.meshv[0] > inc.meshv[1]) {
      VertexGlobalIndex tmp = inc.meshv[0];
      inc.meshv[0] = inc.meshv[1];
      inc.meshv[1] = tmp;
    }
    if (inc.meshv[1] > inc.meshv[2]) {
      VertexGlobalIndex tmp = inc.meshv[1];
      inc.meshv[1] = inc.meshv[2];
      inc.meshv[2] = tmp;
    }
    if (inc.meshv[0] > inc.meshv[1]) {
      VertexGlobalIndex tmp = inc.meshv[0];