Add mesh quality functions from vtk as fluidity diagnostics.

diagnostics/Makefile.dependencies

@@ -60,7 +60,7 @@ Mesh_Diagnostics.o ../include/mesh_diagnostics.mod: Mesh_Diagnostics.F90 \
    ../include/diagnostic_source_fields.mod ../include/fdebug.h \
    ../include/field_options.mod ../include/fields.mod ../include/fldebug.mod \
    ../include/global_parameters.mod ../include/halos_numbering.mod \
-   ../include/state_module.mod
+   ../include/mesh_quality.mod ../include/state_module.mod
 
 ../include/metric_diagnostics.mod: Metric_Diagnostics.o
 	@true

diagnostics/Mesh_Diagnostics.F90

@@ -35,13 +35,15 @@ module mesh_diagnostics
   use fields
   use state_module
   use field_options
+  use mesh_quality
   use diagnostic_source_fields
 
   implicit none
 
   private
 
   public :: calculate_column_ids, calculate_universal_column_ids
+  public :: calculate_mesh_quality
 
 contains
 
@@ -96,4 +98,52 @@ subroutine calculate_universal_column_ids(state, s_field)
 
   end subroutine calculate_universal_column_ids
 
+  subroutine calculate_mesh_quality(state, s_field)
+    type(state_type), intent(in) :: state
+    type(scalar_field), intent(inout) :: s_field
+
+    type(vector_field), pointer :: positions
+
+    integer :: measure
+    character(len=FIELD_NAME_LEN) :: measure_name
+
+    positions => extract_vector_field(state,"Coordinate")
+
+    if (element_count(s_field) /= node_count(s_field)) then
+       FLAbort("Mesh quality measures must be on a P0 mesh")
+    end if
+
+    call get_option(trim(complete_field_path(trim(s_field%option_path))) // &
+                    "/algorithm[0]/quality_function/name", measure_name)
+
+    select case(trim(measure_name))
+    case("radius_ratio")
+       measure=VTK_QUALITY_RADIUS_RATIO
+    case("aspect_ratio")
+       measure=VTK_QUALITY_ASPECT_RATIO
+    case("aspect_frobenius")
+       measure=VTK_QUALITY_ASPECT_FROBENIUS
+    case("edge_ratio")
+       measure=VTK_QUALITY_EDGE_RATIO
+    case("condition")
+       measure=VTK_QUALITY_CONDITION
+    case("min_angle")
+       measure=VTK_QUALITY_MIN_ANGLE
+    case("max_angle")
+       measure=VTK_QUALITY_MAX_ANGLE
+    case("shape")
+       measure=VTK_QUALITY_SHAPE
+    case("min_angl")
+       measure=VTK_QUALITY_SHAPE_AND_SIZE
+    case("area_or_volume")
+       measure=VTK_QUALITY_AREA
+    case default
+        FLAbort("Unknown quality function for "//trim(s_field%name) )
+    end select
+
+    call get_mesh_quality(positions, s_field, measure)
+
+  end subroutine calculate_mesh_quality
+
+
 end module mesh_diagnostics

femtools/Makefile.dependencies

@@ -764,6 +764,13 @@ Mesh_Files.o ../include/mesh_files.mod: Mesh_Files.F90 ../include/elements.mod \
    ../include/state_module.mod ../include/write_gmsh.mod \
    ../include/write_triangle.mod
 
+../include/mesh_quality.mod: Mesh_Quality.o
+	@true
+
+Mesh_Quality.o ../include/mesh_quality.mod: Mesh_Quality.F90 \
+   ../include/element_numbering.mod ../include/fdebug.h ../include/fields.mod \
+   ../include/fldebug.mod
+
 ../include/metric_tools.mod: Metric_tools.o
 	@true
 

femtools/Makefile.in

@@ -106,7 +106,8 @@ OBJS = Dgtools.o Coordinates.o EventCounter.o \
   Fields_Halos.o Profiler.o Profiler_Fortran.o Streamfunction.o \
   GMSH_Common.o Read_GMSH.o Write_GMSH.o \
   Exodusii_C_Interface.o Exodusii_F_Interface.o Exodusii_Common.o Read_Exodusii.o \
-  Mesh_Files.o Vertical_Extrapolation.o
+  Mesh_Files.o Vertical_Extrapolation.o \
+  Mesh_Quality.o Mesh_Quality_C.o
 
 
 # objects to be included in libfemtools:

femtools/Mesh_Quality.F90

@@ -0,0 +1,105 @@
+!    Copyright (C) 2006 Imperial College London and others.
+!    
+!    Please see the AUTHORS file in the main source directory for a full list
+!    of copyright holders.
+!
+!    Prof. C Pain
+!    Applied Modelling and Computation Group
+!    Department of Earth Science and Engineering
+!    Imperial College London
+!
+!    amcgsoftware@imperial.ac.uk
+!    
+!    This library is free software; you can redistribute it and/or
+!    modify it under the terms of the GNU Lesser General Public
+!    License as published by the Free Software Foundation,
+!    version 2.1 of the License.
+!
+!    This library is distributed in the hope that it will be useful,
+!    but WITHOUT ANY WARRANTY; without even the implied warranty of
+!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+!    Lesser General Public License for more details.
+!
+!    You should have received a copy of the GNU Lesser General Public
+!    License along with this library; if not, write to the Free Software
+!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+!    USA
+#include "fdebug.h"
+module mesh_quality
+
+  use iso_c_binding
+  use FLdebug
+  use element_numbering, only: FAMILY_SIMPLEX
+  use fields
+
+
+  interface
+     subroutine mesh_quality_c(dim, n_nodes, n_elements, connectivity_len,&
+       measure, points, connectivity, quality) bind(c)
+       use iso_c_binding
+       integer (c_int) :: dim, n_nodes, n_elements, connectivity_len, measure
+       real (c_double) :: points(dim, n_nodes)
+       integer (c_int) :: connectivity(connectivity_len)
+       real (c_double) :: quality(n_elements)
+     end subroutine mesh_quality_c
+  end interface
+
+  private
+
+  public :: get_mesh_quality
+
+
+  integer, public :: VTK_QUALITY_EDGE_RATIO = 0, &
+       VTK_QUALITY_ASPECT_RATIO = 1, &
+       VTK_QUALITY_RADIUS_RATIO = 2, &
+       VTK_QUALITY_ASPECT_FROBENIUS = 3, &
+       VTK_QUALITY_MED_ASPECT_FROBENIUS = 4, &
+       VTK_QUALITY_MAX_ASPECT_FROBENIUS = 5, &
+       VTK_QUALITY_MIN_ANGLE = 6, &
+       VTK_QUALITY_COLLAPSE_RATIO = 1, &
+       VTK_QUALITY_MAX_ANGLE = 8, &
+       VTK_QUALITY_CONDITION = 9, &
+       VTK_QUALITY_SCALED_JACOBIAN = 10, &
+       VTK_QUALITY_SHEAR = 11, &
+       VTK_QUALITY_RELATIVE_SIZE_SQUARED = 12, &
+       VTK_QUALITY_SHAPE = 13, &
+       VTK_QUALITY_SHAPE_AND_SIZE = 14, &
+       VTK_QUALITY_DISTORTION = 15, &
+       VTK_QUALITY_MAX_EDGE_RATIO = 16, &
+       VTK_QUALITY_SKEW = 17, &
+       VTK_QUALITY_TAPER = 18, &
+       VTK_QUALITY_ASPECT_VOLUME = 19, &
+       VTK_QUALITY_ASPECT_STRETCH = 20, &
+       VTK_QUALITY_ASPECT_DIAGONAL = 21, &
+       VTK_QUALITY_ASPECT_DIMENSION = 22, &
+       VTK_QUALITY_ASPECT_ODDY = 23, &
+       VTK_QUALITY_ASPECT_SHEAR_AND_SIZE = 24, &
+       VTK_QUALITY_ASPECT_JACOBIAN = 25, &
+       VTK_QUALITY_ASPECT_WARPAGE = 26, &
+       VTK_QUALITY_ASPECT_GAMMA = 27, &
+       VTK_QUALITY_AREA = 28, &
+       VTK_QUALITY_ASPECT_BETA = 29
+
+contains
+
+  subroutine get_mesh_quality(positions, s_field, quality_measure)
+    integer, intent(inout) :: quality_measure
+    type(vector_field), intent(in) :: positions 
+    type(scalar_field), intent(inout) :: s_field
+
+    assert(element_count(positions) == element_count(s_field))
+    assert(node_count(s_field) == element_count(s_field))
+
+    if (positions%mesh%shape%numbering%family /= FAMILY_SIMPLEX&
+         .or. positions%mesh%shape%loc /= positions%mesh%shape%dim+1) then
+       FLAbort("Trying to get mesh quality for a mesh which isn't linear simplicial. This isn't currently supported.")
+    endif
+
+    call mesh_quality_c(positions%dim, node_count(positions), ele_count(positions),&
+         size(positions%mesh%ndglno), quality_measure,&
+         positions%val, positions%mesh%ndglno,&
+         s_field%val)
+
+  end subroutine get_mesh_quality
+
+end module mesh_quality

femtools/Mesh_Quality_C.cpp

@@ -0,0 +1,101 @@
+//    Copyright (C) 2006 Imperial College London and others.
+//   
+//    Please see the AUTHORS file in the main source directory for a full list
+//    of copyright holders.
+//
+//    Prof. C Pain
+//    Applied Modelling and Computation Group
+//    Department of Earth Science and Engineering
+//    Imperial College London
+//
+//    amcgsoftware@imperial.ac.uk
+//    
+//    This library is free software; you can redistribute it and/or
+//    modify it under the terms of the GNU Lesser General Public
+//    License as published by the Free Software Foundation,
+//    version 2.1 of the License.
+//
+//    This library is distributed in the hope that it will be useful,
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+//    Lesser General Public License for more details.
+//
+//    You should have received a copy of the GNU Lesser General Public
+//    License along with this library; if not, write to the Free Software
+//    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+//    USA
+
+#include "vtkCell.h"
+#include "vtkDataArray.h"
+#include "vtkPoints.h"
+#include "vtkCellData.h"
+#include "vtkDataSet.h"
+#include "vtkUnstructuredGrid.h"
+#include "vtkMeshQuality.h"
+
+#include <stdio.h>
+
+
+extern "C" {
+void mesh_quality_c(int* dim, int* n_nodes, int* n_elements, int* connectivity_len,
+		    int* measure, double* points, int* connectivity, double* quality) {
+  vtkPoints* pts = vtkPoints::New();
+  vtkUnstructuredGrid* ugrid = vtkUnstructuredGrid::New();
+  pts->SetDataTypeToDouble();
+  pts->Allocate(*n_nodes);
+  double x[3]={0.0,0.0,0.0};
+  for (vtkIdType i=0; i<*n_nodes; ++i) {
+    for (int j=0;j<*dim; j++) {
+      x[j] = points[(*dim)*i+j];
+    }
+    pts->SetPoint(i,x);
+  }
+  ugrid->SetPoints(pts);
+
+  int cell_type;
+  vtkIdType cell[20];
+  if (*dim==2) {
+    cell_type = VTK_TRIANGLE;
+  } else if (*dim==3) {
+    cell_type = VTK_TETRA;
+  }
+
+  ugrid->Allocate(0);
+
+  for(vtkIdType i=0;i<*n_elements;++i) {
+    for (int j=0;j<*dim+1;++j) {
+      // off by one due to fortran indexing
+      cell[j] = connectivity[(*dim+1)*i+j]-1;
+    }
+    ugrid->InsertNextCell(cell_type,*dim+1,cell);
+  }
+
+  vtkMeshQuality* filter = vtkMeshQuality::New();
+
+#if VTK_MAJOR_VERSION <= 5
+  filter->SetInput(ugrid);
+#else
+  filter->SetInputData(ugrid);
+#endif
+
+  filter->SetTriangleQualityMeasure(*measure);
+  if (*measure == VTK_QUALITY_AREA) {
+    filter->SetTetQualityMeasure(VTK_QUALITY_VOLUME);
+  } else {
+    filter->SetTetQualityMeasure(*measure);
+  }
+  filter->Update();
+  vtkDataSet* vgrid=filter->GetOutput();
+  vtkDataArray* data = vgrid->GetCellData()->GetArray(0);
+
+  for (vtkIdType i=0;i<data->GetNumberOfTuples();++i) {
+    quality[i] = data->GetTuple1(i);
+  }
+
+  filter->Delete();
+  ugrid->Delete();
+  pts->Delete();
+
+}
+
+}

schemas/diagnostic_algorithms.rnc

@@ -1382,6 +1382,33 @@ universal_column_id_algorithm =
       }
    )
 
+mesh_quality_algorithm =
+   (
+      ## Mesh Quality function. Basics taken from VTK. See the documentation
+      ## there for further details of specific function choices.
+      ##
+      ## The output mesh for this field must be discontinous P0 and the
+      ## underlying CoordinateMesh must consist of linear triangles
+      ##  or tetrahedra only.
+
+      element algorithm {
+        attribute name { "mesh_quality" },
+        attribute material_phase_support { "single" },
+	(
+	  element quality_function{ attribute name{"radius_ratio" } }|
+	  element quality_function{ attribute name{"aspect_ratio" } }|
+	  element quality_function{ attribute name{"aspect_frobenius" } }|
+	  element quality_function{ attribute name{"edge_ratio" } }|
+	  element quality_function{ attribute name{"condition" } }|
+          element quality_function{ attribute name{"min_angle" } }|		
+          element quality_function{ attribute name{"max_angle" } }|
+          element quality_function{ attribute name{"shape" } }|		
+          element quality_function{ attribute name{"shape_and_size" } }|
+	  element quality_function{ attribute name{"area_or_volume"} }
+	)
+      }
+   )
+
 viscous_dissipation_algorithm =
    (
       ## algorithm for calculating the viscous dissipation.
@@ -1542,6 +1569,7 @@ scalar_diagnostic_algorithms |= projection_scalar_potential_algorithm
 scalar_diagnostic_algorithms |= particle_reynolds_number_algorithm
 scalar_diagnostic_algorithms |= apparent_density_algorithm
 scalar_diagnostic_algorithms |= node_halo_diagnostic_algorithm
+scalar_diagnostic_algorithms |= mesh_quality_algorithm
 scalar_diagnostic_algorithms |= lumped_mass_smoothed_scalar_algorithm
 scalar_diagnostic_algorithms |= l2norm_algorithm
 scalar_diagnostic_algorithms |= helmholtz_smoothed_scalar_algorithm