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

bool PatchPowerMeanP::initialize_block_coordinate_descent( Mesh* mesh, 
                                                      MeshDomain* domain, 
                                                      const Settings* settings,
                                                      PatchSet* patch_set,
                                                      MsqError& err )
{
  clear();
  PatchIterator patches( patch_set );
  
  PatchData pd;
  pd.set_mesh( mesh );
  pd.set_domain( domain );
  if (settings)
    pd.attach_settings(settings);
  
  bool result = true;
  while (patches.get_next_patch( pd, err ) && !MSQ_CHKERR(err))
  {
    double value;
    bool b = evaluate( ObjectiveFunction::ACCUMULATE, pd, value, false, err ); 
    MSQ_ERRZERO(err);
    result = result && b;
  }
  return result;
}

   /*! \fn create_six_quads_patch(PatchData &four_quads, MsqError &err)
     our 2D set up: 6 quads, 1 center vertex outcentered by (0,-0.5), the other centered
      7____6____5___11
      |    |    |    |
      | 2  |  3 | 5  |
      8-_  |  _-4---10       vertex 1 is at (0,0)
      |  -_0_-  |    |       vertex 11 is at (3,2)
      | 0  |  1 | 4  |
      1----2----3----9

      use destroy_patch_with_domain() in sync.
   */
   inline void create_six_quads_patch_with_domain(PatchData &pd, MsqError &err) 
   {
     // associates domain
     Vector3D pnt(0,0,0);
     Vector3D s_norm(0,0,3);
     pd.set_domain( new PlanarDomain(s_norm, pnt) );

     double coords[] = { 1,.5, 0,
                         0, 0, 0,
                         1, 0, 0,
                         2, 0, 0,
                         2, 1, 0,
                         2, 2, 0,
                         1, 2, 0,
                         0, 2, 0,
                         0, 1, 0,
                         3, 0, 0,
                         3, 1, 0,
                         3, 2, 0 };

     size_t indices[] = { 1,  2,  0, 8,
                          2,  3,  4, 0,
                          8,  0,  6, 7,
                          0,  4,  5, 6,
                          3,  9, 10, 4,
                          4, 10, 11, 5 };
     
     bool fixed[] = { false, true, true, true, 
                      false, true, true, true, 
                      true, true, true, true };
     
     pd.fill( 12, coords, 6, QUADRILATERAL, indices, fixed, err );
   }

  // Create patch containing one ideal element, optionally higher-order.
  // For 2D elements, will attach appropriate planar domain.
  inline void create_ideal_element_patch( PatchData& pd,
                                          EntityTopology type, 
                                          size_t num_nodes,
                                          MsqError& err )
  {
     static PlanarDomain zplane(PlanarDomain::XY);
     static Settings settings;
     settings.set_slaved_ho_node_mode( Settings::SLAVE_NONE );
     pd.attach_settings( &settings );
     
  
       // build list of vertex coordinates
     const Vector3D* corners = unit_edge_element( type );
     std::vector<Vector3D> coords( corners, corners+TopologyInfo::corners(type) );
     bool mids[4] = {false};
     TopologyInfo::higher_order( type, num_nodes, mids[1], mids[2], mids[3], err );
     MSQ_ERRRTN(err);
     std::vector<size_t> conn(coords.size());
     for (unsigned i = 0; i < coords.size(); ++i)
       conn[i] = i;
 
     for (unsigned dim = 1; dim <= TopologyInfo::dimension(type); ++dim) {
       if (!mids[dim])
         continue;
       
       int num_side;
       if (dim == TopologyInfo::dimension(type))
         num_side = 1;
       else
         num_side = TopologyInfo::adjacent( type, dim );
       
       for (int s = 0; s < num_side; ++s) {
         unsigned idx = TopologyInfo::higher_order_from_side( type, num_nodes, dim, s, err );
         MSQ_ERRRTN(err);
         conn.push_back(idx);
       
         unsigned n;
         const unsigned* side = TopologyInfo::side_vertices( type, dim, s, n, err );
         MSQ_ERRRTN(err);
         Vector3D avg = coords[side[0]];
         for (unsigned v = 1; v < n; ++v)
           avg += coords[side[v]];
         avg *= 1.0/n;
         coords.push_back(avg);
       }
     }
     
     bool* fixed = new bool[coords.size()];
     std::fill( fixed, fixed+coords.size(), false );
     pd.fill( coords.size(), coords[0].to_array(), 1, &type, 
              &num_nodes, &conn[0], fixed, err );
     delete [] fixed;
     MSQ_ERRRTN(err);
     
     if (TopologyInfo::dimension(type) == 2)
       pd.set_domain( &zplane );
  }

     /*! \fn create_one_tri_patch(PatchData &one_tri_patch, MsqError &err)
            2
           / \      creates a Patch containing an ideal triangle
          /   \
         0-----1
         This Patch also has the normal information. 
     */
   inline void create_one_tri_patch(PatchData &one_tri_patch, MsqError &err)
   {
       /* ************** Creates normal info ******************* */
     Vector3D pnt(0,0,0);
     Vector3D s_norm(0,0,3);
     one_tri_patch.set_domain( new PlanarDomain( s_norm, pnt ) );

       /* *********************FILL tri************************* */
     double coords[] = { 1, 1, 1,
                         2, 1, 1,
                         1.5, 1+sqrt(3.0)/2.0, 1 };
     
     size_t indices[3] = { 0, 1, 2 };
     one_tri_patch.fill( 3, coords, 1, TRIANGLE, indices, 0, err );
   }

   /* Patch used in several quality metric tests.
      Our triangular patch is made of two tris.  tri_1 is a perfect
      equilateral (the ideal for most metrics).  tri_2 is an arbitrary
      triangle.
      Memory allocated in this function must be deallocated with
      destroy_patch_with_domain().
   */
   inline void create_qm_two_tri_patch_with_domain(PatchData &triPatch, MsqError &err) 
   {
     Vector3D pnt(0,0,0);
     Vector3D s_norm(0,0,3);
     triPatch.set_domain( new PlanarDomain(s_norm, pnt) );

     double coords[] = { 0.0, 0.0, 0.0,
                         1.0, 0.0, 0.0,
                         0.5, sqrt(3.0)/2.0, 0.0,
                         2.0, -4.0, 2.0 };
     

     const size_t conn[] = { 0, 1, 2, 0, 3, 1 };
     
     triPatch.fill( 4, coords, 2, TRIANGLE, conn, 0, err );
   }

     /* Patch used in several quality metric tests.
       Our quad patch is made of two quads.  quad_1 is a perfect
       square (the ideal for most metrics).  quad_2 is an arbitrary
       quad.
       Memory allocated in this function must be deallocated with
       destroy_patch_with_domain().
     */
   inline void create_qm_two_quad_patch_with_domain(PatchData &quadPatch, MsqError &err)
   {
     Vector3D pnt(0,0,0);
     Vector3D s_norm(0,0,3);
     quadPatch.set_domain( new PlanarDomain(s_norm, pnt) );

     double coords[] = { 0.0, 0.0, 0.0,
                         1.0, 0.0, 0.0,
                         1.0, 1.0, 0.0,
                         0.0, 1.0, 0.0,
                         2.0, -1.0, .5,
                         1.5, 1.0, 1.0 };
     
     const size_t conn[] = { 0, 1, 2, 3, 1, 4, 5, 2 };
     
     quadPatch.fill( 6, coords, 2, QUADRILATERAL, conn, 0, err );
   }

   /*! \fn create_two_tri_patch(PatchData &one_tri_patch, MsqError &err)
            2
           / \      creates a Patch containing two ideal triangles
          / 0 \
         0-----1
          \ 1 /
           \ /
            3
         This Patch also has the normal information. 
   */
   inline void create_two_tri_patch(PatchData &pd, MsqError &err)
   {
       /* ************** Creates normal info ******************* */
     Vector3D pnt(0,0,1);
     Vector3D s_norm(0,0,3);
     pd.set_domain( new PlanarDomain(s_norm, pnt) );

       // **********************FILL tri*************************

     double coords[] = { 1, 1, 1,
                         2, 1, 1,
                         1.5, 1+sqrt(3.0)/2.0, 1,
                         1.5, 1-sqrt(3.0)/2.0, 1 };

     size_t indices[] = { 0, 1, 2, 0, 3, 1 };
     
     pd.fill( 4, coords, 2, TRIANGLE, indices, 0, err );
  }

bool ObjectiveFunctionTemplate::initialize_block_coordinate_descent( Mesh* mesh, 
                                                      MeshDomain* domain, 
                                                      const Settings* settings,
                                                      PatchSet* ,
                                                      MsqError& err )
{
  std::auto_ptr<PatchSet> patch_set;
  switch (get_quality_metric()->get_metric_type())
  {
    case QualityMetric::VERTEX_BASED:  
      patch_set = std::auto_ptr<PatchSet>(new VertexPatches( 1, false ));
      break;
    case QualityMetric::ELEMENT_BASED: 
      patch_set = std::auto_ptr<PatchSet>(new ElementPatches);
      break;
    default: 
      MSQ_SETERR(err)("Cannot initialize for BCD for unknown metric type", 
                      MsqError::INVALID_STATE);
      return false;
  }

  clear();
  patch_set->set_mesh( mesh );
  PatchIterator patches( patch_set.get() );
  
  PatchData pd;
  pd.set_mesh( mesh );
  pd.set_domain( domain );
  if (settings)
    pd.attach_settings( settings );
  
  bool result = true;
  while (patches.get_next_patch( pd, err ) && !MSQ_CHKERR(err))
  {
    double value;
    bool b = evaluate( ObjectiveFunction::ACCUMULATE, pd, value, false, err ); 
    MSQ_ERRZERO(err);
    result = result && b;
  }
  return result;
}

void PatchDataTestNormals::setUp()
{
  MsqPrintError err(cout);
  
    // Define a mesh on the unit sphere
    // Make six quads corresponding to the six faces
    // of a cube inscribed in the sphere
  const double T = 1.0 / sqrt(3.0);
  double ucoords[] = {  T, -T, -T,
                        T,  T, -T,
                       -T,  T, -T,
                       -T, -T, -T,
                        T, -T,  T,
                        T,  T,  T,
                       -T,  T,  T,
                       -T, -T,  T };
  size_t uconn[] = { 3, 2, 1, 0, // -Z face
                     4, 5, 6, 7, // +Z face
                     0, 1, 5, 4, // +X face
                     1, 2, 6, 5, // +Y face
                     2, 3, 7, 6, // -X face
                     3, 0, 4, 7  // -Y face
                   };
  unboundedMesh.fill( 8, ucoords, 6, QUADRILATERAL, uconn, 0, err );
  CPPUNIT_ASSERT( !err );
  unboundedMesh.set_domain( &unboundedDomain );
  
    // Define a mesh on a cylinder with a radius of
    // one that is capped at z = +/- 2.  Define the
    // mesh as the 8 quads defining the sides of a pair of cubes
    // stacked axially in the cylinder
  const double V = 1.0 / sqrt(2.0);
  double bcoords[] = {  V, -V, -2,
                        V,  V, -2,
                       -V,  V, -2,
                       -V, -V, -2,
                       
                        V, -V,  0,
                        V,  V,  0,
                       -V,  V,  0,
                       -V, -V,  0,
                       
                        V, -V,  2,
                        V,  V,  2,
                       -V,  V,  2,
                       -V, -V,  2};
  size_t bconn[] = { // lower cube side faces
                      0,  1,  5,  4, // +X face
                      1,  2,  6,  5, // +Y face
                      2,  3,  7,  6, // -X face
                      3,  0,  4,  7, // -Y face
                     // upper cube side faces
                      4,  5,  9,  8, // +X face
                      5,  6, 10,  9, // +Y face
                      6,  7, 11, 10, // -X face
                      7,  4,  8, 11, // -Y face
                    };
  boundedMesh.fill( 12, bcoords, 8, QUADRILATERAL, bconn, 0, err );
  CPPUNIT_ASSERT( !err );
  boundedMesh.set_domain( &boundedDomain );
  
    // set element and vertex handles arrays
  size_t i = 0;
  for (i = 0; i < 12; ++i)
    boundedMesh.get_vertex_handles_array()[i] = (Mesh::VertexHandle)i;
  for (i = 0; i < 8; ++i)
    boundedMesh.get_element_handles_array()[i] = (Mesh::ElementHandle)i;
  
    // Bound the unit cylinder at +/- 1 on the z axis  
  std::vector<Mesh::VertexHandle> upper_curve(4), lower_curve(4);
  for (i = 0; i < 4; ++i)
  {
    lower_curve[i] = (Mesh::VertexHandle)i;
    upper_curve[i] = (Mesh::VertexHandle)(i+8);
  }
  boundedDomain.create_curve( -2, lower_curve );
  boundedDomain.create_curve(  2, upper_curve );
}

/*!This function is activated when cullingGlobalPatch is true.  It supplies
  cull_vertices with a single vertex-based patch at a time.  If the patch
  satisfies the culling criterion, it's free vertices are then soft-fixed.
 */
bool TerminationCriterion::cull_vertices_global(PatchData &global_patch,
                                                Mesh *mesh, MeshDomain *domain, const Settings *settings,
                                                OFEvaluator& of_eval,
                                                MsqError &err)
{
  if (!cullingGlobalPatch) return false;

    //PRINT_INFO("CULLING_METHOD FLAG = %i",cullingMethodFlag);
  
    //cull_bool will be changed to true if the criterion is satisfied
  bool cull_bool=false;

  std::vector<Mesh::VertexHandle> mesh_vertices;
  //std::vector<Mesh::VertexHandle> patch_vertices;
  //std::vector<Mesh::ElementHandle> patch_elements;
  //std::vector<Mesh::VertexHandle> fixed_vertices;
  //std::vector<Mesh::VertexHandle> free_vertices;

  // FIXME, verify global_patch is a global patch... how, is this right?
  mesh->get_all_vertices(mesh_vertices, err);
  size_t mesh_num_nodes = mesh_vertices.size();
  size_t global_patch_num_nodes = global_patch.num_nodes() ;
  if (0)  std::cout << "tmp srk mesh_num_nodes= " << mesh_num_nodes << " global_patch_num_nodes= " 
            << global_patch_num_nodes << std::endl;
  if (mesh_num_nodes != global_patch_num_nodes)
    {
      std::cout << "tmp srk cull_vertices_global found non global patch" << std::endl;
      exit(123);
      return false;
    }
  PatchData patch;
  patch.set_mesh( (Mesh*) mesh );
  patch.set_domain( domain );
  patch.attach_settings( settings );

  const MsqVertex* global_patch_vertex_array = global_patch.get_vertex_array( err );
  Mesh::VertexHandle* global_patch_vertex_handles = global_patch.get_vertex_handles_array();

  int num_culled = 0;
  for (unsigned iv=0; iv < global_patch_num_nodes; iv++)
    {
      // form a patch for this vertex; if it is culled, set it to be soft fixed
      Mesh::VertexHandle vert = global_patch_vertex_handles[iv];
      std::vector<Mesh::ElementHandle> elements;
      std::vector<size_t> offsets;
      mesh->vertices_get_attached_elements(&vert, 1, elements, offsets, err);
      
      std::set<Mesh::VertexHandle> patch_free_vertices_set;

      for (unsigned ie=0; ie < elements.size(); ie++)
        {
          std::vector<Mesh::VertexHandle> vert_handles;
          std::vector<size_t> v_offsets;
          mesh->elements_get_attached_vertices(&elements[ie], 1, vert_handles, v_offsets, err);
          for (unsigned jv=0; jv < vert_handles.size(); jv++)
            {
              unsigned char bt;
              mesh->vertex_get_byte(vert_handles[jv], &bt, err);
              MsqVertex v;
              v.set_flags(bt);
              if (v.is_free_vertex())
                patch_free_vertices_set.insert(vert_handles[jv]);
            }
        }

      std::vector<Mesh::VertexHandle> patch_free_vertices_vector(patch_free_vertices_set.begin(), patch_free_vertices_set.end());
      //std::vector<unsigned char> byte_vector(patch_vertices_vector.size());
      //mesh->vertices_get_byte(&vert_handles[0], &byte_vector[0], vert_handles.size(), err);

      patch.set_mesh_entities( elements, patch_free_vertices_vector, err );
      if (cull_vertices(patch, of_eval, err))
        {
          //std::cout << "tmp srk cull_vertices_global found culled patch" << std::endl;
          Mesh::VertexHandle* patch_vertex_handles = patch.get_vertex_handles_array();
          const MsqVertex* patch_vertex_array = patch.get_vertex_array( err );
          for (unsigned jv=0; jv < patch.num_nodes(); jv++)
            {
              if (patch_vertex_handles[jv] == global_patch_vertex_handles[iv])
                {
                  if (patch_vertex_array[jv].is_flag_set(MsqVertex::MSQ_CULLED))
                    {
                      global_patch.set_vertex_culled(iv);
                      ++num_culled;
                      cull_bool = true;
                      //std::cout << "tmp srk cull_vertices_global found culled vertex" << std::endl;
                    }
                }
            }
        }
    }
    if (0)  std::cout << "tmp srk cull_vertices_global found " << num_culled << " culled vertices out of " << global_patch_num_nodes << std::endl;

  return cull_bool;
}

/*! \brief Improves the quality of the MeshSet, calling some
    methods specified in a class derived from VertexMover

    \param const MeshSet &: this MeshSet is looped over. Only the
    mutable data members are changed (such as currentVertexInd).
  */
double VertexMover::loop_over_mesh( ParallelMesh* mesh,
                                    MeshDomain* domain,
                                    const Settings* settings,
                                    MsqError& err )
{
  std::vector<size_t> junk;
  Mesh::VertexHandle vertex_handle;

    // Get the patch data to use for the first iteration
  OFEvaluator& obj_func = get_objective_function_evaluator();
  
  PatchData patch;
  patch.set_mesh( (Mesh*) mesh );
  patch.set_domain( domain );
  patch.attach_settings( settings );

  ParallelHelper* helper = mesh->get_parallel_helper();
  if (!helper) {
    MSQ_SETERR(err)("No ParallelHelper instance", MsqError::INVALID_STATE);
    return 0;
  }

  helper->smoothing_init(err);  MSQ_ERRZERO(err);

  bool did_some, all_culled;
  std::vector<Mesh::VertexHandle> patch_vertices;
  std::vector<Mesh::ElementHandle> patch_elements;
  std::vector<Mesh::VertexHandle> fixed_vertices;
  std::vector<Mesh::VertexHandle> free_vertices;
   
    // Get termination criteria
  TerminationCriterion* outer_crit=this->get_outer_termination_criterion();
  TerminationCriterion* inner_crit=this->get_inner_termination_criterion();
  if(outer_crit == 0){
    MSQ_SETERR(err)("Termination Criterion pointer is Null", MsqError::INVALID_STATE);
    return 0.;
  }
  if(inner_crit == 0){
    MSQ_SETERR(err)("Termination Criterion pointer for inner loop is Null", MsqError::INVALID_STATE);
    return 0.;
  }
  
  PatchSet* patch_set = get_patch_set();
  if (!patch_set) {
    MSQ_SETERR(err)("No PatchSet for QualityImprover!", MsqError::INVALID_STATE);
    return 0.0;
  }
  patch_set->set_mesh( (Mesh*)mesh );

  std::vector<PatchSet::PatchHandle> patch_list;
  patch_set->get_patch_handles( patch_list, err ); MSQ_ERRZERO(err);
  
    // Initialize outer loop
    
  this->initialize(patch, err);        
  if (MSQ_CHKERR(err)) goto ERROR;
  
  obj_func.initialize( (Mesh*)mesh, domain, settings, patch_set, err ); 
  if (MSQ_CHKERR(err)) goto ERROR;
  
  outer_crit->reset_outer( (Mesh*)mesh, domain, obj_func, settings, err); 
  if (MSQ_CHKERR(err)) goto ERROR;
   
   // Loop until outer termination criterion is met
  did_some = true;
  while (!outer_crit->terminate())
  {
    if (!did_some) {
      MSQ_SETERR(err)("Inner termiation criterion satisfied for all patches "
                      "without meeting outer termination criterion.  This is "
                      "an infinite loop.  Aborting.", MsqError::INVALID_STATE);
      break;
    }
    did_some = false;
    all_culled = true;

    ///*** smooth the interior ***////

    // get the fixed vertices (i.e. the ones *not* part of the first independent set)
    helper->compute_first_independent_set(fixed_vertices); 

    // sort the fixed vertices
    std::sort(fixed_vertices.begin(), fixed_vertices.end());

      // Loop over each patch
    std::vector<PatchSet::PatchHandle>::iterator p_iter = patch_list.begin();
    while( p_iter != patch_list.end() )
    {
      // loop until we get a non-empty patch.  patch will be empty
      // for culled vertices with element-on-vertex patches
      do {
	patch_set->get_patch( *p_iter, patch_elements, patch_vertices, err );
	if (MSQ_CHKERR(err)) goto ERROR;
	++p_iter;
//.........这里部分代码省略.........

template <class QMType> inline
void TMPQualityMetricTest<QMType>::test_gradient_common(TargetCalculator* tc)
{
  MsqPrintError err(std::cout);
  
    // check for expected value at center of flattened hex
  
    // construct flattened hex
  const double y = 0.5;
  const double vertices[] = { 0.0, 0.0, 0.0,
                              1.0, 0.0, 0.0,
                              1.0, y  , 0.0,
                              0.0, y  , 0.0 };
  size_t conn[8] = { 0, 1, 2, 3 };
  PatchData pd;
  pd.fill( 4, vertices, 1, QUADRILATERAL, conn, 0, err );
  ASSERT_NO_ERROR(err);
  
    // calculate Jacobian matrix at element center
    
    // derivatives of bilinear map at quad center
  const double deriv_vals[] = { -0.5, -0.5,
                                 0.5, -0.5,
                                 0.5,  0.5,
                                -0.5,  0.5 } ;
  MsqMatrix<4,2> coeff_derivs(deriv_vals);
  MsqMatrix<4,3> coords( vertices );
  MsqMatrix<3,2> J = transpose(coords) * coeff_derivs;
    // calculate expected metric value
  const double expt_val = sqr_Frobenius( J );
    // calculate derivative for each element vertex
  MsqVector<3> expt_grad[4];
  for (int v = 0; v < 4; ++v)
    expt_grad[v] = 2 * J * transpose( coeff_derivs.row(v) );
    
  
    // construct metric
  pd.attach_settings( &settings );
  TestGradTargetMetric< typename TMPTypes<QMType>::MetricType > tm;
  //IdealShapeTarget tc;
  QMType m( tc, &tm );
  PlanarDomain plane( PlanarDomain::XY );
  pd.set_domain( &plane );
  
    // evaluate metric
  double act_val;
  std::vector<size_t> indices;
  std::vector<Vector3D> act_grad;
  size_t h = ElemSampleQM::handle( Sample(2, 0), 0 );
  m.evaluate_with_gradient( pd, h, act_val, indices, act_grad, err );
  ASSERT_NO_ERROR(err);
  
    // compare values
  CPPUNIT_ASSERT_DOUBLES_EQUAL( expt_val, act_val, 1e-10 );
  CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(expt_grad[indices[0]].data()), act_grad[0], 1e-10 );
  CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(expt_grad[indices[1]].data()), act_grad[1], 1e-10 );
  CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(expt_grad[indices[2]].data()), act_grad[2], 1e-10 );
  CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(expt_grad[indices[3]].data()), act_grad[3], 1e-10 );

    // check numerical approx of gradient
  m.QualityMetric::evaluate_with_gradient( pd, h, act_val, indices, act_grad, err );
  ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT_DOUBLES_EQUAL( expt_val, act_val, 1e-10 );
  CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(expt_grad[indices[0]].data()), act_grad[0], 1e-5 );
  CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(expt_grad[indices[1]].data()), act_grad[1], 1e-5 );
  CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(expt_grad[indices[2]].data()), act_grad[2], 1e-5 );
  CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(expt_grad[indices[3]].data()), act_grad[3], 1e-5 );
}

double TargetWriter::loop_over_mesh( MeshDomainAssoc* mesh_and_domain,
                                     const Settings* settings,
                                     MsqError& err )
{
  Mesh* mesh = mesh_and_domain->get_mesh();
  MeshDomain* domain = mesh_and_domain->get_domain();
 
  PatchData patch;
  patch.set_mesh( mesh );
  patch.set_domain( domain );
  if (settings)
    patch.attach_settings( settings );
  
  ElementPatches patch_set;
  patch_set.set_mesh( mesh );
  std::vector<PatchSet::PatchHandle> patches;
  std::vector<PatchSet::PatchHandle>::iterator p;
  std::vector<Mesh::VertexHandle> patch_verts;
  std::vector<Mesh::ElementHandle> patch_elems;
  
  patch_set.get_patch_handles( patches, err ); MSQ_ERRZERO(err);
  
  std::vector< MsqMatrix<3,3> > targets3d;
  std::vector< MsqMatrix<3,2> > targets2dorient;
  std::vector< MsqMatrix<2,2> > targets2d;
  std::vector< double > weights;
  std::vector< Sample > samples;
  for (p = patches.begin(); p != patches.end(); ++p)
  {
    patch_verts.clear();
    patch_elems.clear();
    patch_set.get_patch( *p, patch_elems, patch_verts, err ); MSQ_ERRZERO(err);
    patch.set_mesh_entities( patch_elems, patch_verts, err ); MSQ_ERRZERO(err);
    assert(patch.num_elements() == 1);
    
    MsqMeshEntity& elem = patch.element_by_index(0);
    EntityTopology type = elem.get_element_type();
    patch.get_samples( 0, samples, err ); MSQ_ERRZERO(err);
    if (samples.empty())
      continue;    
    
    if (targetCalc) {
      const unsigned dim = TopologyInfo::dimension(type);
      if (dim == 3) {
        targets3d.resize( samples.size() );
        for (unsigned i = 0; i < samples.size(); ++i) {
          targetCalc->get_3D_target( patch, 0, samples[i], targets3d[i], err ); MSQ_ERRZERO(err);

          if (DBL_EPSILON > det(targets3d[i])) {
            MSQ_SETERR(err)("Inverted 3D target", MsqError::INVALID_ARG);
            return 0.0;
          }
        }
        
        TagHandle tag = get_target_tag( 3, samples.size(), mesh,  err ); MSQ_ERRZERO(err);
        mesh->tag_set_element_data( tag, 1, 
                                    patch.get_element_handles_array(), 
                                    arrptr(targets3d), err ); MSQ_ERRZERO(err);
      }
      else if(targetCalc->have_surface_orient()) {
        targets2dorient.resize( samples.size() );
        for (unsigned i = 0; i < samples.size(); ++i) {
          targetCalc->get_surface_target( patch, 0, samples[i], targets2dorient[i], err ); MSQ_ERRZERO(err);

          MsqMatrix<3,1> cross = targets2dorient[i].column(0) * targets2dorient[i].column(1);
          if (DBL_EPSILON > (cross%cross)) {
            MSQ_SETERR(err)("Degenerate 2D target", MsqError::INVALID_ARG);
            return 0.0;
          }
        }
        
        TagHandle tag = get_target_tag( 2, samples.size(), mesh, err ); MSQ_ERRZERO(err);
        mesh->tag_set_element_data( tag, 1, 
                                    patch.get_element_handles_array(), 
                                    arrptr(targets2dorient), err ); MSQ_ERRZERO(err);
      }
      else {
        targets2d.resize( samples.size() );
        for (unsigned i = 0; i < samples.size(); ++i) {
          targetCalc->get_2D_target( patch, 0, samples[i], targets2d[i], err ); MSQ_ERRZERO(err);

          if (DBL_EPSILON > det(targets2d[i])) {
            MSQ_SETERR(err)("Degenerate/Inverted 2D target", MsqError::INVALID_ARG);
            return 0.0;
          }
        }
        
        TagHandle tag = get_target_tag( 2, samples.size(), mesh, err ); MSQ_ERRZERO(err);
        mesh->tag_set_element_data( tag, 1, 
                                    patch.get_element_handles_array(), 
                                    arrptr(targets2d), err ); MSQ_ERRZERO(err);
      }
    }
      
    if (weightCalc) {
      weights.resize( samples.size() );
      for (unsigned i = 0; i < samples.size(); ++i) {
        weights[i] = weightCalc->get_weight( patch, 0, samples[i], err ); MSQ_ERRZERO(err);
      }
      TagHandle tag = get_weight_tag( samples.size(), mesh, err ); MSQ_ERRZERO(err);
//.........这里部分代码省略.........