gmapply.cpp

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

// ------------------------------------------------------------------
// gmapply.cpp
//
// This file contains the driver and worker function for gmapply, which 
// applies a linear affine transformation to a brep or 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.
// ------------------------------------------------------------------

#ifdef _MSC_VER
#if _MSC_VER == 1200
#include "qmatvec.h"
#endif
#endif

#include "qbrep_constr.h"
#include "qsimpcomp.h"
#include "qfrontend.h"
#include "qmatvec.h"

#define GM_ROUTINE_NAME gmapply
#define GM_ROUTINE_NAME_Q "gmapply"


namespace {
  using namespace QMG;



  //  void mexpause(const char* s) {
  //    mexPrintf("%s\n",s);
  // }

  void matrix_vector_mult(const Matrix& A, 
    const vector<Real>& x,
    vector<Real>& answer) {
    
    
    int m = A.numrows();
    int n = A.numcols();
#ifdef DEBUGGING
    if (n != x.size() || answer.size() != m)
      throw_error("Size mismatch in mvm");
#endif
    
    for (int i = 0; i < m; ++i) {
      Real t = 0.0;
      for (int j = 0; j < n; ++j)
        t += A(i,j)*x[j];
      answer[i] = t;
    }
  }
  
  
  // ------------------------------------------------------------------
  // apply a transformation to a brep -- means applying transformation
  // to control point coordinates.
 
  
  Brep_Under_Construction apply_transformation(const Brep& g, 
    const Matrix& m) {
    int di = g.embedded_dim();
    int gdim = g.gdim();
    if (m.numcols() != di + 1)
      throw_error("Tranformation has the wrong number of columns");
    int newdi = m.numrows();
    if (newdi != di) 
      throw_error("Applying a transform to a brep cannot change the embedded dimension");


    Brep_Under_Construction newg(gdim, newdi);
    Brep_Under_Construction::Propval_inserter pvi(newg);
    for (Brep::Loop_over_propvals_of_brep loop(g);
    loop.notdone();
    ++loop) {
      pvi.insert_propval(loop.prop(),loop.val());
    }

    Brep_Under_Construction::Control_point_inserter cpi(newg);

    vector<Real> destpt(newdi);
    vector<Real> srcpt(di + 1);

    for (int cpind = 0; cpind < g.num_control_points(); ++cpind) {
      for (int j = 0; j < di; ++j)
        srcpt[j] = g.control_point(cpind, j);
      srcpt[di] = -1.0;
      matrix_vector_mult(m, srcpt, destpt);
      cpi.insert_control_point(destpt);
    }

    for (int dim = 0; dim <= gdim; ++dim) {
      Brep_Under_Construction::Top_face_inserter tfi(newg, dim);
      for (Brep::Face_Spec_Loop_Over_Faces_Of_Dim fspec(g, dim);
      fspec.notdone();
      ++fspec) {
        tfi.insert_top_face(g.face_name(fspec));
        Brep_Under_Construction::Top_face_propval_inserter tpvi(newg, fspec);
        for (Brep::Loop_over_propvals_of_face pvloop(g, fspec);
        pvloop.notdone();
        ++pvloop) {
          tpvi.insert_propval(pvloop.prop(), pvloop.val());
        }
        
        Brep_Under_Construction::Top_face_bound_inserter tfbi(newg, fspec);
        Brep_Under_Construction::Top_face_ib_inserter tfii(newg, fspec); 
        for (Brep::Face_Spec_Loop_Over_Face_Subfaces subfspec(g, fspec);
        subfspec.notdone();
        ++subfspec) {
          if (subfspec.is_internal_boundary())
            tfii.insert_ib(subfspec);
          else
            tfbi.insert_boundary(subfspec, subfspec.orientation());
        }
        
        
        vector<Brep::ControlPointIndex> cp;

        Brep_Under_Construction::Top_face_patch_inserter pti(newg, fspec);
        for (Brep::Loop_over_patches_of_face ptloop(g, fspec);
        ptloop.notdone();
        ++ptloop) {
          int deg1 = ptloop.degree1();
          int deg2 = 0;
          PatchType ptype = ptloop.ptype();
          int ncp;
          if (dim == 0) {
            ncp = 1;
          }
          else if (dim == 1) {
            ncp = deg1 + 1;
          }
          else if (dim == 2 && ptype == BEZIER_TRIANGLE) {
            ncp = (deg1 + 1) * (deg1 + 2) / 2;
          }
          else {
            deg2 = ptloop.degree2();
            ncp = (deg1 + 1) * (deg2 + 1);
          }
          cp.resize(0);
          for (int l = 0; l < ncp; ++l)
            cp.push_back(ptloop.control_point(l));
          pti.insert_patch(deg1, deg2, ptype, cp);
        }
      }
    }
    return newg;
  }
  
  // ------------------------------------------------------------------
  // Apply a transformation to a simplicial complex -- means applying
  // the transformation to nodal coordinates.

  SimpComplex_Under_Construction apply_transformation(const SimpComplex& sc,
    const Matrix& m) {

    int di = sc.embedded_dim();
    int gdim = sc.gdim();
    if (di != m.numcols() - 1)
      throw_error("Wrong size matrix in apply_transformation");
    int numv = sc.num_nodes();
    int newdi = m.numrows();
    vector<Real> srcpt(di + 1);
    vector<Real> destpt(newdi);
    vector<Real> paramcoord;

    vector<int> levsize(gdim + 1);
    {
      for (int fdim = 0; fdim <= gdim; ++fdim)
        levsize[fdim] = sc.brep_levelsize(fdim);
    }

    SimpComplex_Under_Construction newsc(gdim, newdi, levsize);

    {
      int npv = sc.num_prop_val();
      for (int i = 0; i < npv; ++i)
        newsc.add_prop_val(sc.prop(i), sc.val(i));
    }

    for (SimpComplex::VertexOrdinalIndex vnumo = 0; vnumo < numv; ++vnumo) {
      SimpComplex::VertexGlobalIndex vnum = sc.ordinal_to_global(vnumo);
      for (int j = 0; j < di; ++j)
        srcpt[j] = sc.real_coord_o(vnumo,j);
      srcpt[di] = -1.0;
      matrix_vector_mult(m, srcpt, destpt);
      newsc.add_vertex(destpt, vnum);
    }

    {
      for (int fdim = 0; fdim <= gdim; ++fdim) {
        vector<Real> param(fdim);
        vector<SimpComplex::VertexGlobalIndex> vlist(fdim + 1);
        for (Brep::FaceIndex faceind = 0; faceind < levsize[fdim]; ++faceind) {
          Brep::Face_Spec fspec(fdim, faceind);
          int nnf = sc.num_node_on_face(fspec);
          for (int seqno = 0; seqno < nnf; ++seqno) {
            SimpComplex::VertexGlobalIndex vnum = sc.node_on_face(fspec, seqno);
            Brep::PatchIndex patchind = sc.patchind_on_face(fspec, seqno);
            for (int j = 0; j < fdim; ++j)
              param[j] = sc.param_coord_on_face(fspec, seqno, j);
            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 j = 0; j < fdim + 1; ++j)
              vlist[j] = sc.node_of_meshface_on_face(fspec, seqno2, j);
            newsc.add_simplex_face(vlist, fspec);
          }
        }
      }
    }
    return newsc;
  }
  
  
  // standard driver routine.

  void worker(const QMG::FrontEnd::ArgValType& argvalhandle,
    QMG::FrontEnd::ReturnValType& returnvalhandle,
    QMG::FrontEnd::Interpreter& interp) {

    using namespace QMG;
    using namespace QMG::FrontEnd;



    argvalhandle.verify_nargin(2, 2, GM_ROUTINE_NAME_Q);
    returnvalhandle.verify_nargout(1, 1, GM_ROUTINE_NAME_Q);

    
    Object_Type_Code code = argvalhandle.determine_argument_type(1);

    if (code == BREP) {
      Brep_Under_Construction newg = 
        apply_transformation(argvalhandle.get_brep(1), 
        argvalhandle.get_matrix(0));
      returnvalhandle.put_brep(0, newg);
    }
    else if (code == SIMPCOMPLEX) {
      SimpComplex_Under_Construction newsc =
        apply_transformation(argvalhandle.get_simpcomp(1), 
        argvalhandle.get_matrix(0));
      returnvalhandle.put_simpcomp(0, newsc);
    }
    else {
      throw_error("Argument 1 to gmapply neither brep nor simpcomplex");
    }
    return;
  }
}


#include "gm_entrypoint.cpp"