gmboundary.cpp

来自「算断裂的」· C++ 代码 · 共 113 行

// ------------------------------------------------------------------
// gmboundary.cpp
//
// This file contains the driver for gmboundary, which 
// applies the boundary operation to a mesh.
// ------------------------------------------------------------------
// 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 "qbrep.h"
#include "qsimpcomp.h"
#include "qfrontend.h"

#define GM_ROUTINE_NAME gmboundary
#define GM_ROUTINE_NAME_Q "gmboundary"


namespace {
  using namespace QMG;
  using namespace QMG::FrontEnd;

  SimpComplex_Under_Construction get_sc_boundary(const SimpComplex& sc) {
    int gdim = sc.gdim();
    if (gdim == 1)
      throw_error("Cannot compute boundary of 1d simplicial complex");
    int di = sc.embedded_dim();
    int newgdim = gdim - 1;
    vector<int> newlevsize(newgdim + 1);
    {
      for (int fdim = 0; fdim <= newgdim; ++fdim) 
        newlevsize[fdim] = sc.brep_levelsize(fdim);
    }
    SimpComplex_Under_Construction newsc(newgdim, di, newlevsize);
    int npv = sc.num_prop_val();
    for (int jj = 0; jj < npv; ++jj)
      newsc.add_prop_val(sc.prop(jj), sc.val(jj));
    map<SimpComplex::VertexGlobalIndex, int> vertex_appears_on_bdry;
    {
      for (int fdim = 0; fdim <= newgdim; ++fdim) {
        for (Brep::FaceIndex faceind = 0; faceind < newlevsize[fdim]; ++faceind) {
          Brep::Face_Spec fspec(fdim, faceind);
          int nnf = sc.num_node_on_face(fspec);
          for (int seqno = 0; seqno < nnf; ++seqno)
            vertex_appears_on_bdry[sc.node_on_face(fspec,seqno)] = 1;
        }
      }
    }
    vector<Real> coord(di);
    int nnode = sc.num_nodes();
    for (SimpComplex::VertexOrdinalIndex vnum_o = 0; vnum_o < nnode; ++vnum_o) {
      SimpComplex::VertexGlobalIndex vnum = sc.ordinal_to_global(vnum_o);
      if (!vertex_appears_on_bdry[vnum]) continue;
      for (int jj = 0; jj < di; ++jj)
        coord[jj] = sc.real_coord_o(vnum_o, jj);
      newsc.add_vertex(coord,vnum);
    }
    {
      vector<Real> param(di);
      vector<SimpComplex::VertexGlobalIndex> simp(di);
      for (int fdim = 0; fdim <= newgdim; ++fdim) {
        param.resize(fdim);
        simp.resize(fdim + 1);
        for (Brep::FaceIndex faceind = 0; faceind < newlevsize[fdim]; ++faceind) {
          Brep::Face_Spec fspec(fdim, faceind);
          int nnf = sc.num_node_on_face(fspec);
          for (int seqno = 0; seqno < nnf; ++seqno) {
            SimpComplex::VertexOrdinalIndex vnum = sc.node_on_face(fspec,seqno);
            for (int jj = 0; jj < fdim; ++jj)
              param[jj] = sc.param_coord_on_face(fspec,seqno,jj);
            Brep::PatchIndex patchind = sc.patchind_on_face(fspec,seqno);
            newsc.add_vertex_bdryinc(vnum, fspec, patchind, param);
          }
          if (fdim == 0) continue;
          int nmf = sc.num_meshface_on_face(fspec);
          for (int seqno2 = 0; seqno2 < nmf; ++seqno2) {
            for (int jj = 0; jj < fdim + 1; ++jj)
              simp[jj] = sc.node_of_meshface_on_face(fspec, seqno2, jj);
            newsc.add_simplex_face(simp,fspec);
          }
        }
      }
    }
    return newsc;
  }




  static void worker(const QMG::FrontEnd::ArgValType& argvalhandle,
    QMG::FrontEnd::ReturnValType& returnvalhandle,
    QMG::FrontEnd::Interpreter& interp) {
    
    argvalhandle.verify_nargin(1, 1, GM_ROUTINE_NAME_Q);
    returnvalhandle.verify_nargout(1, 1, GM_ROUTINE_NAME_Q);
    SimpComplex_Under_Construction newsc = 
      get_sc_boundary(argvalhandle.get_simpcomp(0));
    returnvalhandle.put_simpcomp(0, newsc);
  }
}

#include "gm_entrypoint.cpp"