/*
 * Copyright 2003 Sandia Corporation.
 * Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
 * license for use of this work by or on behalf of the
 * U.S. Government. Redistribution and use in source and binary forms, with
 * or without modification, are permitted provided that this Notice and any
 * statement of authorship are reproduced on all copies.
 */

#include "vtkMeshQuality.h"

#include "vtkDebugLeaks.h"
#include "vtkFieldData.h"
#include "vtkPoints.h"
#include "vtkSmartPointer.h"
#include "vtkTestErrorObserver.h"
#include "vtkTestUtilities.h"
#include "vtkTriangle.h"
#include "vtkUnstructuredGrid.h"
#include "vtkUnstructuredGridReader.h"

static int DumpQualityStats(vtkMeshQuality* iq, const char* arrayname)
{
  cout << "  cardinality: "
       << iq->GetOutput()->GetFieldData()->GetArray(arrayname)->GetComponent(0, 4)
       << "  , range: " << iq->GetOutput()->GetFieldData()->GetArray(arrayname)->GetComponent(0, 0)
       << "  -  " << iq->GetOutput()->GetFieldData()->GetArray(arrayname)->GetComponent(0, 2)
       << endl;

  cout << "  average: " << iq->GetOutput()->GetFieldData()->GetArray(arrayname)->GetComponent(0, 1)
       << "  , standard deviation: "
       << sqrt(fabs(iq->GetOutput()->GetFieldData()->GetArray(arrayname)->GetComponent(0, 3)))
       << endl;

  return 0;
}

int MeshQuality(int argc, char* argv[])
{
  vtkUnstructuredGridReader* mr = vtkUnstructuredGridReader::New();
  vtkUnstructuredGrid* ug;
  vtkMeshQuality* iq = vtkMeshQuality::New();
  char* fname = vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/uGridEx.vtk");

  mr->SetFileName(fname);
  mr->Update();

  ug = mr->GetOutput();
  iq->SetInputConnection(mr->GetOutputPort());
  iq->SaveCellQualityOn();
  cout << "SaveCellQuality: " << iq->GetSaveCellQuality() << endl;

  if (ug->GetNumberOfCells())
  {
    cout << endl;
    cout << "Triangle quality of mesh" << endl;
    cout << mr->GetFileName() << endl;

    iq->SetTriangleQualityMeasureToArea();
    iq->Update();
    cout << " Area:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToEdgeRatio();
    iq->Update();
    cout << " Edge Ratio:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToAspectRatio();
    iq->Update();
    cout << " Aspect Ratio:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToRadiusRatio();
    iq->Update();
    cout << " Radius Ratio:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToAspectFrobenius();
    iq->Update();
    cout << " Frobenius Norm:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToMinAngle();
    iq->Update();
    cout << " Minimal Angle:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToMaxAngle();
    iq->Update();
    cout << " Maximal Angle:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToCondition();
    iq->Update();
    cout << " Condition:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToScaledJacobian();
    iq->Update();
    cout << " ScaledJacobian:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToRelativeSizeSquared();
    iq->Update();
    cout << " RelativeSizeSquared:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToShape();
    iq->Update();
    cout << " Shape:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToShapeAndSize();
    iq->Update();
    cout << " ShapeAndSize:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    iq->SetTriangleQualityMeasureToDistortion();
    iq->Update();
    cout << " Distortion:" << endl;
    DumpQualityStats(iq, "Mesh Triangle Quality");

    cout << endl;
    cout << "Quadrilatedral quality of mesh" << endl;
    cout << mr->GetFileName() << endl;

    iq->SetQuadQualityMeasureToEdgeRatio();
    iq->Update();
    cout << " Edge Ratio:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToAspectRatio();
    iq->Update();
    cout << " Aspect Ratio:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToRadiusRatio();
    iq->Update();
    cout << " Radius Ratio:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToMedAspectFrobenius();
    iq->Update();
    cout << " Average Frobenius Norm:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToMaxAspectFrobenius();
    iq->Update();
    cout << " Maximal Frobenius Norm:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToMaxEdgeRatios();
    iq->Update();
    cout << " Max Edge Ratios:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToSkew();
    iq->Update();
    cout << " Skew:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToTaper();
    iq->Update();
    cout << " Taper:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToWarpage();
    iq->Update();
    cout << " Warpage:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToArea();
    iq->Update();
    cout << " Area:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToStretch();
    iq->Update();
    cout << " Stretch:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToMinAngle();
    iq->Update();
    cout << " Min Angle:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToMaxAngle();
    iq->Update();
    cout << " Max Angle:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToOddy();
    iq->Update();
    cout << " Oddy:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToCondition();
    iq->Update();
    cout << " Condition:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToJacobian();
    iq->Update();
    cout << " Jacobian:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToScaledJacobian();
    iq->Update();
    cout << " Scaled Jacobian:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToShear();
    iq->Update();
    cout << " Shear:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToShape();
    iq->Update();
    cout << " Shape:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToRelativeSizeSquared();
    iq->Update();
    cout << " Relative Size Squared:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToShapeAndSize();
    iq->Update();
    cout << " Shape And Size:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToShearAndSize();
    iq->Update();
    cout << " Shear And Size:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    iq->SetQuadQualityMeasureToDistortion();
    iq->Update();
    cout << " Distortion:" << endl;
    DumpQualityStats(iq, "Mesh Quadrilateral Quality");

    cout << endl;
    cout << "Tetrahedral quality of mesh" << endl;
    cout << mr->GetFileName() << endl;

    iq->SetTetQualityMeasureToEdgeRatio();
    iq->Update();
    cout << " Edge Ratio:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");

    iq->SetTetQualityMeasureToAspectRatio();
    iq->Update();
    cout << " Aspect Ratio:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");

    iq->SetTetQualityMeasureToRadiusRatio();
    iq->Update();
    cout << " Radius Ratio:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");

    iq->SetTetQualityMeasureToAspectFrobenius();
    iq->Update();
    cout << " Frobenius Norm:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");

    iq->SetTetQualityMeasureToMinAngle();
    iq->Update();
    cout << " Minimal Dihedral Angle:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");

    iq->SetTetQualityMeasureToCollapseRatio();
    iq->Update();
    cout << " Collapse Ratio:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");
    cout << endl;

    iq->SetTetQualityMeasureToAspectBeta();
    iq->Update();
    cout << " Aspect Beta:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");
    cout << endl;

    iq->SetTetQualityMeasureToAspectGamma();
    iq->Update();
    cout << " Aspect Gamma:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");
    cout << endl;

    iq->SetTetQualityMeasureToVolume();
    iq->Update();
    cout << " Volume:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");
    cout << endl;

    iq->SetTetQualityMeasureToCondition();
    iq->Update();
    cout << " Condition:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");
    cout << endl;

    iq->SetTetQualityMeasureToJacobian();
    iq->Update();
    cout << " Jacobian:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");
    cout << endl;

    iq->SetTetQualityMeasureToScaledJacobian();
    iq->Update();
    cout << " Scaled Jacobian:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");
    cout << endl;

    iq->SetTetQualityMeasureToShape();
    iq->Update();
    cout << " Shape:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");
    cout << endl;

    iq->SetTetQualityMeasureToRelativeSizeSquared();
    iq->Update();
    cout << " Relative Size Squared:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");
    cout << endl;

    iq->SetTetQualityMeasureToShapeAndSize();
    iq->Update();
    cout << " Shape And Size:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");
    cout << endl;

    iq->SetTetQualityMeasureToDistortion();
    iq->Update();
    cout << " Distortion:" << endl;
    DumpQualityStats(iq, "Mesh Tetrahedron Quality");
    cout << endl;

    cout << "Hexahedral quality of mesh" << endl;
    cout << mr->GetFileName() << endl;

    iq->SetHexQualityMeasureToEdgeRatio();
    iq->Update();
    cout << " Edge Ratio:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToMedAspectFrobenius();
    iq->Update();
    cout << " Med Aspect Frobenius:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToMaxAspectFrobenius();
    iq->Update();
    cout << " Max Aspect Frobenius:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToMaxEdgeRatios();
    iq->Update();
    cout << " Max Edge Ratios:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToSkew();
    iq->Update();
    cout << " Skew:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToTaper();
    iq->Update();
    cout << " Taper:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToVolume();
    iq->Update();
    cout << " Volume:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToStretch();
    iq->Update();
    cout << " Stretch:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToDiagonal();
    iq->Update();
    cout << " Diagonal:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToDimension();
    iq->Update();
    cout << " Dimension:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToOddy();
    iq->Update();
    cout << " Oddy:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToCondition();
    iq->Update();
    cout << " Condition:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToJacobian();
    iq->Update();
    cout << " Jacobian:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToScaledJacobian();
    iq->Update();
    cout << " Scaled Jacobian:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToShear();
    iq->Update();
    cout << " Shear:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToShape();
    iq->Update();
    cout << " Shape:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToRelativeSizeSquared();
    iq->Update();
    cout << " Relative Size Squared:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToShapeAndSize();
    iq->Update();
    cout << " Shape And Size:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToShearAndSize();
    iq->Update();
    cout << " Shear And Size:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;

    iq->SetHexQualityMeasureToDistortion();
    iq->Update();
    cout << " Distortion:" << endl;
    DumpQualityStats(iq, "Mesh Hexahedron Quality");
    cout << endl;
  }

  // Exercise remaining methods for coverage
  iq->Print(cout);

  // Check for warnings
  vtkSmartPointer<vtkTest::ErrorObserver> warningObserver =
    vtkSmartPointer<vtkTest::ErrorObserver>::New();

  iq->AddObserver(vtkCommand::WarningEvent, warningObserver);
  iq->SetTriangleQualityMeasure(100000);
  iq->Update();

  if (warningObserver->GetWarning())
  {
    std::cout << "Caught expected warning: " << warningObserver->GetWarningMessage();
  }
  else
  {
    std::cout << "Failed to catch expected warning" << std::endl;
    return EXIT_FAILURE;
  }
  warningObserver->Clear();

  iq->SetQuadQualityMeasure(100000);
  iq->Update();

  if (warningObserver->GetWarning())
  {
    std::cout << "Caught expected warning: " << warningObserver->GetWarningMessage();
  }
  else
  {
    std::cout << "Failed to catch expected warning" << std::endl;
    return EXIT_FAILURE;
  }
  warningObserver->Clear();

  iq->SetTetQualityMeasure(100000);
  iq->Update();

  if (warningObserver->GetWarning())
  {
    std::cout << "Caught expected warning: " << warningObserver->GetWarningMessage();
  }
  else
  {
    std::cout << "Failed to catch expected warning" << std::endl;
    return EXIT_FAILURE;
  }
  warningObserver->Clear();

  iq->SetHexQualityMeasure(100000);
  iq->Update();

  if (warningObserver->GetWarning())
  {
    std::cout << "Caught expected warning: " << warningObserver->GetWarningMessage();
  }
  else
  {
    std::cout << "Failed to catch expected warning" << std::endl;
    return EXIT_FAILURE;
  }
  warningObserver->Clear();

  iq->Delete();
  mr->Delete();
  delete[] fname;

  return 0;
}
