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

double 
NonGradient::evaluate( double *point,  PatchData &pd, MsqError &err )
{
  if( pd.num_free_vertices() > 1 )
  {
    MSQ_SETERR(err)("Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED);
  }
  const size_t vertexIndex = 0; 
  Vector3D originalVec = pd.vertex_by_index(vertexIndex);
  Vector3D pointVec;
  for( int index = 0; index<3; index++)
  {
    pointVec[index] = point[index];
  }
  pd.set_vertex_coordinates( pointVec, vertexIndex, err ); 
  pd.snap_vertex_to_domain( vertexIndex, err );  

  OFEvaluator& obj_func = get_objective_function_evaluator();
  TerminationCriterion* term_crit=get_inner_termination_criterion();

  double value;
  bool feasible = obj_func.evaluate( pd, value, err ); // MSQ_ERRZERO(err);
  term_crit->accumulate_inner( pd, value, 0, err );  //MSQ_CHKERR(err);
  if (err.error_code() == err.BARRIER_VIOLATED)
    err.clear();   // barrier violation not an error here
  MSQ_ERRZERO(err);
  pd.set_vertex_coordinates( originalVec, vertexIndex, err ); 
  if( !feasible && mUseExactPenaltyFunction )  
  { // "value" undefined btw
    double ensureFiniteRtol= .25;
    value = DBL_MAX * ensureFiniteRtol;
    //std::cout << " Infeasible patch: " << value << std::endl; printPatch( pd, err );
  }
  return value;
}

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 );
}

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

/*!
 */
  void I_DFT_NoBarrierSmoother::optimize_vertex_positions(PatchData &pd, 
                                                          MsqError &err)
  {
    MSQ_FUNCTION_TIMER( "I_DFT_NoBarrierSmoother::optimize_vertex_positions" );
    
      // does the I_DFT_NoBarrier 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();
      //move vertex m
    i_dft_no_barrier_smooth_mesh(pd, m, err); MSQ_ERRRTN(err);
      //snap vertex m to domain
    pd.snap_vertex_to_domain(m,err);
  
  }

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

    // gets the array of coordinates for the patch and print it
    // 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, m, mPercent, err);
    MSQ_ERRRTN(err);
    pd.snap_vertex_to_domain(m,err);
    MSQ_ERRRTN(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);
}

//.........这里部分代码省略.........
  
/*
      if( height[0] == 0.)
      {
         MSQ_SETERR(err)("(B) Zero objective function value", MsqError::INTERNAL_ERROR);
         exit(-1);
      }
*/
      if (ytry <= height[ilo])   
      {
        ytry=amotry(simplex,height,&rowSum[0],ihi,-2.0,pd,err);
        if( mNonGradDebug >= 3 ) 
        {      
         std::cout << "Reflect and Expand from highPt " << ytry << std::endl;
        }      
        //MSQ_PRINT(3)("Reflect and Expand from highPt : %e\n",ytry);
      }
      else 
      {
        if (ytry >= height[inhi]) 
        {
          ysave=height[ihi]; // Contract along highPt
          ytry=amotry(simplex,height,&rowSum[0],ihi,0.5,pd,err);
          if (ytry >= ysave)
          { // contract all directions toward lowPt
            for (int col=0;col<numCol;col++)
            {
              if (col != ilo)
              {
                for (int row=0;row<numRow;row++)
                {
                  rowSum[row]=0.5*(simplex[row+col*numRow]+simplex[row+ilo*numRow]);
                  simplex[row+col*numRow]=rowSum[row];
                }
                height[col] = evaluate(&rowSum[0], pd, err); 
                if( mNonGradDebug >= 3 ) 
                {      
                  std::cout << "Contract all directions toward lowPt value( " << col << " ) = " << height[col] << " ilo = " << ilo << std::endl;
                }      
                //MSQ_PRINT(3)("Contract all directions toward lowPt value( %d ) = %e    ilo = %d\n", col, height[col], ilo);
              }
            }
          }
        }   
      } // ytri > h(ilo) 
    } // after evaluation
    ilo=1; // conditional operator or inline if 
    ihi = height[0] > height[1] ? (inhi=1,0) : (inhi=0,1);
    for (int col=0;col<numCol;col++)
    {
      if (height[col] <= height[ilo])
      {
        ilo=col;  // ilo := argmin height
      }
      if (height[col] > height[ihi])
      {
        inhi=ihi;
        ihi=col;
      } 
      else  // height[ihi] >= height[col]
        if (col != ihi && height[col] > height[inhi] ) inhi=col;
    }
//    rtol=2.0*fabs( height[ihi]-height[ilo] )/
//         ( fabs(height[ihi])+fabs(height[ilo])+threshold );
    afterEvaluation = true;
  } //  while not converged 

  // Always set to current best mesh.
  { 
    if( ilo != 0 )
    {
      double yTemp = height[0];
      height[0] = height[ilo]; // height dimension numCol
      height[ilo] = yTemp;
      for (int row=1;row<numRow;row++)
      { 
          yTemp = simplex[row];
          simplex[row] = simplex[row+ilo*numRow];
          simplex[row+ilo*numRow] = yTemp;
      }
    }
  }
  if( pd.num_free_vertices() > 1 )
  {
    MSQ_SETERR(err)("Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED);
  }

  Vector3D newPoint( &simplex[0] ); 
  size_t vertexIndex = 0; // fix c.f. freeVertexIndex
  pd.set_vertex_coordinates( newPoint, vertexIndex, err ); 
  pd.snap_vertex_to_domain( vertexIndex, err );
  if( term_crit->terminate() )
  {
    if( mNonGradDebug >= 1 ) 
    {      
         std::cout << "Optimization Termination OptStatus: Max Iter Exceeded" << std::endl;
    }      
    //MSQ_PRINT(1)("Optimization Termination OptStatus: Max Iter Exceeded\n"); 
  }
}