C++ PatchData::num_elements方法代码示例

void NonGradient::printPatch(const PatchData &pd, MsqError &err)
{
  if( mNonGradDebug==0 )
  {
    return;
  }
  const size_t numNode = pd.num_nodes();   //27,  27 what?
  MSQ_PRINT(3)("Number of Vertices: %d\n",(int)pd.num_nodes());
  const size_t numVert = pd.num_free_vertices(); // 1
  MSQ_PRINT(3)("Num Free = %d\n",(int)pd.num_free_vertices());
  const size_t numSlaveVert = pd.num_slave_vertices(); //0
  const size_t numCoin = pd.num_corners(); // 64
  const MsqVertex* coord = pd.get_vertex_array(err);
  MSQ_PRINT(3)("Number of Vertices: %d\n",(int)pd.num_nodes());

  std::cout << "Patch " << numNode << "  " << numVert << "  " << numSlaveVert << "  " << numCoin << std::endl;
  MSQ_PRINT(3)("");
  std::cout << "Coordinate ";
  std::cout << "         " << std::endl;
  for( size_t index = 0; index < numVert; index++ )
  {
    std::cout << coord[index][0] << "  " << coord[index][1] << "  " << coord[index][2] << std::endl;
  }
  //const size_t numElt = pd.num_elements();
  if( mNonGradDebug >= 3 ) 
  {      
         std::cout << "Number of Elements: " << pd.num_elements() << std::endl;
  }      
  MSQ_PRINT(3)("Number of Elements: %d\n",(int)pd.num_elements());
}

   void test_elements()
   {
      MsqPrintError err(cout);
     /* Adds TSTT mesh to a MeshSet. */
      MeshSet mesh_set;
      TSTTMesh tstt_mesh;
      tstt_mesh.set_mesh(tri10);
      mesh_set.add_mesh(&tstt_mesh, err); CPPUNIT_ASSERT(!err);

      /* Retrieves a global patch */
      PatchData pd;
      PatchDataParameters pd_params;
      pd_params.set_patch_type(PatchData::ELEMENTS_ON_VERTEX_PATCH, err, 1, 0);
      
      mesh_set.get_next_patch(pd, pd_params, err); CPPUNIT_ASSERT(!err);

      int free_vtx = pd.num_free_vertices(err); CPPUNIT_ASSERT(!err);
      std::cout << "nb of free vertices: " << free_vtx << std::endl;
      CPPUNIT_ASSERT( free_vtx == 1 );
      
      element_array =  pd.get_element_array(err); CPPUNIT_ASSERT(!err);
      num_elements = pd.num_elements();
      CPPUNIT_ASSERT( num_elements == 6 );

      //      for (int i=0; i<num_elements; ++i) {
      //         std::cout << element_array[i];
      //      }
      
      vtx_array = pd.get_vertex_array(err); CPPUNIT_ASSERT(!err);
      num_vertices = pd.num_vertices();
      CPPUNIT_ASSERT( num_vertices == 7 );

      //      for (int i=0; i<num_vertices; ++i) {
      //         std::cout << vtx_array[i];
      //      }

      CPPUNIT_ASSERT( tri_check_validity() == 1 );
      
      mesh_set.get_next_patch(pd, pd_params, err); CPPUNIT_ASSERT(!err);

      element_array =  pd.get_element_array(err); CPPUNIT_ASSERT(!err);
      num_elements = pd.num_elements();
      CPPUNIT_ASSERT( num_elements == 6 );

      //      for (int i=0; i<num_elements; ++i) {
      //         std::cout << element_array[i];
      //      }
      
      vtx_array = pd.get_vertex_array(err); CPPUNIT_ASSERT(!err);
      num_vertices = pd.num_vertices();
      CPPUNIT_ASSERT( num_vertices == 7 );

      //      for (int i=0; i<num_vertices; ++i) {
      //         std::cout << vtx_array[i];
      //      }

      CPPUNIT_ASSERT( tri_check_validity() == 1 );
   }

void NonSmoothDescent::optimize_vertex_positions(PatchData &pd, 
                                                MsqError &err)
{
  MSQ_FUNCTION_TIMER( "NonSmoothDescent" );

  //  cout << "- Executing NonSmoothDescent::optimize_node_positions()\n";
  /* perform the min max smoothing algorithm */
  MSQ_PRINT(2)("\nInitializing the patch iteration\n");

  MSQ_PRINT(3)("Number of Vertices: %d\n",(int)pd.num_nodes());
  MSQ_PRINT(3)("Number of Elements: %d\n",(int)pd.num_elements());
    //Michael: Note: is this a reliable way to get the dimension?
  mDimension = pd.get_mesh()->get_geometric_dimension(err); MSQ_ERRRTN(err);
  MSQ_PRINT(3)("Spatial Dimension: %d\n",mDimension);

  MSQ_PRINT(3)("Num Free = %d\n",(int)pd.num_free_vertices());

  MsqFreeVertexIndexIterator free_iter(pd, err); MSQ_ERRRTN(err);
  free_iter.reset();
  free_iter.next(); 
  freeVertexIndex = free_iter.value();
  MSQ_PRINT(3)("Free Vertex Index = %d\n",freeVertexIndex);

  // TODO - need to switch to validity via metric evaluations should
  // be associated with the compute_function somehow
  /* check for an invalid mesh; if it's invalid return and ask the user 
     to use untangle */
  if (this->validity_check(pd,err)!=1) {
      MSQ_PRINT(1)("ERROR: Invalid mesh\n");
      MSQ_SETERR(err)("Invalid Mesh: Use untangle to create a valid "
                      "triangulation", MsqError::INVALID_MESH);
      return;
  }

  /* assumes one function value per element */
  // TODO - need to include vertex metrics 
  numFunctionValues = pd.num_elements();

  /* initialize the optimization data up to numFunctionValues */
  this->init_opt(err); MSQ_ERRRTN(err);
  this->init_max_step_length(pd,err);  MSQ_ERRRTN(err);
  MSQ_PRINT(3)("Done initializing optimization\n");

  /* compute the initial function values */
  //TODO this should return a bool with the validity
  this->compute_function(&pd, originalFunction, err);  MSQ_ERRRTN(err);
 
  // find the initial active set
  this->find_active_set(originalFunction, mActive, err);   MSQ_ERRRTN(err);

  this->minmax_opt(pd,err);  MSQ_ERRRTN(err);
}

void CachingTargetCalculator::notify_sub_patch( PatchData& ,
                                                CachedTargetData& data,
                                                PatchData& subpatch, 
                                                const size_t* ,
                                                const size_t* element_map,
                                                MsqError& err )
{
    // If no cached data for this patch, just return
  if (data.has_data())
    return;
  
    // Create a new cached data object on the subpatch
  CachedTargetData& sub_data = get_data( subpatch );
  sub_data.clear();
  
    // populate the element offset list, and count the total
    // number of cached target matrices.
  sub_data.elementOffsets.resize( subpatch.num_elements() );
  size_t count_2D = 0, count_3D = 0;
  for (size_t i = 0; i < subpatch.num_elements(); ++i) {
    EntityTopology type = subpatch.element_by_index(i).get_element_type();
    size_t& count = (TopologyInfo::dimension( type ) == 2) ? count_2D : count_3D;
    sub_data.elementOffsets[i] = count;
    NodeSet samples = subpatch.get_samples( i );
    count += samples.num_nodes();
  }
  
  const bool orient = have_surface_orient();
  sub_data.targets3D.resize( count_3D );
  if (orient) 
    sub_data.targetsSurface.resize( count_2D );
  else
    sub_data.targets2D.resize( count_2D );

  for (size_t i = 0; i < subpatch.num_elements(); ++i) {
    EntityTopology type = subpatch.element_by_index(i).get_element_type();
    size_t off = sub_data.elementOffsets[i];
    size_t old_off = data.elementOffsets[element_map[i]];
    NodeSet samples = subpatch.get_samples( i );
    size_t count = samples.num_nodes();
   
    if (TopologyInfo::dimension( type ) == 3) 
      memcpy( &sub_data.targets3D[off], &data.targets3D[old_off], count*sizeof(MsqMatrix<3,3>) );
    else if (orient)
      memcpy( &sub_data.targetsSurface[off], &data.targetsSurface[old_off], count*sizeof(MsqMatrix<3,2>) );
    else
      memcpy( &sub_data.targets2D[off], &data.targets2D[old_off], count*sizeof(MsqMatrix<2,2>) );
  }
}

int NonGradient::getPatchDimension(const PatchData &pd, MsqError &err)
{
  const size_t numElt = pd.num_elements();
  unsigned edimMax = 0;  // in case numElt == 0
  for( size_t elementId = 0; elementId < numElt; elementId++)
  {
    const MsqMeshEntity& element = pd.element_by_index( elementId );
    EntityTopology type = element.get_element_type();
    unsigned edim = TopologyInfo::dimension( type );
    if( elementId == 0 ) 
    {
      edimMax = edim;
    } 
    else
    { 
      if(edimMax != edim)
      {
        MSQ_SETERR(err)("A patch has elements of different dimensions", MsqError::INVALID_MESH);
        std::cout << "A patch has elements of different dimensions" << edimMax << "  "  << edim << std::endl;
        if(edimMax < edim)
        {
          edimMax = edim;
        }
      }
    }
  }
  return( edimMax );
}

void PMeanPMetricTest::test_get_element_evaluations()
{
  MsqError err;
  FauxMetric m;
  ElementPMeanP e( 1.0, &m );
  std::vector<size_t> handles;
    // test that handles array contains all elements
  e.get_evaluations( pd, handles, false, err );
  std::sort( handles.begin(), handles.end() );
  CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
  for (unsigned i = 1; i < handles.size(); ++i)
    CPPUNIT_ASSERT_EQUAL( handles[i-1]+1, handles[i] );
    // test that handles array contains all elements
  e.get_evaluations( pd, handles, true, err );
  std::sort( handles.begin(), handles.end() );
  CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
  for (unsigned i = 1; i < handles.size(); ++i)
    CPPUNIT_ASSERT_EQUAL( handles[i-1]+1, handles[i] );
}

size_t TerminationCriterion::count_inverted( PatchData& pd, MsqError& err )
{
  size_t num_elem = pd.num_elements();
  size_t count=0;
  int inverted, samples;
  for (size_t i = 0; i < num_elem; i++) {
    pd.element_by_index(i).check_element_orientation(pd, inverted, samples, err);
    if (inverted)
      ++count;
  }
  return count;
}

size_t SmartLaplacianSmoother::num_inverted( PatchData& pd, MsqError& err )
{
  size_t result = 0;
  int inverted, junk;
  for (size_t i = 0; i < pd.num_elements(); ++i) {
    pd.element_by_index(i).check_element_orientation( pd, inverted, junk, err ); 
    MSQ_ERRZERO(err);
    if (inverted)
      ++result;
  }
  return result;
}

void PatchDataTestNormals::test_get_element_normals_infinite_domain()
{
  MsqPrintError err(cout);
  Vector3D expected_normals[] = { Vector3D(  0,  0, -1),
                                  Vector3D(  0,  0,  1),
                                  Vector3D(  1,  0,  0),
                                  Vector3D(  0,  1,  0),
                                  Vector3D( -1,  0,  0),
                                  Vector3D(  0, -1,  0) };
                                   
  CPPUNIT_ASSERT( unboundedMesh.num_elements() == 6u );
  for (size_t i = 0; i < 6u; ++i)
  {
    Vector3D norm;
    unboundedMesh.get_domain_normal_at_element( i, norm, err );
    CPPUNIT_ASSERT(!err);
    ASSERT_VECTORS_EQUAL( expected_normals[i], norm );
  }    
}

/*! Writes a gnuplot file directly from the MeshSet.
    This means that any mesh imported successfully into Mesquite
    can be outputed in gnuplot format.

    Within gnuplot, use \b plot 'file1.gpt' w l, 'file2.gpt' w l  
    
    This is not geared for performance, since it has to load a global Patch from
    the mesh to write a mesh file. 
*/
void MeshSet::write_gnuplot(const char* out_filebase,
                   Mesquite::MsqError &err)
{
    // Open the file
  string out_filename = out_filebase;
  out_filename += ".gpt";
  ofstream file(out_filename.c_str());
  if (!file)
  {
    MSQ_SETERR(err)(MsqError::FILE_ACCESS);
    return;
  }

    // loads a global patch
  PatchData pd;
  PatchDataParameters pd_params;
  pd_params.set_patch_type(PatchData::GLOBAL_PATCH, err); MSQ_ERRRTN(err);
  pd_params.no_culling_method();
  get_next_patch(pd, pd_params, err); MSQ_ERRRTN(err);
    
    // Write a header
  file << "\n";
  
  for (size_t i=0; i<pd.num_elements(); ++i)
  {
    std::vector<size_t> vtx_indices;
    pd.elementArray[i].get_node_indices(vtx_indices);
    for (size_t j = 0; j < vtx_indices.size(); ++j)
    {
      file << pd.vertexArray[vtx_indices[j]][0] << ' '
           << pd.vertexArray[vtx_indices[j]][1] << ' '
           << pd.vertexArray[vtx_indices[j]][2] << '\n';
    }
      file << pd.vertexArray[vtx_indices[0]][0] << ' '
           << pd.vertexArray[vtx_indices[0]][1] << ' '
           << pd.vertexArray[vtx_indices[0]][2] << '\n';
    file << '\n';
  }
  
    // Close the file
  file.close();
}

/*! The type of targets computed by this function is selected by the constructor of
    the base classes. */
void WTargetCalculator::compute_target_matrices(PatchData &pd, MsqError &err)
{
  MSQ_FUNCTION_TIMER( "WTargetCalculator::compute_target_matrice" );

  size_t num_elements=pd.num_elements();
  
  PatchData ref_pd, *ref_pd_ptr;
  if ( refMesh != 0 ) {
    // If there is a reference mesh, gets a patch ref_pd equivalent to the patch pd of the main mesh.
    PatchDataParameters ref_pd_params(this->get_all_parameters());
    refMesh->get_next_patch(ref_pd, ref_pd_params, err); MSQ_ERRRTN(err);
    // Make sure topology of ref_pd and pd are equal
    assert( num_elements == ref_pd.num_elements() );
    size_t num_vertices=pd.num_vertices();
    assert( num_vertices == ref_pd.num_vertices() );
    ref_pd_ptr = &ref_pd;
  }
  else {
    // the reference patch is the same as the working patch if there is no reference mesh.
    ref_pd_ptr = &pd;
  }
  
  MsqMeshEntity* elems = pd.get_element_array(err); MSQ_ERRRTN(err);
  MsqMeshEntity* elems_ref = ref_pd_ptr->get_element_array(err); MSQ_ERRRTN(err);

  Matrix3D W_guides[MSQ_MAX_NUM_VERT_PER_ENT];
  TargetMatrix matrices[MSQ_MAX_NUM_VERT_PER_ENT];
  
  for (size_t i=0; i<num_elements; ++i) {
    int nve = elems[i].vertex_count();
    assert( nve = elems_ref[i].vertex_count() );

    compute_guide_matrices(guideMatrix, *ref_pd_ptr, i, W_guides, nve, err); MSQ_ERRRTN(err);
    for (int c = 0; c < nve; ++c)
      matrices[c] = W_guides[c];
    pd.targetMatrices.set_element_corner_tags( &pd, i, matrices, err ); MSQ_ERRRTN(err);
  }

  //if ( refMesh != 0 ) delete ref_pd;
}

/* virtual function reimplemented from QualityMetric. No doxygen doc needed. */
bool LPtoPTemplate::compute_analytical_gradient(PatchData &pd,
                                              Vector3D *const &grad,
                                              double &OF_val,
                                              MsqError &err, size_t array_size)
{
  MSQ_FUNCTION_TIMER( "LPtoPTemplate::compute_analytical_gradient" );
  
    //initialize the scaling value
  double scaling_value=1.0;
 
  size_t num_elements=pd.num_elements();
  size_t num_vertices=pd.num_vertices();
  if( num_vertices!=array_size && array_size>0)
  {
    MSQ_SETERR(err)("Incorrect array size.", MsqError::INVALID_ARG);
    return false;
  }
  
  MsqMeshEntity* elems=pd.get_element_array(err); MSQ_ERRZERO(err);
  MsqVertex* vertices=pd.get_vertex_array(err); MSQ_ERRZERO(err);
  bool qm_bool=true;
  double QM_val;
  OF_val = 0.;
  size_t i;
  int p1;
  
  //Set currentQM to be quality metric (possibly composite) associated with the objective function
  QualityMetric* currentQM = get_quality_metric();
  if(currentQM==NULL) {
    MSQ_SETERR(err)("LPtoPTemplate has NULL QualityMetric pointer.",MsqError::INVALID_STATE);
    return false;
  }
  enum QualityMetric::MetricType qm_type=currentQM->get_metric_type();
  
  if (qm_type!=QualityMetric::ELEMENT_BASED &&
      qm_type!=QualityMetric::VERTEX_BASED) {
    MSQ_SETERR(err)("Make sure MetricType is initialised"
                    "in concrete QualityMetric constructor.",
                    MsqError::INVALID_STATE);
    return false;
  }


  // zeros out objective function gradient
  for (i=0; i<num_vertices; ++i)
    grad[i] =0;
  
  // Computes objective function gradient for an element based metric
  if(qm_type==QualityMetric::ELEMENT_BASED){
      //if scaling, divid by num_elements
    if(dividingByN){
      if(num_elements<=0) {
        MSQ_SETERR(err)("The number of elements should not be zero.",MsqError::INVALID_MESH);
        return false;
      }
      
      scaling_value/=num_elements;
    }
    size_t e, ve;
    size_t nfve; // num free vtx in element
    size_t nve; // num vtx in element
    MsqVertex* ele_free_vtces[MSQ_MAX_NUM_VERT_PER_ENT];
    const size_t *ele_vtces_ind;

    // loops over all elements.
    for (e=0; e<num_elements; ++e) {
      // stores the pointers to the free vertices within the element
      // (using pointer arithmetic).
      nfve = 0;
      nve = elems[e].vertex_count();
      ele_vtces_ind = elems[e].get_vertex_index_array();
      for (ve=0; ve<nve; ++ve) {
        if (vertices[ele_vtces_ind[ve]].is_free_vertex()) {
          ele_free_vtces[nfve] = vertices + ele_vtces_ind[ve];
          ++nfve;
        }
      }

      // Computes q and grad(q)
      Vector3D grad_vec[MSQ_MAX_NUM_VERT_PER_ENT];
      qm_bool = currentQM->compute_element_gradient(
                                     pd, &elems[e],
                                     ele_free_vtces,
                                     grad_vec, nfve, QM_val, err);
      if(MSQ_CHKERR(err) || !qm_bool) return false;

      // computes p*|Q(e)|^{p-1}
      QM_val = fabs(QM_val);
      double QM_pow=1.0;
      double factor;
      if (pVal==1) factor=1;
      else {
        QM_pow=QM_val;
        for (p1=1; p1<pVal-1; ++p1)
          QM_pow*=QM_val;
        factor = QM_pow * pVal;
      }
        //this scales the gradient
      factor *= (scaling_value * get_negate_flag());
//.........这里部分代码省略.........

void MsqMeshEntityTest::test_check_element_orientation( EntityTopology type, 
                                                        int nodes )
{
    // get an ideal element
  MsqError err;
  PatchData pd;
  create_ideal_element_patch( pd, type, nodes, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT_EQUAL( (size_t)1, pd.num_elements() );
  CPPUNIT_ASSERT_EQUAL( (size_t)nodes, pd.num_nodes() );
  MsqMeshEntity& elem = pd.element_by_index(0);
  CPPUNIT_ASSERT_EQUAL( (size_t)nodes, elem.node_count() );
  CPPUNIT_ASSERT_EQUAL( type, elem.get_element_type() );
  const size_t* conn = elem.get_vertex_index_array();
  
    // test that ideal element is not reported as inverted
  int inverted, tested;
  elem.check_element_orientation( pd, inverted, tested, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT_EQUAL( 0, inverted );
  CPPUNIT_ASSERT( tested > 0 );

  bool mids[4] = {false};
  TopologyInfo::higher_order( type, nodes, mids[1], mids[2], mids[3], err );
  MSQ_ERRRTN(err);
  
    // invert element at each vertex and test
  Vector3D centroid;
  elem.get_centroid( centroid, pd, err );
  ASSERT_NO_ERROR(err);
  for (int i = 0; i < nodes; ++i) {
    unsigned dim, num;
    TopologyInfo::side_from_higher_order( type, nodes, i, dim, num, err );
    ASSERT_NO_ERROR(err);
    const Vector3D old_pos = pd.vertex_by_index( conn[i] );
    Vector3D new_pos = old_pos;
    if (dim == TopologyInfo::dimension(type)) { 
      // move mid-element node 3/4 of the way to corner 0
      new_pos += 3*pd.vertex_by_index( conn[0] );
      new_pos *= 0.25;
    }
    else if (dim == 0) { // if a corner vertex
      if (type == TRIANGLE || type == TETRAHEDRON) {
        // move tri/tet vertex past opposite side of element
        new_pos += 2*(centroid - old_pos);
      }
      else if (mids[1]) {
        // if have mid-edge nodes move 3/4 of the way to center vertex
        new_pos += 3*centroid;
        new_pos *= 0.25;
      }
      else {
        // move vertex past centroid
        new_pos += 1.5*(centroid - old_pos);
      }
    }
    else {
      // otherwise move vertex past centroid
      new_pos += 2.5*(centroid - old_pos);
    }
    
    pd.set_vertex_coordinates( new_pos, conn[i], err );
    ASSERT_NO_ERROR(err);
    
      // test that element is inverted
    inverted = tested = 0;
    elem.check_element_orientation( pd, inverted, tested, err );
    ASSERT_NO_ERROR(err);
    std::ostringstream str;
    str << TopologyInfo::short_name(type) << nodes
        << " Vertex " << i 
        << " (Dimension " << dim
        << " Index " << num << ")";
    CppUnit::Message m( "MsqMeshEntity failed to detect inverted element" );
    m.addDetail( str.str() );
    ASSERT_MESSAGE( m, inverted > 0 );
    
      // move vertex back to ideal position
    pd.set_vertex_coordinates( old_pos, conn[i], err );
    ASSERT_NO_ERROR(err);
  }
}

void FauxMetric::get_evaluations( PatchData& pd, std::vector<size_t>& h, bool free, MsqError& err )
{
  h.clear();
  for (size_t i = 0; i < pd.num_elements(); ++i)
    get_element_evaluations( pd, i, h, err );
}

/*! \brief creates a sparse structure for a Hessian, based on the
  connectivity information contained in the PatchData.
  Only the upper triangular part of the Hessian is stored. */
void MsqHessian::initialize(PatchData &pd, MsqError &err)
{
  MSQ_FUNCTION_TIMER( "MsqHession::initialize" );
  delete[] mEntries;
  delete[] mRowStart;
  delete[] mColIndex;
  
  size_t num_vertices = pd.num_free_vertices();
  size_t num_elements = pd.num_elements();
  size_t const * vtx_list;
  size_t e, r, rs, re, c, cs, ce, nz, nnz, nve, i, j;
  MsqMeshEntity* patchElemArray = pd.get_element_array(err); MSQ_CHKERR(err);

  if (num_vertices == 0) {
    MSQ_SETERR( err )( "No vertices in PatchData", MsqError::INVALID_ARG);
    return;
  }

  mSize = num_vertices;

  // Calculate the offsets for a CSC representation of the accumulation
  // pattern.

  size_t* col_start = new size_t[num_vertices + 1];
  //mAccumElemStart = new size_t[num_elements+1];
  //mAccumElemStart[0] = 0;
  
  for (i = 0; i < num_vertices; ++i) {
    col_start[i] = 0;
  }

  for (e = 0; e < num_elements; ++e) {
    nve = patchElemArray[e].node_count();
    vtx_list = patchElemArray[e].get_vertex_index_array();
    int nfe = 0;
    
    for (i = 0; i < nve; ++i) {
      r = vtx_list[i];
      if (r < num_vertices)
        ++nfe;
      
      for (j = i; j < nve; ++j) {
        c = vtx_list[j];

        if (r <= c) {
          if (c < num_vertices)
            ++col_start[c];
        }
        else {
          if (r < num_vertices)
            ++col_start[r];
        }
      }
    }
    //mAccumElemStart[e+1] = mAccumElemStart[e] + (nfe+1)*nfe/2;
  }

  nz = 0;
  for (i = 0; i < num_vertices; ++i) {
    j = col_start[i];
    col_start[i] = nz;
    nz += j;
  }
  col_start[i] = nz;

  // Finished putting matrix into CSC representation

  int* row_instr = new int[5*nz];
  size_t* row_index = new size_t[nz];

  nz = 0;
  for (e = 0; e < num_elements; ++e) {
    nve = patchElemArray[e].node_count();
    vtx_list = patchElemArray[e].get_vertex_index_array();

    for (i = 0; i < nve; ++i) {
      r = vtx_list[i];

      for (j = i; j < nve; ++j) {
        c = vtx_list[j];

        if (r <= c) {
          if (c < num_vertices) {
            row_index[col_start[c]] = r;
            row_instr[col_start[c]] = nz++;
            ++col_start[c];
          }
        }
        else {
          if (r < num_vertices) {
            row_index[col_start[r]] = c;
            //can't use -nz, but can negate row_instr[col_start[r]]
            row_instr[col_start[r]] = nz++;
            row_instr[col_start[r]] = -row_instr[col_start[r]];
            ++col_start[r];
          }
        }
//.........这里部分代码省略.........