C++ MeshImpl::write_vtk方法代码示例

// 输出网格到 VTK 文件
void MeshOpt::writeToVTK(hier::Patch<NDIM>& patch,
                         const double  time,
                         const double  dt,
                         const bool    initial_time)
{
    NULL_USE(dt);
    NULL_USE(time);
    NULL_USE(initial_time);

    const tbox::Pointer< hier::BlockPatchGeometry<NDIM> > pgeom =
            patch.getPatchGeometry();

    int block_index = pgeom->getBlockNumber();

    int patch_index = patch.getPatchNumber();

    std::stringstream bi, pi, df;
    bi << block_index;
    pi << patch_index;
    df << d_flag;

    std::string file_name = df.str() + "_block_ " + bi.str()+ "_patch_" +  pi.str()  + ".vtk";


    MsqError err;
    MeshImpl * mesh = createLocalMesh(patch);
    mesh->write_vtk(file_name.c_str(), err);

    return;
}

   // Test writing VTK unstructured mesh
 void VtkTest::test_write()
 {
   MeshImpl mesh;
   MsqPrintError err(cout);
   
     // Create file containing unstructured mesh test case
   FILE* file = fopen( temp_file_name, "w+" );
   CPPUNIT_ASSERT(file);
   int rval = fputs( mixed_unstructured_data, file );
   fclose( file );
   if (rval == EOF) remove(temp_file_name);
   CPPUNIT_ASSERT(rval != EOF);
   
     // Read unstructured mesh file
   mesh.read_vtk( temp_file_name, err );
   remove(temp_file_name);
   ASSERT_NO_ERROR(err);
   
     // Write unstructured mesh file back out
   mesh.write_vtk( temp_file_name, err );
   if (err)  remove( temp_file_name );
   ASSERT_NO_ERROR(err);
   
     // Check if file contained expected mesh
   test_read_unstructured( temp_file_name );
   remove( temp_file_name );
 }

int main( int argc, char* argv[] )
{
  parse_options( argv, argc );
  
  MeshImpl mesh;
  XYRectangle domain( max_x - min_x, max_y - min_y, min_x, min_y );
  MsqError err;
  
  create_input_mesh( input_x, mesh, err );
  if (MSQ_CHKERR(err)) { std::cerr << err << std::endl; return 2; }
  
  domain.setup( &mesh, err );
  if (MSQ_CHKERR(err)) { std::cerr << err << std::endl; return 2; }
  
  QualityMetric* metric = 0;
  if (mMetric == 'c')
    metric = new ConditionNumberQualityMetric;
  else
    metric = new IdealWeightInverseMeanRatio;
  
  LPtoPTemplate function( 1, metric );
  
  VertexMover* solver = 0;
  if (mSolver == 'j')
    solver = new ConjugateGradient( &function );
  else
    solver = new FeasibleNewton( &function );
    
  if (PatchSetUser* psu = dynamic_cast<PatchSetUser*>(solver)) 
    psu->use_global_patch();
  
  TerminationCriterion inner;
  inner.add_absolute_vertex_movement( 1e-4 );
  inner.write_mesh_steps( "synchronous", TerminationCriterion::GNUPLOT );
  solver->set_inner_termination_criterion( &inner );  
  
  InstructionQueue q;
  QualityAssessor qa( metric, 10 );
  q.add_quality_assessor( &qa, err );
  q.set_master_quality_improver( solver, err );
  q.add_quality_assessor( &qa, err );
  MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &domain);
  q.run_instructions( &mesh_and_domain, err );
  delete solver;
  delete metric;
  
  if (MSQ_CHKERR(err)) 
    { std::cerr << err << std::endl; return 3; }
    
  mesh.write_vtk( outputFile, err );
  if (MSQ_CHKERR(err)) 
    { std::cerr << err << std::endl; return 2; }
  
  return 0;
}

void cond_write_file( MeshImpl& mesh, const char* filename ) 
{
  if (filename) {
    MsqPrintError err(std::cerr);
    mesh.write_vtk( filename, err );
    if (MSQ_CHKERR(err)) {
      std::cerr << filename << ": failed to write file" << std::endl;
      exit(1);
    }
    std::cout << "Wrote file: " << filename << std::endl;
  }
}

int run_example( const Example& e, bool write_output_file )
{
  MsqPrintError err(std::cerr);
  MeshImpl mesh;
  DomainClassifier domain;
  HexLagrangeShape hex27;
  
  std::cout << std::endl 
            << "--------------------------------------------------------------------"
            << std::endl
            << e.desc << std::endl
            << "--------------------------------------------------------------------"
            << std::endl;
   
  std::string name = e.func( domain, mesh, err );
  if (MSQ_CHKERR(err)) return 2;
  std::cout << "Loaded mesh from: " << name << std::endl;
  
  UntangleWrapper untangler;
  untangler.set_slaved_ho_node_mode( Settings::SLAVE_NONE );
  untangler.set_mapping_function( &hex27 );
  MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &domain, false, true);
  untangler.run_instructions( &mesh_and_domain, err ); 
  if (MSQ_CHKERR(err)) return 1;
  ShapeImprover smoother;
  smoother.set_slaved_ho_node_mode( Settings::SLAVE_NONE );
  smoother.set_mapping_function( &hex27 );
  smoother.set_vertex_movement_limit_factor( 0.05 );
  smoother.run_instructions( &mesh_and_domain, err );
  if (MSQ_CHKERR(err)) return 1;
  
  if (write_output_file) {
    size_t idx = name.find( ".vtk" );
    if (idx != std::string::npos) {
      std::string newname( name.substr(0, idx) );
      newname += ".out";
      newname += name.substr(idx);
      name.swap(newname);
    }
    else {
      name += ".out";
    }
    
    mesh.write_vtk( name.c_str(), err ); MSQ_CHKERR(err);
    std::cout << "Write mesh to file: " << name << std::endl;
  }
  
  return smoother.quality_assessor().invalid_elements();
}

   // Test writing quadtratic elements
 void VtkTest::test_write_field_attrib()
 {
   MeshImpl mesh;
   MsqPrintError err(cout);
   
     // Create file containing unstructured mesh test case
   FILE* file = fopen( temp_file_name, "w+" );
   fputs( structured_3d_points_data, file );
   fputs( simple_field_attrib, file );
   fclose( file );
   
     // Read unstructured mesh file
   mesh.read_vtk( temp_file_name, err );
   remove(temp_file_name);
   ASSERT_NO_ERROR(err);
   
     // Write unstructured mesh file back out
   mesh.write_vtk( temp_file_name, err );
   if (err)  remove( temp_file_name );
   ASSERT_NO_ERROR(err);
   
     // Check if file contained expected mesh
   check_field_attrib( temp_file_name );
 }

static int do_smoother( const char* input_file, 
                        const char* output_file, 
                        const char* ref_mesh_file,
                        double of_power, 
                        unsigned metric_idx,
                        AveragingScheme avg_scheme )
{
  MsqPrintError err(cerr);
  
  TMetric *const target_metric = metrics[metric_idx].u;
  cout << "Input file:  " << input_file << endl;
  cout << "Metric:      ";
  if (avg_scheme != NONE)
    cout << averaging_names[avg_scheme] << " average of ";
  cout << metrics[metric_idx].n << endl;
  cout << "Of Power:    " << of_power << endl;
  
  
  auto_ptr<TargetCalculator> tc;
  auto_ptr<MeshImpl> ref_mesh_impl;
  auto_ptr<ReferenceMesh> ref_mesh;
  if (ref_mesh_file) {
    ref_mesh_impl.reset(new MeshImpl);
    ref_mesh_impl->read_vtk( ref_mesh_file, err );
    if (MSQ_CHKERR(err)) return 2;
    ref_mesh.reset( new ReferenceMesh( ref_mesh_impl.get() ));
    tc.reset( new RefMeshTargetCalculator( ref_mesh.get() ) );
  }
  else {
    tc.reset( new IdealShapeTarget( ) );
  }
    
  TQualityMetric jacobian_metric( tc.get(), target_metric );
  ElementPMeanP elem_avg( of_power, &jacobian_metric );
  VertexPMeanP vtx_avg( of_power, &jacobian_metric );
  QualityMetric* mmetrics[] = { &jacobian_metric, &elem_avg, &vtx_avg, &jacobian_metric };
  QualityMetric* metric = mmetrics[avg_scheme];

  TerminationCriterion outer, inner;
  outer.add_iteration_limit( 1 );
  inner.add_absolute_vertex_movement( 1e-4 );
  inner.add_iteration_limit( 100 );
  
  PMeanPTemplate obj1( of_power, metric );
  PatchPowerMeanP obj2( of_power, metric );
  ObjectiveFunction& objective = *((avg_scheme == PATCH) ? (ObjectiveFunction*)&obj2 : (ObjectiveFunction*)&obj1);
  
  ConjugateGradient solver( &objective, err );
  if (MSQ_CHKERR(err)) return 1;
  solver.set_inner_termination_criterion( &inner );
  solver.set_outer_termination_criterion( &outer );
  solver.use_global_patch();
  
  ConditionNumberQualityMetric qm_metric;
  QualityAssessor before_assessor;
  QualityAssessor after_assessor;
  before_assessor.add_quality_assessment( metric, 10);
  before_assessor.add_quality_assessment( &qm_metric );
  after_assessor.add_quality_assessment( metric, 10 );

  InstructionQueue q;
  q.add_quality_assessor( &before_assessor, err );
  q.set_master_quality_improver( &solver, err );
  q.add_quality_assessor( &after_assessor, err );
  
  MeshImpl mesh;
  mesh.read_vtk( input_file, err );
  if (MSQ_CHKERR(err)) return 2;
  PlanarDomain geom = make_domain( &mesh, err );
  if (MSQ_CHKERR(err)) return 1;
  q.run_instructions( &mesh, &geom, err );
  if (MSQ_CHKERR(err)) return 3;
  mesh.write_vtk( output_file, err );
  if (MSQ_CHKERR(err)) return 2;
  cout << "Wrote: " << output_file << endl;

  before_assessor.scale_histograms(&after_assessor);
  
  return 0;
}

//.........这里部分代码省略.........
    EntityTopology type;
    mesh.elements_get_topologies( &elem, &type, 1, err );
    std::copy( conn.begin() + TopologyInfo::corners(type), conn.end(), 
               std::back_inserter( higher_order ) );
  }
  std::sort( higher_order.begin(), higher_order.end() );
  higher_order.erase( std::unique( higher_order.begin(), higher_order.end() ), 
                      higher_order.end() );

    // build depth map for all vertices
  while (!next.empty()) {
    curr.swap( next );
    next.clear();
    while (!curr.empty()) {
      Mesh::ElementHandle elem = curr.back();
      curr.pop_back();
      
      conn.clear();
      mesh.elements_get_attached_vertices( &elem, 1, conn, off, err );
      if (err) return 1;
      
      int min = std::numeric_limits<int>::max();
      for (size_t i = 0; i < conn.size(); ++i) {
        d_iter = depth.find( conn[i] );
        if (d_iter != depth.end() && d_iter->second < min)
          min = d_iter->second;
      }
      
      if (min == std::numeric_limits<int>::max()) {
        next.push_back( elem );
        continue;
      }
      
      for (size_t i = 0; i < conn.size(); ++i) {
        d_iter = depth.find( conn[i] );
      
        if (d_iter == depth.end() || d_iter->second > min+1)
          depth[conn[i]] = min+1;
      }
    }
  }
  
    // write depth map to tag for debugging purposes
  std::vector<int> depth_vals(verts.size());
  for (size_t i = 0; i < verts.size(); ++i)
    depth_vals[i] = depth[verts[i]];
  TagHandle tag = mesh.tag_create( "depth", Mesh::INT, 1, 0, err );
  if (err) return 1;
  mesh.tag_set_vertex_data( tag, verts.size(), arrptr(verts), arrptr(depth_vals), err );
  if (err) return 1;
  
  
    // set tag specifying slaved vertices
  for (size_t i = 0; i < verts.size(); ++i)
    if (std::binary_search( higher_order.begin(), higher_order.end(), verts[i] ))
      depth_vals[i] = depth[verts[i]] > n;
    else
      depth_vals[i] = 0;
  tag = mesh.tag_create( "slaved", Mesh::INT, 1, 0, err );
  if (err) return 1;
  mesh.tag_set_vertex_data( tag, verts.size(), arrptr(verts), arrptr(depth_vals), err );
  if (err) return 1;
  
    // perturb mid-edge nodes along boundary
  std::vector<MsqVertex> coords;
  for (size_t i = 0; i < skin.size(); ++i) {
    if (!std::binary_search( higher_order.begin(), higher_order.end(), skin[i]))
      continue;
  
    curr.clear();
    mesh.vertices_get_attached_elements( &skin[i], 1, curr, off, err );
    if (err) return 1;
    assert(curr.size() == 1);
    conn.clear();
    mesh.elements_get_attached_vertices( arrptr(curr), 1, conn, off, err );
    if (err) return 1;
    
    // estimate element center
    coords.resize( conn.size() );
    mesh.vertices_get_coordinates( arrptr(conn), arrptr(coords), conn.size(), err );
    if (err) return 1;
    
    Vector3D mean(0.0);
    for (size_t j = 0; j < coords.size(); ++j)
      mean += coords[j];
    mean /= coords.size();
    
    size_t idx = std::find( conn.begin(), conn.end(), skin[i] ) - conn.begin();
    assert(idx < conn.size());
    Vector3D init = coords[idx];
    Vector3D pos = (1 - PERTURB_FRACT) * init + PERTURB_FRACT * mean;
    mesh.vertex_set_coordinates( skin[i], pos, err );
    if (err) return 1;
  }
  
  mesh.write_vtk( argv[3], err );
  if (err) return 1;
  
  return 0;
}

int uwt( bool skip,
         UntangleWrapper::UntangleMetric metric,
         const char* input_file_base,
         int expected,
         bool flip_domain )
{
  if (skip)
    return 0;

  if (!brief_output)
    std::cout << std::endl << "**********************************************" << std::endl;
  std::cout << "Running \"" << input_file_base << "\" for " << tostr(metric) << std::endl;
  if (!brief_output)
    std::cout << "**********************************************" << std::endl << std::endl;

    // get mesh
  MsqError err;
  MeshImpl mesh;
  std::string input_file( VTK_2D_DIR );
  input_file += input_file_base;
  mesh.read_vtk( input_file.c_str(), err );
  if (err) {
    std::cerr << err << std::endl;
    std::cerr << "ERROR: " << input_file << " : failed to read file" << std::endl;
    return 1;
  }
    // get domain
  std::vector<Mesh::VertexHandle> verts;
  mesh.get_all_vertices( verts, err );
  if (err || verts.empty()) abort();
  MsqVertex coords;
  mesh.vertices_get_coordinates( arrptr(verts), &coords, 1, err );
  if (err) abort();
  Vector3D norm(0,0,flip_domain ? -1 : 1);
  PlanarDomain domain( norm, coords );
    // run wrapper
  UntangleWrapper wrapper( metric );
  wrapper.set_vertex_movement_limit_factor( 0.005 );
  double constant = (metric == UntangleWrapper::BETA) ? beta : mu_sigma;
  if (constant > 0)
    wrapper.set_metric_constant( constant );
  if (brief_output)
    wrapper.quality_assessor().disable_printing_results();
  MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &domain);
  wrapper.run_instructions( &mesh_and_domain, err );
  if (err) {
    std::cerr << err << std::endl;
    std::cerr << "ERROR: optimization failed" << std::endl;
    return 1;
  }
    // write output file
  if (write_output) {
    std::string result_file(tostr(metric));
    result_file += "-";
    result_file += input_file_base;
    mesh.write_vtk( result_file.c_str(), err );
    if (err) {
      std::cerr << err << std::endl;
      std::cerr << "ERROR: " << result_file << " : failed to write file" << std::endl;
      err.clear();
    }
    else {
      std::cerr << "Wrote file: " << result_file << std::endl;
    }
  }
  
    // test number of inverted elements
  int count, junk;
  wrapper.quality_assessor().get_inverted_element_count( count, junk, err );
  if (err) abort();
  if (count < expected) {
    std::cout << "WARNING: expected " << expected 
              << " inverted elements but finished with only " 
              << count << std::endl
              << "Test needs to be updated?" << std::endl << std::endl;
    return 0;
  }
  else if (count == expected) {
    std::cout << "Completed with " << count << " inverted elements remaining" 
              << std::endl << std::endl;
    return 0;
  }
  else {
    std::cerr << "ERROR: expected " << expected 
              << " inverted elements but finished with " 
              << count << std::endl << std::endl;
    return 1;
  }
}

static int test(std::string filename_prefix, std::string mesh_topology_name, MeshDomain *domain=0)
{
  int rank, nprocs;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);

  if (nprocs > 2) { cerr << "parallel_untangle_shape::test(" << mesh_topology_name << " can only be run with 1 or 2 processors" << std::endl; return 0; }

  /* create processor-specific file names */
  ostringstream in_name, out_name, gold_name;
  in_name << filename_prefix << "par_untangle_original_" << mesh_topology_name << "_mesh." << nprocs << "." << rank << ".vtk";
  gold_name << filename_prefix << "par_untangle_smoothed_" << mesh_topology_name << "_mesh." << nprocs << "." << rank << ".vtk";
  out_name << "par_untangle_smoothed_" << mesh_topology_name << "_mesh." << nprocs << "." << rank << ".vtk";

  cout << "in_name= " << in_name.str() << " gold_name= " << gold_name.str() << " out_name= " << out_name.str() << std::endl;

  /* load different mesh files on each processor */
  MsqError err;
  MeshImpl mesh;
  mesh.read_vtk(in_name.str().c_str(), err);
  if (err) {cerr << err << endl; return 1;}

  /* create parallel mesh instance, specifying tags 
   * containing parallel data */
  ParallelMeshImpl parallel_mesh(&mesh, "GLOBAL_ID", "PROCESSOR_ID");
  ParallelHelperImpl helper;
  helper.set_communicator(MPI_COMM_WORLD);
  helper.set_parallel_mesh(&parallel_mesh);
  parallel_mesh.set_parallel_helper(&helper);

  /* do Laplacian smooth */
  //LaplaceWrapper optimizer;
  //optimizer.run_instructions(&parallel_mesh, err);

  int  msq_debug             = 0; // 1,2,3 for more debug info
  bool always_smooth         = true;
  int innerIter = 100;
  double gradNorm = 1.e-9;

  ParShapeImprover si(innerIter, gradNorm);
  //Mesh *pmesh = &parallel_mesh;
  si.run(parallel_mesh, domain, err, always_smooth, msq_debug);
  if (err) {cerr << err << endl; return 1; }

  /* write mesh */
  mesh.write_vtk(out_name.str().c_str(),err);
  if (err) {cerr << err << endl; return 1;}

  //std::cout << "P[ " << rank <<"] reading gold..." << std::endl;

  /* compare mesh with gold copy */
  MeshImpl gold;
  gold.read_vtk(gold_name.str().c_str(),err);
  if (err) {cerr << err << endl; return 1;}

  //std::cout << "P[ " << rank <<"] read gold, checking mesh diff..." << std::endl;
  bool do_print=true;
  double tol = 1.e-4;
  bool diff = MeshUtil::meshes_are_different(mesh, gold, err, tol, do_print);
  if (err) {cerr << err << endl; return 1;}
  //std::cout << "P[ " << rank <<"] read gold, checking mesh diff...done" << std::endl;

  if (diff) {cerr << "Error, computed mesh different from gold copy" << std::endl; return 1;}
  
  print_timing_diagnostics(cout);

  return 0;
}

int main( int argc, char* argv[] )
{
  const char* input_file = 0;
  const char* output_file = 0;
  for (int i = 1; i < argc; ++i) {
    if (!strcmp("-h", argv[i]))
      usage(argv[0],true);
    else if (!input_file)
      input_file = argv[i];
    else if (!output_file)
      output_file = argv[i];
    else
      usage(argv[0]);
  }
  if (!input_file)
    input_file = DEFAULT_INPUT;
  
  MsqError err;
  MeshImpl mesh;
  mesh.read_vtk( input_file, err );
  if (err) {
    std::cerr << err << std::endl
              << input_file << ": failed to read file" << std::endl;
    return 3;
  }
  
  PlanarDomain plane(PlanarDomain::XY);
#ifdef TEST_OLD_WRAPPER
  ShapeImprovementWrapper smoother;
#else
  ShapeImprover smoother;
#endif
  IdealWeightInverseMeanRatio extra_metric;
  smoother.quality_assessor().add_quality_assessment(&extra_metric);
  MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &plane);
  smoother.run_instructions( &mesh_and_domain, err );
  if (err) {
    std::cerr << err << std::endl
              << input_file << ": smoother failed" << std::endl;
    return 2;
  }
  
  if (output_file) {
    mesh.write_vtk( output_file, err );
    if (err) {
      std::cerr << err << std::endl
                << output_file << ": failed to write file" << std::endl;
      return 3;
    }
  }

  if (smoother.quality_assessor().invalid_elements()) {
    std::cerr << "Resulting mesh contains invalid elements: untangler did not succeed" << std::endl;
    return 4;
  }
  
  const QualityAssessor::Assessor* quality = 
    smoother.quality_assessor().get_results(&extra_metric);
  if (!quality) {
    std::cerr << "Failed to get quality stats for IMR metric" << std::endl;
    return 2;
  }
  
  if (fabs(1 - quality->get_average()) > 1e-3) {
    std::cerr << "Average quality is not optimal." << std::endl;
    return 4;
  }
  
  if (quality->get_stddev() > 1e-3) {
    std::cerr << "Not all elements have optimal quality." << std::endl;
    return 4;
  }
  
  return 0;
}

int main( int argc, char* argv[] )
{
  MeshTransform xform;
  RotateArg rotate_arg( &xform );
  ScaleArg scale_arg( &xform );
  TranslateArg translate_arg( &xform );
  CLArgs::ToggleArg freeonly, skin;

  CLArgs args( "vtkxform", "Transform a mesh",
               "Apply one or more transformations to vertex coordinates "
               "in a mesh read from a VTK file." );
  
  const char* ROTATE_VALS[] = { "a", "i", "j", "k" };
  args.double_list_flag( ROTATE_FLAG, "Specify a rotation as an angle in degrees counter-clockwise about a vector", &rotate_arg );
  args.limit_list_flag( ROTATE_FLAG, 4, ROTATE_VALS );
  
  const char* SCALE_VALS[] = { "s", "sx", "sy", "sz" };
  args.double_list_flag( SCALE_FLAG, "Specify factor(s) by which to scale mesh about origin", &scale_arg );
  args.limit_list_flag( SCALE_FLAG, 1, SCALE_VALS );
  args.limit_list_flag( SCALE_FLAG, 3, SCALE_VALS + 1 );
  
  const char* TRANSLATE_VALS[] = { "dx", "dy", "dz" };
  args.double_list_flag( TRANSLATE_FLAG, "Specify translation of vertex coordinates.", &translate_arg );
  args.limit_list_flag( TRANSLATE_FLAG, 3, TRANSLATE_VALS );

  args.toggle_flag( 'f', "Do not move fixed vertices.", &freeonly );
  args.toggle_flag( 'k', "Mark boundary vertices as fixed", &skin );

  args.add_required_arg( "input_file" );
  args.add_required_arg( "output_file" );
  
  
  std::vector<std::string> files;
  if (!args.parse_options( argc, argv, files, std::cerr )) {
    args.print_usage( std::cerr );
    exit(1);
  }
  std::string input_file = files[0];
  std::string output_file = files[1];

  MeshImpl mesh;
  MsqError err;
  mesh.read_vtk( input_file.c_str(), err );
  if (err) {
    std::cerr << err << std::endl 
              << "Failed to read file: " << input_file << std::endl;
    return 1;
  }
  
  if (skin.value()) {
    mesh.mark_skin_fixed( err, false );
    if (err) {
      std::cerr << err << std::endl
                << "Failed to skin mesh from file: " << input_file << std::endl;
      return 1;
    }
  }
  
  xform.skip_fixed_vertices( freeonly.value() );
  MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, 0);
  xform.loop_over_mesh( &mesh_and_domain, 0, err );
  if (err) {
    std::cerr << err << std::endl ;
    return 2;
  }
  
  mesh.write_vtk( output_file.c_str(), err );
  if (err) {
    std::cerr << err << std::endl 
                    << "Failed to write file: " << output_file << std::endl;
    return 1;
  }
  
  return 0;
}

void BCDTest::compare_bcd( ObjectiveFunction* OF, string name, const char* mesh_file )
{
  MsqPrintError err(cout);
  size_t i;
  vector<MsqVertex> initial_coords, global_coords, bcd_coords;
  vector<Mesh::VertexHandle> vertex_list;
  
    // set up a smoother
  TerminationCriterion iterations, vertex_movement;
  iterations.add_iteration_limit( 2 );
  vertex_movement.add_absolute_vertex_movement( 1e-3 );

  SolverType global_solver( OF );
  SolverType bcd_solver( OF );
  global_solver.use_global_patch();
  bcd_solver.use_element_on_vertex_patch();
  bcd_solver.do_block_coordinate_descent_optimization();
  global_solver.set_inner_termination_criterion( &vertex_movement );
  bcd_solver.set_inner_termination_criterion( &iterations );
  bcd_solver.set_outer_termination_criterion( &vertex_movement );

  QualityAssessor qa;
  qa.add_quality_assessment( &mMetric );

  InstructionQueue global_q, bcd_q;
  global_q.add_quality_assessor( &qa, err );
  global_q.set_master_quality_improver( &global_solver, err );
  global_q.add_quality_assessor( &qa, err );
  bcd_q.set_master_quality_improver( &bcd_solver, err );
  bcd_q.add_quality_assessor( &qa, err );
  
    // read mesh
  MeshImpl mesh;
  mesh.read_vtk( mesh_file, err ); ASSERT_NO_ERROR(err);
  mesh.get_all_vertices( vertex_list, err ); ASSERT_NO_ERROR(err);
  CPPUNIT_ASSERT(!vertex_list.empty());
  initial_coords.resize( vertex_list.size() );
  mesh.vertices_get_coordinates( arrptr(vertex_list), arrptr(initial_coords), vertex_list.size(), err );
  ASSERT_NO_ERROR(err);
  
    // run global smoother
  global_q.run_instructions( &mesh, err ); 
  ASSERT_NO_ERROR(err);
  mesh.write_vtk( (name + "-gbl.vtk").c_str(), err );
  global_coords.resize( vertex_list.size() );
  mesh.vertices_get_coordinates( arrptr(vertex_list), arrptr(global_coords), vertex_list.size(), err );
  ASSERT_NO_ERROR(err);
  
    // restore initial vertex positions
  for (i = 0; i < vertex_list.size(); ++i) {
    mesh.vertex_set_coordinates( vertex_list[i], initial_coords[i], err );
    ASSERT_NO_ERROR(err);
  }
  
    // run local smoother
  bcd_q.run_instructions( &mesh, err );
  ASSERT_NO_ERROR(err);
  mesh.write_vtk( (name + "-bcd.vtk").c_str(), err );
  bcd_coords.resize( vertex_list.size() );
  mesh.vertices_get_coordinates( arrptr(vertex_list), arrptr(bcd_coords), vertex_list.size(), err );
  ASSERT_NO_ERROR(err);
  
    // compare results
  for (i = 0; i < bcd_coords.size(); ++i)
    CPPUNIT_ASSERT_VECTORS_EQUAL( global_coords[i], bcd_coords[i], 1e-2 );
}

int main( int argc, char* argv[] )
{
  const double default_fraction = 0.05;
  const double zero = 0.0;
  int one = 1;
  CLArgs::ToggleArg allow_invalid( false );
  CLArgs::DoubleRangeArg rand_percent( default_fraction, &zero, 0 );
  CLArgs::IntRangeArg unoptimize( 0, &one, 0 );
  
  CLArgs args( "vtkrandom",
               "Randomize mesh vertex locations.",
               "Read VTK file, randomize locations of containded vertices, and re-write file." );
  args.toggle_flag( INVALID_FLAG, "Allow inverted elements in output", &allow_invalid );
  args.double_flag( PERCENT_FLAG, "fract", "Randomize fraction", &rand_percent );
  args.int_flag( UNOPTIMIZE_FLAG, "N", "Use UnOptimizer with N passes rather than Randomize", &unoptimize );
  add_domain_args( args );
  args.add_required_arg( "input_file" );
  args.add_required_arg( "output_file" );

  std::vector<std::string> files;
  if (!args.parse_options( argc, argv, files, std::cerr )) {
    args.print_usage( std::cerr );
    exit(1);
  }
  std::string input_file = files[0];
  std::string output_file = files[1];
  
  MsqError err;
  MeshImpl mesh;
  mesh.read_vtk( input_file.c_str(), err );
  if (err) {
    std::cerr << "ERROR READING FILE: " << input_file << std::endl
                    << err << std::endl;
    return 2;
  }
  MeshDomain* domain = process_domain_args( &mesh );

  TerminationCriterion tc;
  QualityAssessor qa( false );
  InstructionQueue q;
  Randomize op( rand_percent.value() );
  IdealWeightInverseMeanRatio metric;
  PMeanPTemplate of( 1, &metric );
  UnOptimizer op2( &of );
  if (unoptimize.seen()) {
    tc.add_iteration_limit( unoptimize.value() );
    op2.set_outer_termination_criterion( &tc );
    q.add_preconditioner( &op, err );
    q.set_master_quality_improver( &op2, err );
  }
  else {
    q.set_master_quality_improver( &op, err );
  }
  q.add_quality_assessor( &qa, err );
  q.run_instructions( &mesh, domain, err );
  if (err) {
    std::cerr << err << std::endl;
    return 3;
  }

  int inverted, junk;
  if (qa.get_inverted_element_count( inverted, junk, err ) && inverted ) {
    if (allow_invalid.value())
      std::cerr << "Warning: output mesh contains " << inverted << " inverted elements" << std::endl;
    else {
      std::cerr << "Error: output mesh contains " << inverted << " inverted elements" << std::endl;
      return 4;
    }
  }
  
  mesh.write_vtk( output_file.c_str(), err );
  if (err) {
    std::cerr << "ERROR WRITING FILE: " << output_file << std::endl
                    << err << std::endl;
    return 2;
  }
  
  return 0;
}