C++ PatchData::get_vertex_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 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 SmartLaplacianSmoother::optimize_vertex_positions( PatchData &pd, 
                                                        MsqError &err )
{
  assert(pd.num_free_vertices() == 1);
  const size_t center_vtx_index = 0;
  const size_t init_inverted = num_inverted( pd, err ); MSQ_ERRRTN(err);
  
  adjVtxList.clear();
  pd.get_adjacent_vertex_indices( center_vtx_index, adjVtxList, err );
  MSQ_ERRRTN(err);
  
  if (adjVtxList.empty())
    return;
  
  const MsqVertex* verts = pd.get_vertex_array(err);
  const size_t n = adjVtxList.size();
  
  const Vector3D orig_pos = verts[center_vtx_index];
  Vector3D new_pos = verts[ adjVtxList[0] ];
  for (size_t i = 1; i < n; ++i)
    new_pos += verts[ adjVtxList[i] ];
  new_pos *= 1.0/n;
  pd.set_vertex_coordinates( new_pos, center_vtx_index, err );
  pd.snap_vertex_to_domain( center_vtx_index, err );  MSQ_ERRRTN(err);
  const size_t new_inverted = num_inverted( pd, err ); MSQ_ERRRTN(err);
  if (new_inverted > init_inverted)
    pd.set_vertex_coordinates( orig_pos, center_vtx_index, err );
}

bool EdgeLengthQualityMetric::evaluate_vertex(PatchData &pd, MsqVertex* vert,
                                             double &fval, MsqError &err)
{
  fval=0.0;
  size_t this_vert = pd.get_vertex_index(vert);
  size_t other_vert;
  vector<size_t> adj_verts;
  Vector3D edg;
  pd.get_adjacent_vertex_indices(this_vert,adj_verts,err);  MSQ_ERRZERO(err);
  int num_sample_points=adj_verts.size();
  double *metric_values=new double[num_sample_points];
  MsqVertex* verts = pd.get_vertex_array(err);  MSQ_ERRZERO(err);
  int point_counter=0;
  while(!adj_verts.empty()){
    other_vert=adj_verts.back();
    adj_verts.pop_back();
    edg[0]=verts[this_vert][0]-verts[other_vert][0];
    edg[1]=verts[this_vert][1]-verts[other_vert][1];
    edg[2]=verts[this_vert][2]-verts[other_vert][2];    
    metric_values[point_counter]=edg.length();
    ++point_counter;
  }
  fval=average_metrics(metric_values,num_sample_points,err);  MSQ_ERRZERO(err);
  delete[] metric_values;
  
  return true;
  
}

void NonGradient::initialize_mesh_iteration(PatchData &pd, MsqError &err)
{
  int elementDimension = getPatchDimension( pd, err );  // to do: react to error
  int dimension = elementDimension * pd.num_free_vertices();
  //printPatch( pd, err );
  setDimension(dimension);
  int maxNumEval = 100*dimension;  // 1. Make this a user parameter
  setMaxNumEval(maxNumEval);
  double threshold = 1.e-10; // avoid division by zero
  setThreshold(threshold);
  double minEdgeLen = 0.0;
  double maxEdgeLen = 0.0;
//  double ftol = 0.;
  if( dimension > 0 )
  {
    pd.get_minmax_edge_length( minEdgeLen, maxEdgeLen );
    //ftol = minEdgeLen * 1.e-4; // Turn off Amoeba convergence criterion
    if( mNonGradDebug >= 1 ) 
    {      
         std::cout << "minimum edge length " << minEdgeLen << " maximum edge length " << maxEdgeLen << std::endl;
    }      
    MSQ_PRINT(3)("minimum edge length %e    maximum edge length %e\n", minEdgeLen,  maxEdgeLen);
  }
//  setTolerance(ftol);
  int numRow = dimension;
  int numCol = numRow+1;  
  if( numRow*numCol <= simplex.max_size() )
  { 
    simplex.assign(numRow*numCol, 0.);  // guard against previous simplex value
    double scale = minEdgeLen * mScaleDiameter;; 
    const MsqVertex* coord = pd.get_vertex_array(err);
    if( pd.num_free_vertices() > 1 )
    {
      MSQ_SETERR(err)("Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED);
    }
    size_t index = 0;
    for( int col = 0; col < numCol; col++ )
    {
      for (int row=0;row<numRow;row++)
      {
        simplex[ row + col*numRow ] = coord[index][row];
        if( row == col-1 )
        {
          simplex[ row + col*numRow ] += scale/ static_cast<double>(numCol);
        }
      }
    }
  }
  else
  {
    MSQ_SETERR(err)("Use patch with fewer free vertices", MsqError::OUT_OF_MEMORY);
    if( mNonGradDebug >= 1 ) 
    {      
      std::cout << "ERROR: Too many free vertices in patch" << std::endl;
    }      
    //MSQ_PRINT(1)("ERROR: Too many free vertices in patch\n");
  }
}

   //!Calculate the area of the triangle formed by the three vertices.
static inline double compute_corner_area( PatchData &pd,
                                          size_t vert_1,
                                          size_t vert_2,
                                          size_t vert_3,
                                          MsqError &err)
{
  const MsqVertex* verts = pd.get_vertex_array(err);
  Vector3D vec_1=verts[vert_2]-verts[vert_1];
  Vector3D vec_2=verts[vert_3]-verts[vert_1];
  Vector3D cross_vec=vec_1*vec_2;
  return (cross_vec.length()/2.0);
}

   //!Calculate the volume of the tetrahedron formed by the four vertices.
static inline double compute_corner_volume( PatchData &pd,
                                            size_t vert_1,
                                            size_t vert_2,
                                            size_t vert_3,
                                            size_t vert_4,
                                            MsqError &err)
{
  const MsqVertex* verts = pd.get_vertex_array(err);
  Vector3D vec_1=verts[vert_2]-verts[vert_1];
  Vector3D vec_2=verts[vert_3]-verts[vert_1];
  Vector3D vec_3=verts[vert_4]-verts[vert_1];
  return fabs((vec_3%(vec_1*vec_2))/6.0);

}  

/*! \todo Michael:  optimize_vertex_position is probably not implemented
  in an optimal way.  We used to use all of the vertices in
  the patch as 'adjacent' vertices.  Now we call get_adjacent_vertex_indices.
  We could use a VERTICES_ON_VERTEX type of patch or a global patch?
*/
void SmartLaplacianSmoother::optimize_vertex_positions(PatchData &pd, 
                                                MsqError &err)
{
    //default the laplacian smoother to 3 even for 2-d elements.
    //int dim = get_mesh_set()->space_dim();
  size_t dim = 3;
  MsqVertex* verts=pd.get_vertex_array(err);  MSQ_ERRRTN(err);
    //the objective function evaluator
  OFEvaluator& obj_func = get_objective_function_evaluator();
    //variables for the function values.
  double orig_val=0;
  double mod_val=0;
    //compute the original function value and check validity
  bool valid_flag = obj_func.evaluate(pd,orig_val,err);  MSQ_ERRRTN(err);
  // does the Laplacian smoothing
  MsqFreeVertexIndexIterator free_iter(pd, err);  MSQ_ERRRTN(err);
  free_iter.reset();
  free_iter.next();
    //m is the free vertex.
  size_t m=free_iter.value();
  vector<size_t> vert_indices;
  vert_indices.reserve(25);
    //get vertices adjacent to vertex m
  pd.get_adjacent_vertex_indices(m,vert_indices,err);  MSQ_ERRRTN(err);
    //move vertex m
    //save the original position of the free vertex
  Vector3D orig_position(verts[m]);
    //smooth the patch
  centroid_smooth_mesh(pd, vert_indices.size(), vert_indices,
                       m, dim, err); MSQ_ERRRTN(err);
    //snap vertex m to domain
  pd.snap_vertex_to_domain(m,err);  MSQ_ERRRTN(err);
    //if the original function val was invalid, then we allow the move
    //But, if it wasn valid, we need to decide.
  if(valid_flag){
      //compute the new value
    valid_flag = obj_func.evaluate(pd,mod_val,err);  MSQ_ERRRTN(err);
      //if the new value is worse the original OR if the new value is not
      //valid (we already know the original value was valid by above) then
      //we don't allow the move.
    if(!valid_flag || mod_val>orig_val){
        //move the vert back to where it was.
      verts[m]=orig_position;
        //PRINT_INFO("\norig = %f, new = %f, new valid = %d",orig_val,mod_val,valid_flag);
    }
    
  }
  
}

 inline void move_vertex( PatchData& pd, 
                          const Vector3D& position,
                          const Vector3D& delta,
                          MsqError& err )
 {
   const MsqVertex* array = pd.get_vertex_array( err ); 
   if (err) return;
   
   int idx = 0, cnt = 0;
   for (size_t i = 0; i < pd.num_nodes(); ++i)
     if ((array[i] - position).length_squared() < DBL_EPSILON) 
       { idx = i; ++cnt; }
   
   CPPUNIT_ASSERT_EQUAL( cnt, 1 );
   
   pd.move_vertex( delta, idx, err );
 }

void Randomize::optimize_vertex_positions(PatchData &pd, 
                                                MsqError &err)
{
    //cout << "- Executing Randomize::optimize_vertex_position()\n";

  int num_local_vertices = pd.num_vertices();
  // gets the array of coordinates for the patch and print it 
  MsqVertex *patch_coords = pd.get_vertex_array(err); MSQ_ERRRTN(err);
  // does the randomize smooth
  MsqFreeVertexIndexIterator free_iter(&pd, err); MSQ_ERRRTN(err);
  free_iter.reset();
  free_iter.next();
    //find the free vertex.
  int m=free_iter.value();
  randomize_vertex(pd, num_local_vertices,
                   patch_coords[m], err);  MSQ_ERRRTN(err);
  pd.snap_vertex_to_domain(m,err); MSQ_ERRRTN(err);
}

  void test_hex_vertices()
  {
    MsqPrintError err(cout);
    // prints out the vertices.
    const MsqVertex* ideal_vertices = oneHexPatch.get_vertex_array(err); CPPUNIT_ASSERT(!err);
    size_t num_vtx = oneHexPatch.num_nodes();
    CPPUNIT_ASSERT_EQUAL(size_t(8), num_vtx);
    
    MsqVertex vtx;

    vtx.set(1,1,1);
    CPPUNIT_ASSERT_EQUAL(vtx, ideal_vertices[0]);
    
    vtx.set(2,2,2);
    CPPUNIT_ASSERT_EQUAL(vtx, ideal_vertices[6]);
    
    vtx.set(1,2,2);
    CPPUNIT_ASSERT_EQUAL(vtx, ideal_vertices[7]);
  }

/*!For the given vertex, vert, with connected edges of lengths l_j for
  j=1...k, the metric value is the average (where the default average
  type is SUM) of
        u_j = ( | l_j - lowVal | - (l_j - lowVal) )^2 +
              ( | highVal - l_j | - (highVal - l_j) )^2.
*/
bool EdgeLengthRangeQualityMetric::evaluate_common(PatchData &pd, 
                                             size_t this_vert,
                                             double &fval, 
                                             std::vector<size_t>& adj_verts,
                                             MsqError &err)
{
  fval=0.0;
  Vector3D edg;
  pd.get_adjacent_vertex_indices(this_vert,adj_verts,err);  MSQ_ERRZERO(err);
  int num_sample_points=adj_verts.size();
  double *metric_values=new double[num_sample_points];
  const MsqVertex* verts = pd.get_vertex_array(err);  MSQ_ERRZERO(err);
    //store the length of the edge, and the first and second component of
    //metric values, respectively.
  double temp_length=0.0;
  double temp_first=0.0;
  double temp_second=0.0;
    //PRINT_INFO("INSIDE ELR, vertex = %f,%f,%f\n",verts[this_vert][0],verts[this_vert][1],verts[this_vert][2]);
    //loop while there are still more adjacent vertices.
  for (unsigned i = 0; i < adj_verts.size(); ++i) 
  {
    edg = verts[this_vert] - verts[adj_verts[i]];
      //compute the edge length
    temp_length=edg.length();
      //get the first component
    temp_first = temp_length - lowVal;
    temp_first = fabs(temp_first) - (temp_first);
    temp_first*=temp_first;
      //get the second component
    temp_second = highVal - temp_length;
    temp_second = fabs(temp_second) - (temp_second);
    temp_second*=temp_second;
      //combine the two components
    metric_values[i]=temp_first+temp_second;
  }
    //average the metric values of the edges
  fval=average_metrics(metric_values,num_sample_points,err);
    //clean up
  delete[] metric_values;
    //always return true because mesh is always valid wrt this metric.
  return !MSQ_CHKERR(err);
  
}

void TerminationCriterion::accumulate_patch( PatchData& pd, MsqError& err )
{
  if (terminationCriterionFlag & MOVEMENT_FLAGS)
  {
    double patch_max_dist = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
    if (patch_max_dist > maxSquaredMovement)
      maxSquaredMovement = patch_max_dist;
    pd.recreate_vertices_memento( previousVerticesMemento, err );  MSQ_ERRRTN(err);
  }
    
    //if terminating on bounded vertex movement (a bounding box for the mesh)
  if(terminationCriterionFlag & BOUNDED_VERTEX_MOVEMENT)
  {
    const MsqVertex* vert = pd.get_vertex_array(err);
    int num_vert = pd.num_free_vertices();
    int i=0;
      //for each vertex
    for(i=0;i<num_vert;++i)
    {
        //if any of the coordinates are greater than eps
      if( (vert[i][0]>boundedVertexMovementEps) ||
          (vert[i][1]>boundedVertexMovementEps) ||
          (vert[i][2]>boundedVertexMovementEps) )
      {
        ++vertexMovementExceedsBound;
      }
    }
  }


  if ((terminationCriterionFlag|cullingMethodFlag) & UNTANGLED_MESH) {
    size_t new_count = count_inverted( pd, err );
      // be careful here because size_t is unsigned
    globalInvertedCount += new_count;
    globalInvertedCount -= patchInvertedCount;
    patchInvertedCount = new_count;
    //if (innerOuterType==TYPE_OUTER) 
    //  MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Num Patch Inverted: " << " " << patchInvertedCount << " globalInvertedCount= " << globalInvertedCount << std::endl;
      
    MSQ_ERRRTN(err);
  }
}

/*!Function calculates a scale factor for the patch, then moves
  the incident vertex randomly in each of the three coordinate
  directions (relative to the scale factor multiplied by mPercent).
*/
static inline void randomize_vertex(PatchData &pd,
                                    size_t free_ind,
                                    double percent,
                                    MsqError &err)
{
    size_t i;
    short j;
    const MsqVertex* verts = pd.get_vertex_array(err);
    MSQ_ERRRTN(err);
    const size_t num_vtx = pd.num_nodes();
    //a scale w.r.t. the patch size
    double scale_factor=0.0;
    //a "random" number between -1 and 1
    double rand_double=0.0;
    //a "random" int
    int rand_int=0;
    if (num_vtx<=1) {
        MSQ_PRINT(1)("WARNING: Number of incident vertex is zero.  Returning.\n");
        return;
    }

    for (i=0; i<num_vtx; ++i) {
        if(i != free_ind)
            scale_factor+=(verts[i]-verts[free_ind]).length();
    }
    scale_factor/=( (double) num_vtx - 1.0 );
    Vector3D delta;
    for (j=0; j<3; ++j) {
        rand_int = rand();
        //number between 0 and 1000
        rand_int = rand_int%1000;
        //number between -1 and 1
        rand_double = (((double) rand_int)/500.0)-1.0;
        delta[j] = scale_factor * rand_double * percent;
    }
    pd.move_vertex( delta, free_ind, err );

    return;
}

double TetDihedralWeight::get_weight( PatchData& pd, 
                                      size_t element,
                                      Sample ,
                                      MsqError& err )
{
  const double eps = 1e-10;

  MsqMeshEntity &elem = pd.element_by_index( element );
  if (elem.get_element_type() != TETRAHEDRON) {
    MSQ_SETERR(err)(MsqError::UNSUPPORTED_ELEMENT);
    return 0.0;
  }

  const size_t *indices = elem.get_vertex_index_array();
  Vector3D v01, v02, v31, v32;
  if (refMesh) {
    const Mesh::VertexHandle* vtx_hdl = pd.get_vertex_handles_array();
    Mesh::VertexHandle handles[] = { vtx_hdl[indices[0]],
                                     vtx_hdl[indices[1]],
                                     vtx_hdl[indices[2]],
                                     vtx_hdl[indices[3]] };
    Vector3D coords[4];
    refMesh->get_reference_vertex_coordinates( handles, 4, coords, err );
    MSQ_ERRZERO(err);

    v01 = coords[1] - coords[0];
    v02 = coords[2] - coords[0];
    v31 = coords[1] - coords[3];
    v32 = coords[2] - coords[3];
  }
  else {  
    const MsqVertex* coords = pd.get_vertex_array();
    v01 = coords[indices[1]] - coords[indices[0]];
    v02 = coords[indices[2]] - coords[indices[0]];
    v31 = coords[indices[1]] - coords[indices[3]];
    v32 = coords[indices[2]] - coords[indices[3]];
  }
  
  Vector3D n012 = v02 * v01;
  Vector3D n013 = v31 * v01;
  Vector3D n023 = v02 * v32;
  Vector3D n123 = v31 * v32;
  
    // normalize face vectors.
  double l012 = n012.length();
  double l013 = n013.length();
  double l023 = n023.length();
  double l123 = n123.length();
  n012 *= (l012 < eps) ? 0.0 : 1.0/l012;
  n013 *= (l013 < eps) ? 0.0 : 1.0/l013;
  n023 *= (l023 < eps) ? 0.0 : 1.0/l023;
  n123 *= (l123 < eps) ? 0.0 : 1.0/l123;
  
    // calculate dihedral handles for each edge
  double ds[] = { da(n012 % n013),
                  da(n012 % n123),
                  da(n012 % n023),
                  da(n013 % n023),
                  da(n013 % n123),
                  da(n023 % n123) };
  
    // calculate weight from max dihedral handle
  double d = *std::max_element( ds, ds+6 );
  return 1/(1 + exp(-mA*(d - mCutoff)));
}