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

   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 test_unsigned_area()
    {
      MsqPrintError err(cout);
      MsqMeshEntity* tri = oneTriPatch.get_element_array(err);
      CPPUNIT_ASSERT(!err);
      CPPUNIT_ASSERT(fabs(tri->compute_unsigned_area(oneTriPatch,err)
                          -(sqrt(3.0)/4.0)) < tolEps);
      MsqMeshEntity* quad = oneQuadPatch.get_element_array(err);
      CPPUNIT_ASSERT(!err);
      CPPUNIT_ASSERT(fabs(quad->compute_unsigned_area(oneQuadPatch,err)
                          -1.0) < tolEps);
    }

  //! test the centroid of the first element in the Patch
  void test_centroid(PatchData &pd, Vector3D &correct)
  {
    MsqPrintError err(cout);
    double eps = 1e-6;
    Vector3D centroid;

    MsqMeshEntity* elem = pd.get_element_array(err); CPPUNIT_ASSERT(!err);
    elem->get_centroid(centroid, pd, err); CPPUNIT_ASSERT(!err);

//     cout << "centroid: "<< centroid <<endl; 
//     cout << "correct: "<< correct <<endl; 

    for (int i=0; i<3; ++i)
      CPPUNIT_ASSERT_DOUBLES_EQUAL( centroid[i], correct[i], eps);
  }

/*! 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;
}

bool VertexConditionNumberQualityMetric::evaluate_with_indices( 
                                              PatchData& pd,
                                              size_t this_vert,
                                              double& value,
                                              std::vector<size_t>& indices,
                                              MsqError& err )
{
  bool rval = evaluate( pd, this_vert, value, err ); MSQ_ERRFALSE(err);
  
  indices.clear();

  MsqMeshEntity* elems = pd.get_element_array(err);
  size_t num_elems;
  const size_t *v_to_e_array = pd.get_vertex_element_adjacencies( this_vert, num_elems, err );
  MSQ_ERRZERO(err);  
  
    //vector to hold the other verts which form a corner.
  vector<size_t> other_vertices;
  other_vertices.reserve(4);
  size_t i=0;

  //loop over the elements attached to this vertex
  for(i=0;i<num_elems;++i){
      //get the vertices connected to this vertex for this element
    elems[v_to_e_array[i]].get_connected_vertices(this_vert,
                                                  other_vertices,
                                                  err);  MSQ_ERRZERO(err);
    for (unsigned j = 0; j < other_vertices.size(); ++j) {
      if (other_vertices[j] < pd.num_free_vertices())
        indices.push_back(other_vertices[j]);
    }
  }
  
  std::sort( indices.begin(), indices.end() );
  indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
  if (this_vert < pd.num_free_vertices())
    indices.push_back( this_vert );
  return rval;
}

/*! \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];
          }
        }
//.........这里部分代码省略.........

/*!For the given vertex, vert, with connected elements, e_i for i=1...K,
  the LocalSizeQualityMetric computes the corner volumes (or areas) of
  each e_i at the corner defined by vert.  The corner volume is defined
  as the volume of the tet defined by the edges of an element which contain
  the common vertex, vert.  That volume is then diveded by the average corner
  volume of all the element corners connected to this vertex.  For
  vertices attached to pyramid elements, this metric is undefined.
*/
bool LocalSizeQualityMetric::evaluate( PatchData &pd, size_t this_vert,
                                       double &fval, MsqError &err )
{
  fval=0.0;
    //get the element array
  MsqMeshEntity* elems = pd.get_element_array(err);  MSQ_ERRZERO(err);
    //get the vertex to element array and the offset array
  //const size_t* elem_offset = pd.get_vertex_to_elem_offset(err);  MSQ_ERRZERO(err);
  //const size_t* v_to_e_array = pd.get_vertex_to_elem_array(err);  MSQ_ERRZERO(err);
    //find the offset for this vertex
  //size_t this_offset = elem_offset[this_vert];
    //get the number of elements attached to this vertex (given by the
    //first entry in the vertex to element array)
  //size_t num_elems = v_to_e_array[this_offset];
    //PRINT_INFO("\nIN LOCAL SIZE CPP, num_elements = %i",num_elems);
  size_t num_elems;
  const size_t *v_to_e_array = pd.get_vertex_element_adjacencies( this_vert, num_elems, err ); 
  MSQ_ERRZERO(err);
  
  if(num_elems <= 0){
    return true;
  }
  
    //create an array to store the local metric values before averaging
    //Can we remove this dynamic allocatio?
  double* met_vals = new double[num_elems];
    //vector to hold the other verts which form a corner.
  std::vector<size_t> other_vertices;
  other_vertices.reserve(4);
  double total_val=0.0;
  size_t i=0;
    //loop over the elements attached to this vertex
  for(i=0;i<num_elems;++i){
      //get the vertices which (with this_vert) form the corner of
      //the ith element.
    elems[v_to_e_array[i]].get_connected_vertices(this_vert,
                                                              other_vertices,
                                                              err);  MSQ_ERRZERO(err);
      ////PRINT_INFO("\nINSIDE LOCAL SIZE CPP other_vertices size = %i",other_vertices.size());
    
    switch(other_vertices.size()){
        //if a surface element, compute the corner area
      case 2:
        met_vals[i] = compute_corner_area(pd, this_vert, other_vertices[0],
                                          other_vertices[1], err);  MSQ_ERRZERO(err);
        break;
          //if a volume element, compute the corner volume 
      case 3:
        met_vals[i] = compute_corner_volume(pd, this_vert, other_vertices[0],
                                            other_vertices[1],
                                            other_vertices[2], err);  MSQ_ERRZERO(err);
        break;
      default:
          //otherwise, there is was an error.  Either the wrong number
          //of vertices were returned fom get_connected_vertices or
          //the element does not have the correct number of edges
          //connected to this vertex (possibly a pyramid element).
        met_vals[i]=0.0;
        MSQ_SETERR(err)("Incorrect number of vertices returned from "
                        "get_connected_vertices.", MsqError::UNSUPPORTED_ELEMENT);
        return false;
    };
      //keep track of total so that we can compute the linear average
    total_val+=met_vals[i];
    //PRINT_INFO("\nIN LOCAL SIZE CPP, total_val = %f, i = %i",total_val,i);
      //clear the vector of other_vertices for re-use.
    other_vertices.clear();
    //PRINT_INFO("\nIN LOCAL SIZE CPP, after clean size = %f",other_vertices.size());
    
  }
    //calculate the linear average... num_elems is non-zero here.
  total_val /= (double) num_elems;
  //PRINT_INFO("\nIN LOCAL SIZE CPP, average = %f",total_val);
    //if the average is non-zero
    //divide each entry by the linear average
  if(total_val!=0){
    for(i=0;i<num_elems;++i){
      met_vals[i]/=total_val;
    }
      //calculate fval by averaging the corner values
    fval = average_metrics(met_vals, num_elems, err);  MSQ_ERRZERO(err);
    //PRINT_INFO("\nIN LOCAL SIZE CPP, inside if statement");
  }
  //PRINT_INFO("\nIN LOCAL SIZE CPP, fval = %f",fval);
    //clean up the dynamically allocated array
  delete []met_vals;
    //always return true... the vertex position is never invalid
  return true;
  
}

bool MeshSet::get_next_global_patch( PatchData& pd,
                                     PatchDataParameters& pd_params,
                                     MsqError& err )
{
    // We only support global patches for a single Mesh
  if (meshSet.size() != 1)
  {
    MSQ_SETERR(err)( 
       "Global patches only supported for single-Mesh MeshSets.",
       MsqError::NOT_IMPLEMENTED );
    return false;
  }

  pd.mType = PatchData::GLOBAL_PATCH;
  pd.domainHint = NO_DOMAIN_HINT;
  if (mDomain)
    pd.domainHint = mDomain->hint();

    // for a global patch, we always reset to start of the mesh.
  reset(err);
  if (MSQ_CHKERR(err)) return false;

  size_t i;

    // Get sizes for mesh data
  size_t num_verts, num_elems, num_uses;
  (*currentMesh)->get_all_sizes( num_verts, num_elems, num_uses, err ); MSQ_ERRZERO(err);

    // Get handles and connectivity
  pd.vertexHandlesArray.resize( num_verts );
  pd.elementHandlesArray.resize( num_elems );
  pd.elemConnectivityArray.resize( num_uses );
  msq_std::vector<size_t> offsets(num_elems+1);
  (*currentMesh)->get_all_mesh( &pd.vertexHandlesArray[0], num_verts,
                                &pd.elementHandlesArray[0], num_elems,
                                &offsets[0], offsets.size(),
                                &pd.elemConnectivityArray[0], 
                                pd.elemConnectivityArray.size(),
                                err ); MSQ_ERRZERO(err);
  
    // Get element topologies
  pd.elementArray.resize( num_elems );
  msq_std::vector<EntityTopology> elem_topologies(num_elems);
  (*currentMesh)->elements_get_topologies( &pd.elementHandlesArray[0], 
                                           &elem_topologies[0], 
                                           num_elems, err );MSQ_ERRZERO(err);

    // Put them into the patch
  MsqMeshEntity* pd_elem_array = pd.get_element_array(err);MSQ_ERRZERO(err);
  for (i = 0; i < num_elems; ++i)
    pd_elem_array[i].set_element_type( elem_topologies[i] );
    
    // Complete connectivity data in patch
  pd.initialize_data( &offsets[0], err ); MSQ_ERRZERO(err);


    // Get vertex coordinates
  pd.vertexArray.resize( num_verts );
  MsqVertex* pd_vert_array = pd.get_vertex_array(err);MSQ_ERRZERO(err);
  (*currentMesh)->vertices_get_coordinates(&pd.vertexHandlesArray[0],
                                           pd_vert_array,
                                           num_verts,
                                           err); MSQ_ERRZERO(err);

    // Get vertex boundary flag
  if (vertArraySize < num_verts)
  {
    delete [] vertexOnBoundary;
    vertArraySize = num_verts;
    vertexOnBoundary = new bool[vertArraySize];
  }
  (*currentMesh)->vertices_are_on_boundary( &pd.vertexHandlesArray[0],
                                            vertexOnBoundary,
                                            num_verts,
                                            err );MSQ_ERRZERO(err);

  for (i = 0; i < num_verts; i++)
  {
      // Get its flags
    /*(*currentMesh)->vertex_get_byte(vertArray[i],
                                    &(pd_vert_array[i].vertexBitFlags),
                                    err); MSQ_CHKERR(err);*/
      // Set its hard-fixed flag
    if (/*(*currentMesh)->vertex_is_fixed(vertArray[i], err) ||*/
       vertexOnBoundary[i])
    {
      pd_vert_array[i].vertexBitFlags |= MsqVertex::MSQ_HARD_FIXED;
    }
    else
    {
      pd_vert_array[i].vertexBitFlags &= ~(MsqVertex::MSQ_HARD_FIXED);
    }
  }

  return true;
}

//.........这里部分代码省略.........
      // Get the number of elements in this vertex
  size_t num_elems =
    (*currentMesh)->vertex_get_attached_element_count(vertex, err);
  if (MSQ_CHKERR(err)) return false;
  pd.elementHandlesArray.resize( num_elems );
 
    // Get the elements attached to this vertex
  if (elemArraySize < num_elems)
  {
    delete [] elemTopologies;
    elemTopologies = new EntityTopology[num_elems];
    elemArraySize = num_elems;
  }

  (*currentMesh)->vertex_get_attached_elements(vertex,
                                               &pd.elementHandlesArray[0],
                                               num_elems, err);
  if (MSQ_CHKERR(err)) return false;

    // Get the topologies of those elements
  (*currentMesh)->elements_get_topologies(&pd.elementHandlesArray[0],
                                          elemTopologies,
                                          num_elems, err);
  if (MSQ_CHKERR(err)) return false;

    // Figure out how many vertices we need to allocate
  //size_t num_vert_uses = 1;
  //size_t i;
  //for (i = 0; i < num_elems; ++i)
  //  num_vert_uses += vertices_in_topology(elemTopologies[i]);
  size_t num_vert_uses = (*currentMesh)->
    get_vertex_use_count( &pd.elementHandlesArray[0], num_elems, err ); 
  MSQ_ERRZERO(err);

    // All elems share at least 1 vertex (the center vertex).  The
    // center vertex is used 1 time, but it was counted num_elems times.
  size_t num_verts = num_vert_uses - num_elems + 1;
  pd.vertexHandlesArray.resize( num_verts );
  pd.elementArray.resize( num_elems );
  pd.elemConnectivityArray.resize( num_vert_uses );

    // Get the vertices attached to those elements
  
  if (csrOffsetsSize < num_elems + 1)
  {
    delete [] csrOffsets;
    csrOffsets = new size_t[num_elems + 1];
    csrOffsetsSize = num_elems + 1;
  }
  (*currentMesh)->elements_get_attached_vertices(&pd.elementHandlesArray[0],
                                                 num_elems,
                                                 &pd.vertexHandlesArray[0],
                                                 num_verts,
                                                 &pd.elemConnectivityArray[0],
                                                 num_vert_uses,
                                                 csrOffsets,
                                                 err); 
  if (MSQ_CHKERR(err)) return false;
  
    // Update with actual vertex count
  pd.vertexHandlesArray.resize( num_verts );

    // Put the elements into the PatchData
  MsqMeshEntity* pd_elem_array = pd.get_element_array(err);
  for (i = 0; i < num_elems; ++i)
    pd_elem_array[i].set_element_type( elemTopologies[i] );
  pd.initialize_data( csrOffsets, err ); MSQ_ERRZERO(err);

    // Get the coordinates of the vertices and its flags.
  pd.vertexArray.resize( num_verts );
  MsqVertex* pd_vert_array = pd.get_vertex_array(err);
    //get the coordinates
  (*currentMesh)->vertices_get_coordinates(&pd.vertexHandlesArray[0],
                                           pd_vert_array,
                                           num_verts,
                                           err);
  if (MSQ_CHKERR(err)) return false;
  for (i = 0; i < num_verts; i++)
  {

      // If it's not the center vertex, mark it as hard fixed
    if (pd.vertexHandlesArray[i] != vertex)
    {
        // Get its flags
      (*currentMesh)->vertex_get_byte(pd.vertexHandlesArray[i],
                                      &(pd_vert_array[i].vertexBitFlags),
                                      err);
      if (MSQ_CHKERR(err)) return false;
      pd_vert_array[i].vertexBitFlags |= MsqVertex::MSQ_HARD_FIXED;
    }
      //else it is the center vertex.  We therefore already have
      //the fixed flag stored center_fixed_byte.  The hard fixed
      //flag has already been removed (when flag was retreived).
    else{
      pd_vert_array[i].vertexBitFlags = (center_fixed_byte);
    }
  }

  return true;
}

  bool I_DFT_NoBarrierSmoother::i_dft_no_barrier_smooth_mesh(
    PatchData &pd, size_t free_index, MsqError &err)
  {
      //get vertex and element arrays
    MsqVertex* verts=pd.get_vertex_array(err); MSQ_ERRFALSE(err);
    MsqMeshEntity* elems=pd.get_element_array(err);MSQ_ERRFALSE(err);
    size_t num_elements = pd.num_elements();
      //matrix to store a map for the given corner
      // (element/vertex combination) that maps the
      //vertex to 0 and the other vertices to 1,2,3
    int local_matrix_map[MSQ_MAX_NUM_VERT_PER_ENT];
    size_t local_matrix_map_length=MSQ_MAX_NUM_VERT_PER_ENT;
    size_t local_matrix_map_used=0;
    Vector3D new_position;
      //will be used as the denominator
    double scale_value = 0.0;
    size_t i, j, k;
      //Z and Z^(t)
    Matrix3D z_mat_transpose;
    Matrix3D z_mat[4];
      //array of Y's... one for each vertex in this corner.
    Matrix3D y_mat[4];
      //c_(j,k) for the vertices
    double c_scalars[4];
      //zeta for this corner
    Vector3D zeta_vector;
    vector<size_t> elems_verts;
      //loop over the two or three dimensions
    for(i=0;i<num_elements;++i) {
        //actually get the map for this corner...
      local_matrix_map_used = elems[i].get_local_matrix_map_about_vertex(
        pd, &verts[free_index],local_matrix_map_length, local_matrix_map,err);
      MSQ_ERRFALSE(err);
      
        //get a vector of the vertex indices for this element
      elems[i].get_vertex_indices(elems_verts);
        //get the W array for this elements
      int elem_idx = pd.get_element_index(&elems[i]);
      const TargetMatrix* W = pd.targetMatrices.get_element_corner_tags(
        &pd, elem_idx, err );      MSQ_ERRFALSE(err);

        //initial c_scalars and zeta_vector to 0.0 and calculate Y for
        //each vertex in the this corner.
      zeta_vector.set( 0.0, 0.0, 0.0);
      for(j=0;j<local_matrix_map_used;++j){
        c_scalars[j]=0.0;
        if(local_matrix_map[j] < 0){
          MSQ_SETERR(err)("Invalid index returned from MsqMeshEntity.\n",
                          MsqError::INVALID_STATE);
          return false;
        }
        inv(z_mat[j],W[local_matrix_map[j]]);
          //inv(z_mat[j], I);
        z_mat_transpose=transpose(z_mat[j]);
        matmult( y_mat[j], z_mat[j], z_mat_transpose);
          //std::cout<<"\n j "<<j<<"\n";
          //std::cout<<"\n"<<y_mat[j]<<"\n";
      }
        //0 case (c_scalars[0])
      for(j=0;j<3;++j){
        for(k=0;k<3;++k){
          c_scalars[0] += ((y_mat[0])[j][k]);
        }
        c_scalars[0] += ((y_mat[j+1])[j][j]);
      }
        //1, 2, and if need 3
      for(j=0;j<3;++j){
        for(k=0;k<3;++k){
          c_scalars[j+1] -= (((y_mat[0])[j][k]) + ((y_mat[j+1])[k][j]));
        }
        c_scalars[j+1] += ((y_mat[((j+2)%3)+1])[((j+1)%3)][((j+2)%3)]);
        c_scalars[j+1] += ((y_mat[((j+1)%3)+1])[((j+2)%3)][((j+1)%3)]);
      }
        //now do zeta_vector
      for(j=0;j<3;++j){
        for(k=0;k<3;++k){
          zeta_vector[j] += ((z_mat[0])[k][j] - (z_mat[k+1])[k][j]);
        }
      }
    
        // accumulate this corner's contribution to the numerator
        // and the denominator...
      for(j=1;j<local_matrix_map_used;++j){
          //std::cout<<"c_k ("<<j<<") ="<<c_scalars[j]<<"\n";
          //std::cout<<"x_k ("<<j<<") ="<<verts[elems_verts[ (size_t) local_matrix_map[j]]]<<"\n";
        new_position+=verts[elems_verts[ (size_t) local_matrix_map[j]]]
          *c_scalars[j];
      }
      new_position+=zeta_vector;
    
        //std::cout<<"c_k (0) ="<<c_scalars[0]<<"\n";
        //std::cout<<"zeta ="<<zeta_vector<<"\n";
      scale_value+=c_scalars[0];
    
    }
    if(scale_value==0.0){
      MSQ_SETERR(err)("Invalid accummulation.\n",
                      MsqError::INVALID_STATE);
      return false;
    }
//.........这里部分代码省略.........

bool VertexConditionNumberQualityMetric::evaluate( PatchData& pd, 
                                                   size_t this_vert, 
                                                   double& fval, 
                                                   MsqError& err )
{
    //pd.generate_vertex_to_element_data();
  bool return_flag;
  fval=MSQ_MAX_CAP;
    //get the element array
  MsqMeshEntity* elems = pd.get_element_array(err);
    //get the vertex to element array and the offset array
  //const size_t* elem_offset = pd.get_vertex_to_elem_offset(err);  MSQ_ERRZERO(err);
  //const size_t* v_to_e_array = pd.get_vertex_to_elem_array(err);  MSQ_ERRZERO(err);
    //find the offset for this vertex
  //size_t this_offset = elem_offset[this_vert];
    //get the number of elements attached to this vertex (given by the
    //first entry in the vertex to element array)
  //size_t num_elems = v_to_e_array[this_offset];
    //PRINT_INFO("\nIN LOCAL SIZE CPP, num_elements = %i",num_elems);
    //if no elements, then return true
  size_t num_elems;
  const size_t *v_to_e_array = pd.get_vertex_element_adjacencies( this_vert, num_elems, err );
  MSQ_ERRZERO(err);  
    
  if(num_elems <= 0){
    return true;
  }
  
    //create an array to store the local metric values before averaging
    //Can we remove this dynamic allocatio?
  std::vector<double> met_vals(num_elems);
    //vector to hold the other verts which form a corner.
  vector<size_t> other_vertices;
  other_vertices.reserve(4);
  size_t i=0;
    //only 3 temp_vec will be sent to cond-num calculator, but the
    //additional vector3Ds may be needed during the calculations
  size_t elem_index;
  Vector3D temp_vec[6];
  const MsqVertex *vertices=pd.get_vertex_array(err);
  //loop over the elements attached to this vertex
  for(i=0;i<num_elems;++i){
      //get the vertices connected to this vertex for this element
    elem_index = v_to_e_array[i];
    elems[elem_index].get_connected_vertices(this_vert,
                                             other_vertices,
                                             err);  MSQ_ERRZERO(err);
      //switch over the element type of this element
    switch(elems[v_to_e_array[i]].get_element_type()){
    
      case TRIANGLE:
        temp_vec[0]=vertices[other_vertices[0]]-vertices[this_vert];
        temp_vec[2]=vertices[other_vertices[1]]-vertices[this_vert];
          //make relative to equilateral
        temp_vec[1]=((2*temp_vec[2])-temp_vec[0])*MSQ_SQRT_THREE_INV;
        return_flag=condition_number_2d(temp_vec,elem_index,pd,met_vals[i],err);  MSQ_ERRZERO(err);
        if(!return_flag)
          return return_flag;
        break;
      case QUADRILATERAL:
        temp_vec[0]=vertices[other_vertices[0]]-vertices[this_vert];
        temp_vec[1]=vertices[other_vertices[1]]-vertices[this_vert];
        return_flag=condition_number_2d(temp_vec,elem_index,pd,met_vals[i],err);  MSQ_ERRZERO(err);
        if(!return_flag)
          return return_flag;
        break;
      case TETRAHEDRON:
        temp_vec[0]=vertices[other_vertices[0]]-vertices[this_vert];
        temp_vec[3]=vertices[other_vertices[1]]-vertices[this_vert];
        temp_vec[4]=vertices[other_vertices[2]]-vertices[this_vert];
          //transform to equilateral tet
        temp_vec[1]=((2*temp_vec[3])-temp_vec[0])/MSQ_SQRT_THREE;
        temp_vec[2]=((3*temp_vec[4])-temp_vec[0]-temp_vec[3])/
          (MSQ_SQRT_THREE*MSQ_SQRT_TWO);
        return_flag=condition_number_3d(temp_vec,pd,met_vals[i],err);  MSQ_ERRZERO(err);
        if(!return_flag)
          return return_flag;
        break;
      case HEXAHEDRON:
        temp_vec[0]=vertices[other_vertices[0]]-vertices[this_vert];
        temp_vec[1]=vertices[other_vertices[1]]-vertices[this_vert];
        temp_vec[2]=vertices[other_vertices[2]]-vertices[this_vert];
        return_flag=condition_number_3d(temp_vec,pd,met_vals[i],err);  MSQ_ERRZERO(err);
        if(!return_flag)
          return return_flag;
        break;
      default:
        MSQ_SETERR(err)(MsqError::UNSUPPORTED_ELEMENT,
          "Element type (%d) not uspported in VertexConditionNumberQM.\n",
          (int)(elems[v_to_e_array[i]].get_element_type()));
        fval=MSQ_MAX_CAP;
        return false;
        
    }// end switch over element type
    other_vertices.clear();
  }//end loop over elements
  fval = average_metrics(arrptr(met_vals), num_elems, err);  MSQ_ERRZERO(err);
  return true;
}

bool LPtoPTemplate::concrete_evaluate(PatchData &pd, double &fval,
                                      MsqError &err){
  size_t index=0;
  MsqMeshEntity* elems=pd.get_element_array(err);
  bool obj_bool=true;
    //double check for pVal=0;
  if(pVal==0){
    MSQ_SETERR(err)("pVal equal zero not allowed.  L_0 is not a valid norm.",
                    MsqError::INVALID_STATE);
    return false;
  }
  
    //Michael:  this may not do what we want
    //Set currentQM to be the first quality metric* in the list
  QualityMetric* currentQM = get_quality_metric();
  if(currentQM==NULL)
    currentQM=get_quality_metric_list().front();
  if(currentQM==NULL) {
    MSQ_SETERR(err)("NULL QualityMetric pointer in LPtoPTemplate",
                    MsqError::INVALID_STATE);
    return false;
  }
  size_t num_elements=pd.num_elements();
  size_t num_vertices=pd.num_vertices();
  size_t total_num=0;
  if(currentQM->get_metric_type()==QualityMetric::ELEMENT_BASED)   
    total_num=num_elements;
  else if (currentQM->get_metric_type()==QualityMetric::VERTEX_BASED)
    total_num=num_vertices;
  else {
    MSQ_SETERR(err)("Make sure MetricType is initialised in concrete "
                    "QualityMetric constructor.", MsqError::INVALID_STATE);
    return false;
  }
  
  msq_std::vector<double> metric_values(total_num);
  if(currentQM->get_metric_type()==QualityMetric::ELEMENT_BASED)
  {
    for (index=0; index<num_elements;index++)
    {
        //if invalid return false after clean-up
      obj_bool=currentQM->evaluate_element(pd, (&elems[index]),
                                           metric_values[index], err);
      if(MSQ_CHKERR(err) || !obj_bool){
        fval=0.0;
        return false;
      }
      
      metric_values[index]=fabs(metric_values[index]);
      MSQ_DBGOUT(3) << "      o  Quality metric value for element "
                    << index << "\t: " << metric_values[index] << "\n";
    }
  }
  else if(currentQM->get_metric_type()==QualityMetric::VERTEX_BASED)
  {
    MsqVertex* vertices=pd.get_vertex_array(err); MSQ_ERRZERO(err);
    for (index=0; index<num_vertices;index++)
    {
        //evaluate metric for this vertex
      obj_bool=currentQM->evaluate_vertex(pd, (&vertices[index]),
                                          metric_values[index], err);
        //if invalid return false after clean-up
      if(MSQ_CHKERR(err) || !obj_bool){
        fval=0.0;
        return false;
      }
      
      metric_values[index]=fabs(metric_values[index]);
    }
  }
  fval=compute_function(&metric_values[0], total_num, err);
  return !MSQ_CHKERR(err);
}

/*! \fn LPtoPTemplate::compute_analytical_hessian(PatchData &pd, MsqHessian &hessian, MsqError &err)

    For each element, each entry to be accumulated in the Hessian for
    this objective function (\f$ \sum_{e \in E} Q(e)^p \f$ where \f$ E \f$
    is the set of all elements in the patch) has the form:
    \f$ pQ(e)^{p-1} \nabla^2 Q(e) + p(p-1)Q(e)^{p-2} \nabla Q(e) [\nabla Q(e)]^T \f$.

    For \f$ p=2 \f$, this simplifies to
    \f$ 2Q(e) \nabla^2 Q(e) + 2 \nabla Q(e) [\nabla Q(e)]^T \f$.

    For \f$ p=1 \f$, this simplifies to \f$ \nabla^2 Q(e) \f$.

    The \f$ p=1 \f$ simplified version is implemented directly
    to speed up computation. 

    This function does not support vertex-based metrics.
    
    \param pd The PatchData object for which the objective function
           hessian is computed.
    \param hessian: this object must have been previously initialized.
*/
bool LPtoPTemplate::compute_analytical_hessian(PatchData &pd,
                                               MsqHessian &hessian,
                                               Vector3D *const &grad,
                                               double &OF_val,
                                               MsqError &err)
{
  double scaling_value=1.0;
  
  MSQ_FUNCTION_TIMER( "LPtoPTemplate::compute_analytical_hessian" );

  MsqMeshEntity* elements = pd.get_element_array(err); MSQ_ERRZERO(err);
  MsqVertex* vertices = pd.get_vertex_array(err); MSQ_ERRZERO(err);
  size_t num_elems = pd.num_elements();
    //if scaling divide by the number of elements.
  if(dividingByN){
    if(num_elems<=0) {
      MSQ_SETERR(err)("LPtoP is attempting to divide by zero in analytical Hessian.",
                      MsqError::INVALID_MESH);
      return false;
    }
    scaling_value/=num_elems;
  }
  
  size_t num_vertices = pd.num_vertices();
  Matrix3D elem_hessian[MSQ_MAX_NUM_VERT_PER_ENT*(MSQ_MAX_NUM_VERT_PER_ENT+1)/2];
  Matrix3D elem_outer_product;
  Vector3D grad_vec[MSQ_MAX_NUM_VERT_PER_ENT];
  double QM_val;
  double fac1, fac2;
  Matrix3D grad_outprod;
  bool qm_bool;
  QualityMetric* currentQM = get_quality_metric();
  
  MsqVertex* ele_free_vtces[MSQ_MAX_NUM_VERT_PER_ENT];
  short i;
  for (i=0; i<MSQ_MAX_NUM_VERT_PER_ENT; ++i) ele_free_vtces[i]=NULL;
  
  const size_t* vtx_indices;
    
  size_t e, v;
  size_t nfve; // number of free vertices in element
  short j,n;

  hessian.zero_out();
  for (v=0; v<num_vertices; ++v) grad[v] = 0.;
  OF_val = 0.;
  
  // Loops over all elements in the patch.
  for (e=0; e<num_elems; ++e) {
    short nve = elements[e].vertex_count();
    
    // Gets a list of free vertices in the element.
    vtx_indices = elements[e].get_vertex_index_array();
    nfve=0;
    for (i=0; i<nve; ++i) {
      if ( vertices[vtx_indices[i]].is_free_vertex() ) {
        ele_free_vtces[nfve] = vertices + vtx_indices[i];
        ++nfve;
      }
    }
    
    // Computes \nabla^2 Q(e). Only the free vertices will have non-zero entries. 
    qm_bool = currentQM->compute_element_hessian(pd,
                                    elements+e, ele_free_vtces,
                                    grad_vec, elem_hessian,
                                    nfve, QM_val, err);
    if (MSQ_CHKERR(err) || !qm_bool) return false;


    // **** Computes Hessian ****
    double QM_pow=1.;
    if (pVal == 1) {
      n=0;
      for (i=0; i<nve; ++i) {
        for (j=i; j<nve; ++j) {
            //negate if necessary
          elem_hessian[n] *= (scaling_value * get_negate_flag());
          ++n;
        }
//.........这里部分代码省略.........

/* 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());
//.........这里部分代码省略.........