C++ mesquite::MeshImpl类代码示例

   /* Automatically called by CppUnit before each test function. */
  void setUp()
  {
      // Read a VTK file -- 1 triangle flanked by 1 quad on each side (1 tri + 3 quads)
    mMesh = new Mesquite::MeshImpl;
    mMesh->read_vtk(MESH_FILES_DIR "2D/VTK/hybrid_3quad_1tri.vtk", mErr);
    CPPUNIT_ASSERT(!mErr);

      // Get mesh data
    mMesh->get_all_elements( mElements, mErr );
    CPPUNIT_ASSERT(!mErr);
    mMesh->elements_get_attached_vertices( &mElements[0],
                                           mElements.size(),
                                           mConnectivity,
                                           mOffsets,
                                           mErr );
    CPPUNIT_ASSERT(!mErr);
    
      // Construct list of vertices w/out duplicates from
      // connectivity list.
    std::vector<Mesquite::Mesh::VertexHandle>::iterator new_end;
    mVertices = mConnectivity;
    std::sort( mVertices.begin(), mVertices.end() );
    new_end = std::unique( mVertices.begin(), mVertices.end() );
    mVertices.resize( new_end - mVertices.begin() );
  }

  void test_plane_tri_xz()
     {
       MsqPrintError err(cout); 
       Mesquite::MeshImpl mesh;
       mesh.read_vtk(MESH_FILES_DIR "2D/VTK/tri_5_xz.vtk", err);
       CPPUNIT_ASSERT(!err);

         //create geometry: plane y=5, normal (0,1,0)
       Vector3D pnt(0,-5,0);
       Vector3D s_norm(0,-1,0);
       Mesquite::PlanarDomain msq_geom(s_norm, pnt);
       
         // creates an intruction queue
       InstructionQueue queue1;
       
         //creates a asm quality metric ...
       ConditionNumberQualityMetric smooth;
       
         // ... and builds an objective function with it (untangle)
       LPtoPTemplate smooth_func(&smooth,1,err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         // creates the cg optimization procedures
       ConjugateGradient pass1( &smooth_func, err );
         //pass1->set_patch_type(PatchData::ELEMENTS_ON_VERTEX_PATCH, err,1 ,1);
       pass1.use_global_patch();
       pass1.set_debugging_level(1);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
       QualityAssessor qa=QualityAssessor( &smooth );
       
         //**********Set stopping criterion  5 iterates ****************
       TerminationCriterion sc5;
       sc5.add_iteration_limit( 5 );
       pass1.set_inner_termination_criterion(&sc5);
         //StoppingCriterion sc5(StoppingCriterion::NUMBER_OF_PASSES,5);
         //pass1->set_stopping_criterion(&sc5);
       TerminationCriterion sc_inner;
       sc_inner.add_iteration_limit( 5 );
       pass1.set_inner_termination_criterion(&sc_inner);
         //pass1->set_maximum_iteration(5);
       
       queue1.set_master_quality_improver(&pass1, err); CPPUNIT_ASSERT(!err);
         //********************UNTANGLE*******************************
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         // launches optimization on mesh_set1
       double orig_qa_val=qa.loop_over_mesh(&mesh, &msq_geom, 0, err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
       queue1.run_instructions(&mesh, &msq_geom, err); CPPUNIT_ASSERT(!err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
       double fin_qa_val=qa.loop_over_mesh(&mesh, &msq_geom, 0, err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         //make sure 'quality' improved
       CPPUNIT_ASSERT( (fin_qa_val-orig_qa_val) <= 0.0 );
       print_timing_diagnostics(cout);
     }

int main()
{
  Mesquite::MeshImpl mesh;
  MsqPrintError err(cout);
  mesh.read_vtk(MESH_FILES_DIR "2D/vtk/quads/untangled/square_quad_2.vtk", err);
  if (err) return 1;
  
     //create geometry: plane z=0, normal (0,0,1)
  Vector3D pnt(0,0,5);
  Vector3D s_norm(0,0,1);
  Mesquite::PlanarDomain msq_geom(s_norm, pnt);
  
    // creates an intruction queue
  LaplaceWrapper laplacian_smoother;
  
  mesh.write_vtk("original_mesh.vtk", err); 
  if (err) return 1;
  
    // launches optimization on mesh_set1
  MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &msq_geom);
  laplacian_smoother.run_instructions(&mesh_and_domain, err); 
  if (err) return 1;
 
  mesh.write_vtk("smoothed_mesh.vtk", err); 
  if (err) return 1;
}

  void test_cg_mesh_cond_sphere()
  {
    Mesquite::MeshImpl mesh;
    Mesquite::MsqPrintError err(cout);
    mesh.read_vtk(MESH_FILES_DIR "2D/vtk/quads/untangled/quads_on_sphere_529.vtk", err);
    CPPUNIT_ASSERT(!err);
    
      //create geometry: sphere, center (2,2,0), radius 3
    Vector3D center(2,2,0);
    SphericalDomain msq_geom(center, 3.0);
    
      // creates an intruction queue
    InstructionQueue queue1;
    
      // creates a mean ratio quality metric ...
    ConditionNumberQualityMetric shape;
    UntangleBetaQualityMetric untan;
    
      // ... and builds an objective function with it
    LPtoPTemplate obj_func(&shape, 2, err);
      //Make sure no errors
    CPPUNIT_ASSERT(!err);
      // creates the steepest descent optimization procedures
    ConjugateGradient pass1( &obj_func, err );
      //SteepestDescent* pass2 = new SteepestDescent( obj_func );
    pass1.use_global_patch();
      //Make sure no errors
    CPPUNIT_ASSERT(!err);
    QualityAssessor qa=QualityAssessor( &shape );
   
      //**********Set stopping criterion  5 iterates ****************
      //StoppingCriterion sc5(StoppingCriterion::NUMBER_OF_PASSES,5);
      //pass1->set_stopping_criterion(&sc5);
    TerminationCriterion sc5;
    sc5.add_iteration_limit( 5 );
    pass1.set_inner_termination_criterion(&sc5);
      //CG's debugging print, increase integer to get more print info
    pass1.set_debugging_level(0);
 
      //  queue1.add_preconditioner(pass2, err); CPPUNIT_ASSERT(!err);
    queue1.set_master_quality_improver(&pass1, err); CPPUNIT_ASSERT(!err);
      //Make sure no errors
    CPPUNIT_ASSERT(!err);
      // launches optimization on mesh_set1
    MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &msq_geom);
    double orig_qa_val=qa.loop_over_mesh(&mesh_and_domain, 0, err);
      //Make sure no errors
    CPPUNIT_ASSERT(!err);
    queue1.run_instructions(&mesh_and_domain, err); CPPUNIT_ASSERT(!err);
      //Make sure no errors
    CPPUNIT_ASSERT(!err);
    double fin_qa_val=qa.loop_over_mesh(&mesh_and_domain, 0, err);
      //Make sure no errors
    CPPUNIT_ASSERT(!err);
      //make sure 'quality' improved
    CPPUNIT_ASSERT( fin_qa_val <= orig_qa_val );
  }

  void test_lapl_geo_sphere()
     {
       Mesquite::MeshImpl mesh;
       Mesquite::MsqPrintError err(cout);
       
       mesh.read_vtk(MESH_FILES_DIR "2D/vtk/tris/untangled/Mesquite_geo_10242.vtk", err);
       
         //create geometry sphere:  ratius 1, centered at (0,0,0)
       Vector3D center(0,0,0);
       Mesquite::SphericalDomain msq_geom(center, 1.0);
  
         // creates an intruction queue
       InstructionQueue queue1;

         // creates an edge length metric ...
       EdgeLengthQualityMetric edg_len;
      
         //create the laplacian smoother
       LaplacianSmoother lapl;
         //Make sure no errors
       CPPUNIT_ASSERT(!err);

         //create a quality assessor
       QualityAssessor qa=QualityAssessor( &edg_len );

         //*******Set stopping criterion 10 iterates  ***********
         //StoppingCriterion sc10(StoppingCriterion::NUMBER_OF_PASSES,10);
         //lapl->set_stopping_criterion(&sc10);
       TerminationCriterion sc10;
       sc10.add_iteration_limit( 10 );
       lapl.set_outer_termination_criterion(&sc10);
         //qa, qi, qa
       queue1.set_master_quality_improver(&lapl, err); 
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         // launches optimization on mesh_set1
       MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &msq_geom);
       double orig_qa_val=qa.loop_over_mesh(&mesh_and_domain, 0, err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
       queue1.run_instructions(&mesh_and_domain, err); CPPUNIT_ASSERT(!err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
       double fin_qa_val=qa.loop_over_mesh(&mesh_and_domain, 0, err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         //make sure 'quality' improved
       CPPUNIT_ASSERT( fin_qa_val <= orig_qa_val );
     }

   void test_smart_lapl_sphere()
     {
       Mesquite::MeshImpl mesh;
       Mesquite::MsqPrintError err(cout); 
       mesh.read_vtk(MESH_FILES_DIR "2D/vtk/quads/untangled/quads_on_sphere_529.vtk", err);
       
         //create geometry sphere:  ratius 1, centered at (0,0,0)
       Vector3D center(2,2,0);
       SphericalDomain msq_geom(center, 3.0);
  
         // creates an intruction queue
       InstructionQueue queue1;

         // creates an edge length metric ...
       IdealWeightInverseMeanRatio shape_metric(err);
       LInfTemplate shape_func(&shape_metric);
       
         //create the smart laplacian smoother
       SmartLaplacianSmoother s_lapl(&shape_func);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);

         //*******Set stopping criterion 5 iterates  ***********
       TerminationCriterion sc5;
       sc5.add_iteration_limit( 5 );
       s_lapl.set_outer_termination_criterion(&sc5);
         //qa, qi, qa
       queue1.set_master_quality_improver(&s_lapl, err); CPPUNIT_ASSERT(!err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         // launches optimization on mesh_set1
       QualityAssessor qa=QualityAssessor( &shape_metric );
       MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &msq_geom);
       double orig_val=qa.loop_over_mesh(&mesh_and_domain, 0, err);
       
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
       queue1.run_instructions(&mesh_and_domain, err); CPPUNIT_ASSERT(!err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);

       double final_val= qa.loop_over_mesh(&mesh_and_domain, 0, err);
  
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         //make sure 'quality' improved
       CPPUNIT_ASSERT( final_val < orig_val );
     }

int main()
{
  MsqPrintError err(std::cout);
  Mesquite::MeshImpl mesh;
  mesh.read_vtk(MESH_FILES_DIR "3D/vtk/hexes/untangled/hexes_4by2by2.vtk", err);
  
    // creates an intruction queue
  InstructionQueue queue1;
  
    // creates a mean ratio quality metric ...
  IdealWeightInverseMeanRatio mean_ratio(err);
  if (err) return 1;
  ConditionNumberQualityMetric cond_num;
  mean_ratio.set_averaging_method(QualityMetric::LINEAR);
  
    // ... and builds an objective function with it
    //LInfTemplate* obj_func = new LInfTemplate(mean_ratio);
  LPtoPTemplate obj_func(&mean_ratio, 2, err);
  if (err) return 1;
   // creates the steepest descent optimization procedures
  SteepestDescent pass1( &obj_func );
  pass1.use_global_patch();
  //if (err) return 1;
  //pass1.set_maximum_iteration(6);
  
  QualityAssessor stop_qa=QualityAssessor(&mean_ratio);
  if (err) return 1;
  stop_qa.add_quality_assessment(&cond_num);
  if (err) return 1;
  
  
   //**************Set stopping criterion****************
// StoppingCriterion sc1(&stop_qa,1.0,1.8);
    //StoppingCriterion sc2(StoppingCriterion::NUMBER_OF_PASSES,1);
  TerminationCriterion tc2;
  tc2.add_iteration_limit( 1 );
// CompositeAndStoppingCriterion sc(&sc1,&sc2);
  pass1.set_inner_termination_criterion(&tc2);

  // adds 1 pass of pass1 to mesh_set1
//  queue1.add_preconditioner(pass1, err); 
//  if (err) return 1;
  queue1.add_quality_assessor(&stop_qa,err);
  queue1.set_master_quality_improver(&pass1, err); 
  if (err) return 1;
  queue1.add_quality_assessor(&stop_qa,err);
  if (err) return 1;
  // adds 1 passes of pass2 to mesh_set1
//  mesh_set1.add_quality_pass(pass2);

  mesh.write_vtk("original_mesh.vtk", err); 
  if (err) return 1;
  
    // launches optimization on mesh_set1
  queue1.run_instructions(&mesh, err);
  if (err) return 1;
  
  mesh.write_vtk("smoothed_mesh.vtk", err); 
  if (err) return 1;
  
  return 0;
}

int main()
{     
    /* Reads a Mesh file */
  const char *file_name = 
//      MESH_FILES_DIR "2D/vtk/tris/untangled/equil_tri2.vtk";
//      MESH_FILES_DIR "2D/vtk/tris/untangled/tri_20258.vtk";
//      MESH_FILES_DIR "3D/vtk/tets/untangled/tet_1.vtk";
//      MESH_FILES_DIR "3D/vtk/hexes/untangled/cube_tet_2.vtk";
     MESH_FILES_DIR "3D/vtk/tets/untangled//tire.vtk";
  printf("Loading mesh set 1\n");
  MsqPrintError err( cout );
  Mesquite::MeshImpl mesh;
  mesh.read_vtk(file_name, err);
  if (err) return 1;
  
    // Creates an intruction queue
    //  printf("Creating instruction queue\n");
  InstructionQueue queue1;

  // Creates a condition number quality metric 
  //  printf("Creating quality metric\n");
  ConditionNumberQualityMetric cond_no;
  
  // Create the NonSmooth Steepest Descent procedures
  //  printf("creating optimizer\n");
  NonSmoothDescent minmax_method( &cond_no );

  // Set a termination criterion
  TerminationCriterion tc2;
  tc2.add_iteration_limit( 1 );
  minmax_method.set_outer_termination_criterion(&tc2);
  // Set up the quality assessor
  //  printf("Setting up the quality assessor\n");
  QualityAssessor quality_assessor=QualityAssessor(&cond_no);

  // assess the quality of the initial mesh
  queue1.add_quality_assessor(&quality_assessor, err); 
  if (err) return 1;

  // Set the max min method to be the master quality improver
  queue1.set_master_quality_improver(&minmax_method, err); 
  if (err) return 1;

  // assess the quality of the final mesh
  queue1.add_quality_assessor(&quality_assessor, err); 
  if (err) return 1;

  // write out the original mesh
  //  printf("Writing out the original mesh\n");
  mesh.write_vtk("original_mesh.vtk", err); 
  if (err) return 1;

  // launches optimization on mesh_set1
  //  printf("Running the instruction queue\n");
  queue1.run_instructions(&mesh, err); 
  if (err) return 1;

  // write out the smoothed mesh
  //  printf("Writing out the final mesh\n");
  mesh.write_vtk("smoothed_mesh.vtk", err); 
  if (err) return 1;
  
  return 0;
}

bool smooth_mixed_mesh( const char* filename )
{
  Mesquite::MsqPrintError err(cout);
  
  // print a little output so we know when we died
  std::cout << 
  "**************************************************************************" 
  << std::endl << 
  "* Smoothing: " << filename
  << std::endl  <<
  "**************************************************************************" 
  << std::endl;
  
  // The instruction queue to set up
  InstructionQueue Q;
  
  // Use numeric approx of derivitives until analytic solutions
  // are working for pyramids
  IdealWeightInverseMeanRatio mr_metric(err);
  //sRI_DFT dft_metric;
  UntangleBetaQualityMetric un_metric(0);
  CPPUNIT_ASSERT(!err);
  
    // Create Mesh object
  Mesquite::MeshImpl mesh;
  mesh.read_vtk(filename, err);
  CPPUNIT_ASSERT(!err);

  // Set up a preconditioner
  LInfTemplate pre_obj_func( &un_metric );
  ConjugateGradient precond( &pre_obj_func, err ); CPPUNIT_ASSERT(!err);
  precond.use_element_on_vertex_patch();
  TerminationCriterion pre_term, pre_outer;
  //pre_term.add_relative_quality_improvement( 0.1 );
  pre_term .add_iteration_limit( 3 );
  pre_outer.add_iteration_limit( 1 );
  CPPUNIT_ASSERT(!err);
  precond.set_inner_termination_criterion( &pre_term );
  precond.set_outer_termination_criterion( &pre_outer );
  //precond.use_element_on_vertex_patch();

  // Set up objective function
  LPtoPTemplate obj_func(&mr_metric, 1, err);
  CPPUNIT_ASSERT(!err);

  // Create solver
  FeasibleNewton solver( &obj_func, true );
  CPPUNIT_ASSERT(!err);
  solver.use_global_patch();
  CPPUNIT_ASSERT(!err);

  // Set stoping criteria for solver
  TerminationCriterion tc_inner;
  tc_inner.add_relative_quality_improvement( 0.25 );
  solver.set_inner_termination_criterion(&tc_inner);
   
  TerminationCriterion tc_outer;
  tc_outer.add_iteration_limit( 1 );
  CPPUNIT_ASSERT(!err);
  solver.set_outer_termination_criterion(&tc_outer);

  // Create a QualityAssessor
  Mesquite::QualityAssessor qa;
  qa.add_quality_assessment( &mr_metric );
  qa.add_quality_assessment( &un_metric );
  Q.add_quality_assessor( &qa, err ); 
  CPPUNIT_ASSERT(!err);
 
  // Add untangler to queue
  Q.add_preconditioner( &precond, err ); CPPUNIT_ASSERT(!err);
  Q.add_quality_assessor( &qa, err ); 
  CPPUNIT_ASSERT(!err);
 
  // Add solver to queue
  Q.set_master_quality_improver(&solver, err); 
  CPPUNIT_ASSERT(!err);
  Q.add_quality_assessor( &qa, err ); 
  CPPUNIT_ASSERT(!err);
 
  // And smooth...
  Q.run_instructions(&mesh, err); 
  CPPUNIT_ASSERT(!err);

  return false;
}

int main( int argc, char* argv[] )
{
  unsigned i;
  const char* input_file = MESH_FILES_DIR "3D/VTK/mixed-hex-pyr-tet.vtk";
  if (argc == 2)
    input_file = argv[1];
  else if (argc != 1)
  {
    std::cerr << "Invalid arguments.\n";
    return 2;
  }
  
  
  Mesquite::MsqPrintError err(cout);
  IdealWeightMeanRatio m1;
  IdealWeightInverseMeanRatio m2(err);
  ConditionNumberQualityMetric m3;
  QualityMetric* metrics[] = { &m1, &m2, &m3, 0 };

    // Read Mesh
  Mesquite::MeshImpl mesh;
  mesh.read_vtk(MESH_FILES_DIR "3D/VTK/12-pyramid-unit-sphere.vtk", err);
  CPPUNIT_ASSERT(!err);
  Mesquite::MeshImpl ideal_mesh;
  ideal_mesh.read_vtk(MESH_FILES_DIR "3D/VTK/12-pyramid-unit-sphere.vtk", err);
  CPPUNIT_ASSERT(!err);

    // Check that the mesh read correctly, and contains what is
    // expected later.

    // Get mesh data
    // Expecting file to contain 12 pyramid elements constructed
    // from 15 vertices.
  std::vector<Mesh::VertexHandle> vert_array;
  std::vector<Mesh::ElementHandle> elem_array;
  std::vector<size_t> conn_offsets;
  mesh.get_all_elements( elem_array, err ); 
  CPPUNIT_ASSERT(!err);
  CPPUNIT_ASSERT( elem_array.size() == 12 );
  mesh.elements_get_attached_vertices( &elem_array[0],
                                        elem_array.size(),
                                        vert_array,
                                        conn_offsets,
                                        err );
  CPPUNIT_ASSERT(!err);
  CPPUNIT_ASSERT(vert_array.size() == 60);
  CPPUNIT_ASSERT(conn_offsets.size() == 13);
  EntityTopology type_array[12];
  mesh.elements_get_topologies( &elem_array[0], type_array, 12, err );
  CPPUNIT_ASSERT(!err);
  
    // Verify element types and number of vertices
  for (i = 0; i < 12; ++i)
  {
    CPPUNIT_ASSERT( type_array[i] == PYRAMID );
    CPPUNIT_ASSERT( conn_offsets[i] == 5*i );
  }
  
    // All pyramids should share a common apex, at the
    // center of the sphere
  Mesh::VertexHandle apex_handle = vert_array[4];
  for (i = 1; i < 12; ++i)
  {
    CPPUNIT_ASSERT( vert_array[5*i+4] == apex_handle );
  }
  
    // Verify that apex is at origin and all other vertices are
    // on unit sphere
  MsqVertex vertices[60];
  mesh.vertices_get_coordinates( &vert_array[0], vertices, 60, err );
  CPPUNIT_ASSERT(!err);
  for (i = 0; i < 60; ++i)
  {
    if (vert_array[i] == apex_handle)
      CPPUNIT_ASSERT( vertices[i].within_tolerance_box( Vector3D(0,0,0), 1e-6 ) );
    else
      CPPUNIT_ASSERT( fabs(1.0 - vertices[i].length()) < 1e-6 );
  }
  
    // Try smoothing w/out moving the free vertex and verify that
    // the smoother didn't move the vertex
  Vector3D position(0,0,0);
  for (i = 0; metrics[i] != NULL; ++i)
    CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, apex_handle, position, metrics[i] ) );
  
    // Now try moving the vertex and see if the smoother moves it back
    // to the origin
  position.set( 0.1, 0.1, 0.1 );
  for (i = 0; metrics[i] != NULL; ++i)
    CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, apex_handle, position, metrics[i] ) );
  
    // Now try moving the vertex further and see if the smoother moves it back
    // to the origin
  position.set( 0.3, 0.3, 0.3 );
  for (i = 0; metrics[i] != NULL; ++i)
    CPPUNIT_ASSERT( !smooth_mesh( &mesh, &ideal_mesh, apex_handle, position, metrics[i] ) );

    // Now try smoothing a real mixed mesh
  CPPUNIT_ASSERT( !smooth_mixed_mesh( input_file ) );

//.........这里部分代码省略.........

int main( )
{
  Mesquite::MeshImpl mesh;
  MsqPrintError err(cout);
  mesh.read_vtk(VTK_2D_DIR "quads/tangled/tangled_quad.vtk", err);
  if (err) return 1;
  
  // Set Domain Constraint
  Vector3D pnt(0,0,5);
  Vector3D s_norm(0,0,1);
  PlanarDomain msq_geom(s_norm, pnt);
                                                                              
    // creates an intruction queue
  InstructionQueue queue1;
  
    // creates a mean ratio quality metric ...
  ConditionNumberQualityMetric shape_metric;
  UntangleBetaQualityMetric untangle(2);
  Randomize pass0(.05);
    // ... and builds an objective function with it
    //LInfTemplate* obj_func = new LInfTemplate(shape_metric);
  LInfTemplate obj_func(&untangle);
  LPtoPTemplate obj_func2(&shape_metric, 2, err);
  if (err) return 1;
    // creates the steepest descent optimization procedures
  ConjugateGradient pass1( &obj_func, err );
  if (err) return 1;
  
    //SteepestDescent* pass2 = new SteepestDescent( obj_func2 );
  ConjugateGradient pass2( &obj_func2, err );
  if (err) return 1;
  pass2.use_element_on_vertex_patch();
  if (err) return 1;
  pass2.use_global_patch();
  if (err) return 1;
  QualityAssessor stop_qa=QualityAssessor(&shape_metric);
  QualityAssessor stop_qa2=QualityAssessor(&shape_metric);
  if (brief_output) {
    stop_qa.disable_printing_results();
    stop_qa2.disable_printing_results();
  }
  
  stop_qa.add_quality_assessment(&untangle);
    // **************Set stopping criterion**************
    //untangle beta should be 0 when untangled
  TerminationCriterion sc1;
  sc1.add_relative_quality_improvement( 0.000001 );
  TerminationCriterion sc3;
  sc3.add_iteration_limit( 10 );
  TerminationCriterion sc_rand;
  sc_rand.add_iteration_limit( 1 );
  
    //StoppingCriterion sc1(&stop_qa,-1.0,.0000001);
    //StoppingCriterion sc3(&stop_qa2,.9,1.00000001);
    //StoppingCriterion sc2(StoppingCriterion::NUMBER_OF_PASSES,10);
    //StoppingCriterion sc_rand(StoppingCriterion::NUMBER_OF_PASSES,1);
    //either until untangled or 10 iterations
  pass0.set_outer_termination_criterion(&sc_rand);
  pass1.set_outer_termination_criterion(&sc1);
  pass2.set_inner_termination_criterion(&sc3);
  
    // adds 1 pass of pass1 to mesh_set1
  queue1.add_quality_assessor(&stop_qa,err); 
  if (err) return 1;
    //queue1.add_preconditioner(pass0,err);MSQ_CHKERR(err);
    //queue1.add_preconditioner(pass1,err);MSQ_CHKERR(err);
    //queue1.set_master_quality_improver(pass2, err); MSQ_CHKERR(err);
  queue1.set_master_quality_improver(&pass1, err);
  if (err) return 1;
  queue1.add_quality_assessor(&stop_qa2,err);
  if (err) return 1;
  if (write_output)
    mesh.write_vtk("original_mesh.vtk", err);
  if (err) return 1;
  
    // launches optimization on mesh_set1
  MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &msq_geom);
  queue1.run_instructions(&mesh_and_domain, err);
  if (err) return 1;
  
  if (write_output)
    mesh.write_vtk("smoothed_mesh.vtk", err); 
  if (err) return 1;
  
  if (!brief_output)
    print_timing_diagnostics(cout);
  return 0;
}

int main(int argc, char* argv[])
{
  const char* input_file = DEFAULT_INPUT;
  const char* output_file = NULL;
  switch (argc) {
    default:
      help(argv[0]);
    case 3:
      if (!strcmp(argv[2],"-h"))
        help(argv[0]);
      output_file = argv[2];
    case 2:
      if (!strcmp(argv[1],"-h"))
        help(argv[0]);
      input_file = argv[1];
    case 1:
      ;
  }

    /* Read a VTK Mesh file */
  MsqPrintError err(cout);
  Mesquite::MeshImpl mesh;
  mesh.read_vtk( input_file, err);
  if (err) {
    std::cerr << input_file << ": file read failed" << endl;
    return 1;
  }
  
  double min_ini, avg_ini, max_ini;
  tet_dihedral_angle_ratios( mesh, min_ini, avg_ini, max_ini, err );
  
    /* Run optimizer */
  ViscousCFDTetShapeWrapper smoother(1e-2);
  smoother.run_instructions( &mesh, err);
  if (err) return 1;
  
  double min_fin, avg_fin, max_fin;
  tet_dihedral_angle_ratios( mesh, min_fin, avg_fin, max_fin, err );
  
  cout << "\tMinimum \tAverage \tMaximum "<<endl
       << "Init\t"<<min_ini<<'\t'<<avg_ini<<'\t'<<max_ini<<endl
       << "Fini\t"<<min_fin<<'\t'<<avg_fin<<'\t'<<max_fin<<endl;
  
  if (output_file) {
    mesh.write_vtk( output_file, err );
    if (err) {
      std::cerr << output_file << ": file write failed" << endl;
      return 1;
    }
    else {
      std::cout << "Wrote file: " << output_file << endl;
    }
  }
      
  if (min_ini > min_fin || avg_ini > avg_fin || max_ini > max_fin) {
    std::cerr << "Dihedral handle ratio decreased" << endl;
    return 1;
  }
  
  return 0;
}

  void test_plane_quad_tangled()
     {
       Mesquite::MeshImpl mesh;
       MsqPrintError err(cout); 
       mesh.read_vtk(MESH_FILES_DIR "2D/VTK/tangled_quad.vtk", err);
       CPPUNIT_ASSERT(!err);

         //create geometry: plane z=5, normal (0,0,1)
       Vector3D pnt(0,0,5);
       Vector3D s_norm(0,0,1);
       Mesquite::PlanarDomain msq_geom(s_norm, pnt);
       
         // creates an intruction queue
       InstructionQueue queue1, queue2;

         //creates a mean ratio quality metric ...
       ConditionNumberQualityMetric shape;
       UntangleBetaQualityMetric untan(.1);
  
         // ... and builds an objective function with it (untangle)
       LInfTemplate untan_func(&untan);
       LPtoPTemplate shape_func(&shape,2,err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         // creates the cg optimization procedures
       ConjugateGradient pass1( &untan_func, err );
       ConjugateGradient pass2( &shape_func, err );
       pass1.use_element_on_vertex_patch();
       pass2.use_global_patch();
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
       QualityAssessor stop_qa=QualityAssessor( &untan );
       QualityAssessor qa=QualityAssessor( &shape );
         //turn off printing if print flag not set.
       if(pF==0){
         stop_qa.disable_printing_results();
         qa.disable_printing_results();
       }
       //**********Set stopping criterion  untangle ver small ********
       //StoppingCriterion sc_qa(&stop_qa,-100,MSQ_MIN);
       //pass1->set_stopping_criterion(&sc_qa);
       TerminationCriterion sc_of;
       sc_of.add_iteration_limit( 10 );
       pass1.set_outer_termination_criterion(&sc_of);
       
         //**********Set stopping criterion  5 iterates ****************
         //StoppingCriterion sc5(StoppingCriterion::NUMBER_OF_PASSES,5);
         //pass2->set_stopping_criterion(&sc5);
       TerminationCriterion sc5;
       sc5.add_iteration_limit( 5 );
       pass2.set_inner_termination_criterion(&sc5);
         //pass2->set_maximum_iteration(5);
         //TerminationCriterion sc_inner;
         //sc_inner.add_iteration_limit( 5 );
         //pass2->set_inner_termination_criterion(&sc_inner);
       queue1.set_master_quality_improver(&pass1, err); CPPUNIT_ASSERT(!err);
       queue2.set_master_quality_improver(&pass2, err); CPPUNIT_ASSERT(!err);
         //********************UNTANGLE*******************************
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         // launches optimization on mesh_set1
       double orig_qa_val=stop_qa.loop_over_mesh(&mesh, &msq_geom, 0, err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
       queue1.run_instructions(&mesh, &msq_geom, err); CPPUNIT_ASSERT(!err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
       double fin_qa_val=stop_qa.loop_over_mesh(&mesh, &msq_geom, 0, err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         //make sure 'quality' improved
       CPPUNIT_ASSERT( (fin_qa_val-orig_qa_val) <= 0.0 );
         //make sure sc_qa really was satisfied
       CPPUNIT_ASSERT( fin_qa_val <= MSQ_MIN );
       
         //********************SMOOTH*******************************
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         // launches optimization on mesh_set1
       orig_qa_val=qa.loop_over_mesh(&mesh, &msq_geom, 0, err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
       queue2.run_instructions(&mesh, &msq_geom, err); CPPUNIT_ASSERT(!err);
       //Make sure no errors
       CPPUNIT_ASSERT(!err);
       fin_qa_val=qa.loop_over_mesh(&mesh, &msq_geom, 0, err);
         //Make sure no errors
       CPPUNIT_ASSERT(!err);
         //make sure 'quality' improved
       CPPUNIT_ASSERT( (fin_qa_val-orig_qa_val) <= 0.0 );
       print_timing_diagnostics(cout);
     }

  void test_vertices()
  {
    size_t nbVert = mVertices.size();
    CPPUNIT_ASSERT_EQUAL(9,(int)nbVert);

    Mesquite::MsqVertex correct_coords[9], coords[9];
    correct_coords[0].set(1,0,0);
    correct_coords[1].set(0,1.732,0);
    correct_coords[2].set(-1,0,0);
    correct_coords[3].set(-1,-2,0);
    correct_coords[4].set(1,-2,0);
    correct_coords[5].set(2.732,1,0);
    correct_coords[6].set(1.732,2.732,0);
    correct_coords[7].set(-1.732,2.732,0);
    correct_coords[8].set(-2.732,1,0);

    mMesh->vertices_get_coordinates(&mVertices[0], coords, nbVert, mErr);
    CPPUNIT_ASSERT(!mErr);
    for (size_t i=0; i<nbVert; ++i) {
      for (int j=0; j<3; ++j)
        CPPUNIT_ASSERT_DOUBLES_EQUAL(coords[i][j], correct_coords[i][j], .01);
    }

    coords[3].set(2.,3.,4.);
    mMesh->vertex_set_coordinates(mVertices[3], coords[3], mErr);
    CPPUNIT_ASSERT(!mErr);
    Mesquite::MsqVertex coords_2;
    mMesh->vertices_get_coordinates(&mVertices[3], &coords_2, 1, mErr);
    CPPUNIT_ASSERT(!mErr);
    for (int j=0; j<3; ++j)
      CPPUNIT_ASSERT_DOUBLES_EQUAL(coords[3][j], coords_2[j], 1e-6);
  }

  void test_vertex_byte()
  {
    size_t nbVert = mVertices.size();
    size_t i;
	unsigned char* bytes = new unsigned char[nbVert];
    mMesh->vertices_get_byte(&mVertices[0], bytes, nbVert, mErr); 
    CPPUNIT_ASSERT(!mErr);

    // Asserts all vertex bytes are initialised to 0. 
    for (i=0; i<nbVert; ++i)
      CPPUNIT_ASSERT(bytes[i] == 0);

    // Test various vertex byte read / write routines.
    bytes[3] |= 4;
    mMesh->vertices_set_byte(&mVertices[0], bytes, nbVert, mErr); 
    CPPUNIT_ASSERT(!mErr);
    mMesh->vertex_set_byte(mVertices[5], 8, mErr); 
    CPPUNIT_ASSERT(!mErr);
    unsigned char byte;
    mMesh->vertex_get_byte(mVertices[3], &byte, mErr);
    CPPUNIT_ASSERT(!mErr);
    CPPUNIT_ASSERT(byte == 4);
    mMesh->vertices_get_byte(&mVertices[0], bytes, nbVert, mErr);
    CPPUNIT_ASSERT(!mErr);
    for (i=0; i<nbVert; ++i) {
      if (i==3)
        CPPUNIT_ASSERT(bytes[i] == 4);
      else if (i==5)
        CPPUNIT_ASSERT(bytes[i] == 8);
      else
        CPPUNIT_ASSERT(bytes[i] == 0);
    }

    delete [] bytes;
  }