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

bool
TMPQualityMetric::evaluate_surface_common( PatchData& pd,
        Sample s,
        size_t e,
        const NodeSet& bits,
        size_t* indices,
        size_t& num_indices,
        MsqVector<2>* derivs,
        MsqMatrix<2,2>& W,
        MsqMatrix<2,2>& A,
        MsqMatrix<3,2>& S_a_transpose_Theta,
        MsqError& err )
{
    EntityTopology type = pd.element_by_index( e ).get_element_type();

    const MappingFunction2D* mf = pd.get_mapping_function_2D( type );
    if (!mf) {
        MSQ_SETERR(err)( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
        return false;
    }

    MsqMatrix<3,2> J;
    mf->jacobian( pd, e, bits, s, indices, derivs, num_indices, J, err );

    // If we have a 3x2 target matrix
    if (targetCalc->have_surface_orient()) {
        MsqVector<3> u, u_perp;
        MsqMatrix<3,2> W_hat;
        targetCalc->get_surface_target( pd, e, s, W_hat, err );
        MSQ_ERRZERO(err);
        // Use the cross product of the columns of W as the normal of the
        // plane to work in (i.e. u.).  W should have been constructed such
        // that said cross product is in the direction of (n_s)_init.  And if
        // for some reason it as not, then using something other than said
        // cross product is likely to produce very wrong results.
        project_to_matrix_plane( W_hat, W, u, u_perp );
        // Do the transforms on A to align it with W and project into the plane.
        if (!project_to_perp_plane( J, u, u_perp, A, S_a_transpose_Theta ))
            return false;
    }
    // Otherwise if we have a 2x2 target matrix (i.e. the target does
    // not contain orientation information), project into the plane
    // tangent to J.
    else {
        MsqVector<3> u, u_perp;
        targetCalc->get_2D_target( pd, e, s, W, err );
        MSQ_ERRZERO(err);
        project_to_matrix_plane( J, A, u, u_perp );
        S_a_transpose_Theta.set_column(0, u_perp);
        S_a_transpose_Theta.set_column(1, u*u_perp);
        // If the domain is set, adjust the sign of things correctly
        // for the case where the element is inverted with respect
        // to the domain.
        if (pd.domain_set()) {
            Vector3D n;
            pd.get_domain_normal_at_sample( e, s, n, err );
            MSQ_ERRZERO(err);
            // if sigma == -1
            if (Vector3D(u.data()) % n < 0.0) {
                // flip u
                u = -u;
                // S_a_transpose_Theta == Theta, because S_a == I here.
                // u_perp is unaffected by flipping u, so only the second
                // column of S_a_transpose_Theta and the second row of A
                // are flipped because u x u_perp will be flipped.
                S_a_transpose_Theta.set_column(1, -S_a_transpose_Theta.column(1) );
                A.set_row( 1, -A.row(1) );
            }
        }
    }

    return true;
}

void FeasibleNewton::optimize_vertex_positions(PatchData &pd, 
                                               MsqError &err)
{
  MSQ_FUNCTION_TIMER( "FeasibleNewton::optimize_vertex_positions" );
  MSQ_DBGOUT(2) << "\no  Performing Feasible Newton optimization.\n";
  
  //
  // the only valid 2D meshes that FeasibleNewton works for are truly planar which 
  // lie in the X-Y coordinate plane.
  //

  XYPlanarDomain *xyPlanarDomainPtr = dynamic_cast<XYPlanarDomain*>(pd.get_domain());
    // only optimize if input mesh is a volume or an XYPlanarDomain
  if (!pd.domain_set() || xyPlanarDomainPtr != NULL)  
  {
    const double sigma   = 1e-4;
    const double beta0   = 0.25;
    const double beta1   = 0.80;
    const double tol1    = 1e-8;
    const double tol2    = 1e-12;
    const double epsilon = 1e-10;
    double original_value, new_value;
    double beta;
  
    int nv = pd.num_free_vertices();
    std::vector<Vector3D> grad(nv), d(nv);
    bool fn_bool=true;// bool used for determining validity of patch
  
    OFEvaluator& objFunc = get_objective_function_evaluator();

    int i;

    // TODD -- Don't blame the code for bad things happening when using a
    //         bad termination test or requesting more accuracy than is
    //	     possible.
    //
    //         Also, 

    // 1.  Allocate a hessian and calculate the sparsity pattern.
    mHessian.initialize(pd, err); MSQ_ERRRTN(err);
  
    // does the Feasible Newton iteration until stopping is required.
    // Terminate when inner termination criterion signals.

    /* Computes the value of the stopping criterion*/
    TerminationCriterion* term_crit=get_inner_termination_criterion();
    while ( !term_crit->terminate() ) {
      fn_bool = objFunc.update( pd, original_value, grad, mHessian,  err );
      MSQ_ERRRTN(err);
      if (!fn_bool) {
        MSQ_SETERR(err)("invalid patch for hessian calculation", MsqError::INTERNAL_ERROR);
        return; 
      }

      if (MSQ_DBG(3)) { // avoid expensive norm calculations if debug flag is off
        MSQ_DBGOUT(3) << "  o  objective function: " << original_value << std::endl;
        MSQ_DBGOUT(3) << "  o  gradient norm: " << length(grad) << std::endl;
        MSQ_DBGOUT(3) << "  o  Hessian norm: " << mHessian.norm() << std::endl;
      }
    
      // Prints out free vertices coordinates. 
      //
      // Comment out the following because it is way to verbose for larger
      // problems.  Consider doing:
      //  inner_term_crit->write_mesh_steps( "filename.vtk" );
      // instead.
      // - j.kraftcheck 2010-11-17
//      if (MSQ_DBG(3)) {
//        MSQ_DBGOUT(3) << "\n  o Free vertices ("<< pd.num_free_vertices()
//                  <<")original coordinates:\n ";
//        MSQ_ERRRTN(err);
//        const MsqVertex* toto1 = pd.get_vertex_array(err); MSQ_ERRRTN(err);
//        MsqFreeVertexIndexIterator ind1(pd, err); MSQ_ERRRTN(err);
//        ind1.reset();
//        while (ind1.next()) {
//          MSQ_DBGOUT(3) << "\t\t\t" << toto1[ind1.value()];
//        }
//      }
      
      // 4. Calculate a direction using preconditionned conjugate gradients
      //    to find a zero of the Newton system of equations (H*d = -g)
      //    (a) stop if conjugate iteration limit reached
      //    (b) stop if relative residual is small
      //    (c) stop if direction of negative curvature is obtained

      mHessian.cg_solver(arrptr(d), arrptr(grad), err); MSQ_ERRRTN(err);

      // 5. Check for descent direction (inner produce of gradient and
      //    direction is negative.
      double alpha = inner( grad, d );
      // TODD -- Add back in if you encounter problems -- do a gradient
      //         step if the direction from the conjugate gradient solver
      //         is not a descent direction for the objective function.  We
      //         SHOULD always get a descent direction from the conjugate
      //         method though, unless the preconditioner is not positive
      //         definite.
      // If direction is positive, does a gradient (steepest descent) step.

      if (alpha > -epsilon) {
    
//.........这里部分代码省略.........