CG3 triangle mesh support and VTK visualization of CG2 and CG3 meshes…

…/functions (#532)

* Added vtkfile to python layer

* Added high order VTK-write support for triangles, still tons of work to do to make this robust

* More flake8 fixes

* Some tmp fixes in c++ to compile without warnings

* Hopefully last warning fix of the day. Not smart to work while tired

* Debatable CPP demo fix

Loads of similar code with slight variations, since the dolfin::generation::xxxMesh is initialized in a very different way than meshes from points

* Added demo to toctree and made headline

* Added visualization of third order cg trianglemeshes

* Flake8fix

* Added high order VTK-write support for triangles, still tons of work to do to make this robust

* Some tmp fixes in c++ to compile without warnings

* Hopefully last warning fix of the day. Not smart to work while tired

* Add a partitioning to python interface and some tests  (#510)

* test read mesh data

* add partition functions to python layer

* add distributed mesh test

* parametrize subset_comm in test

* use enum for graph partitioners

* some small fix and docstring improvement

* reorder parameters for consistency

* use type hint for return

* Changes using existing infrastructure for cell coordinates

* Fixed remainder of rebase

* Bigger rewrite of VTK writer, and mesh initialization

Moved all vtk permutations to separate function. Changing how the file writers (also xdmf and hdf5) access the perturbation map. Moved mesh degree determination to separate function

* Some minor changes after looking at diff

* Removed unused variable

* Fixed visualization for quad mesh, and higher order triangle mesh

* Removed unused variable

* More work on CG-3. Single cell works, demo does not. Unsure why

* Revert change caused by order in vtk::mapping being wrong initially

* Changed test as I mean that this is the wrong orientation

* Added unit test for cg2 and cg3 meshes

* Added some documentation to the demo

* Revert ordering of quads to lexiographic from counter-clockwise. Added note for future PR, and cleanup before PR of this branch

* Changed cross product in compute_quadrilateral

as the old implementation assumed that the permutation was (0,1,2,3) in Mesh.cpp, while the vtk_mapping was (0,1,3,2) (Due to lexiographic ordering in quads). Removing the duplication of mappings caused this change

* Skip parallel testing of higher order meshes

* VTK writer can now write functions for higher order meshes

* Removed demo

* Actually removing file

* Removed demo from docs

* Fix typo

* Some clang-format fixes and other requests by Garth

* More clang-format

* Removed unused import and added missing space

* Removed unusedd c++ variable

* Re-expose VTKFile and add some additional notes on how support differs from XDMF and VTk

* MIssed something in merge conflict

* More remainders of merge

* New attempt on getting clang-format to work

cpp/dolfin/io/HDF5File.cpp

@@ -20,6 +20,7 @@
 #include <dolfin/function/FunctionSpace.h>
 #include <dolfin/la/PETScVector.h>
 #include <dolfin/la/utils.h>
+#include <dolfin/mesh/CoordinateDofs.h>
 #include <dolfin/mesh/DistributedMeshTools.h>
 #include <dolfin/mesh/Mesh.h>
 #include <dolfin/mesh/MeshEntity.h>
@@ -291,7 +292,8 @@ void HDF5File::write(const mesh::Mesh& mesh, int cell_dim,
     const auto& global_vertices = mesh.topology().global_indices(0);
 
     // Permutation to VTK ordering
-    const std::vector<std::int8_t> perm = mesh::vtk_mapping(cell_type);
+    const std::vector<std::uint8_t> perm
+        = mesh.coordinate_dofs().cell_permutation();
 
     if (cell_dim == tdim or !mpi_io)
     {

cpp/dolfin/io/VTKFile.cpp

@@ -159,12 +159,14 @@ std::string init(const mesh::Mesh& mesh, const std::string filename,
                                       counter, filename, ".vtu");
   clear_file(vtu_filename);
 
-  // Number of cells and vertices
+  // Number of cells
   const std::size_t num_cells = mesh.topology().ghost_offset(cell_dim);
-  const std::size_t num_vertices = mesh.topology().ghost_offset(0);
+
+  // Number of points in mesh (can be more than the number of vertices)
+  const int num_nodes = mesh.geometry().points().rows();
 
   // Write headers
-  vtk_header_open(num_vertices, num_cells, vtu_filename);
+  vtk_header_open(num_nodes, num_cells, vtu_filename);
 
   return vtu_filename;
 }
@@ -390,22 +392,24 @@ void write_point_data(const function::Function& u, const mesh::Mesh& mesh,
   std::ostringstream ss;
   ss << std::scientific;
   ss << std::setprecision(16);
-  for (auto& vertex : mesh::MeshRange(mesh, 0))
+  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& points
+      = mesh.geometry().points();
+  for (int i = 0; i < points.rows(); ++i)
   {
-    if (rank == 1 && dim == 2)
+    if (rank == 1 and dim == 2)
     {
       // Append 0.0 to 2D vectors to make them 3D
-      for (std::size_t i = 0; i < 2; i++)
-        ss << values(vertex.index(), i) << " ";
+      for (std::size_t j = 0; j < 2; j++)
+        ss << values(i, j) << " ";
       ss << 0.0 << "  ";
     }
-    else if (rank == 2 && dim == 4)
+    else if (rank == 2 and dim == 4)
     {
       // Pad 2D tensors with 0.0 to make them 3D
-      for (std::size_t i = 0; i < 2; i++)
+      for (std::size_t j = 0; j < 2; j++)
       {
-        ss << values(vertex.index(), (2 * i + 0)) << " ";
-        ss << values(vertex.index(), (2 * i + 1)) << " ";
+        ss << values(i, (2 * j + 0)) << " ";
+        ss << values(i, (2 * j + 1)) << " ";
         ss << 0.0 << " ";
       }
       ss << 0.0 << " ";
@@ -415,8 +419,8 @@ void write_point_data(const function::Function& u, const mesh::Mesh& mesh,
     else
     {
       // Write all components
-      for (std::size_t i = 0; i < dim; i++)
-        ss << values(vertex.index(), i) << " ";
+      for (std::size_t j = 0; j < dim; j++)
+        ss << values(i, j) << " ";
       ss << " ";
     }
   }

cpp/dolfin/io/VTKFile.h

@@ -35,8 +35,8 @@ namespace io
 
 /// Output of meshes and functions in VTK format
 
-/// XML format for visualisation purposes. It is not suitable to
-/// checkpointing as it may decimate some data.
+/// XML format is suitable for visualisation of higher order geometries.
+/// It is not suitable to checkpointing as it may decimate some data.
 
 class VTKFile
 {

cpp/dolfin/io/VTKWriter.cpp

@@ -13,6 +13,8 @@
 #include <dolfin/function/FunctionSpace.h>
 #include <dolfin/la/PETScVector.h>
 #include <dolfin/la/utils.h>
+#include <dolfin/mesh/Connectivity.h>
+#include <dolfin/mesh/CoordinateDofs.h>
 #include <dolfin/mesh/Mesh.h>
 #include <dolfin/mesh/MeshEntity.h>
 #include <dolfin/mesh/MeshFunction.h>
@@ -31,31 +33,47 @@ namespace
 {
 //-----------------------------------------------------------------------------
 // Get VTK cell type
-std::uint8_t vtk_cell_type(const mesh::Mesh& mesh, std::size_t cell_dim)
+std::uint8_t vtk_cell_type(const mesh::Mesh& mesh, std::size_t cell_dim,
+                           std::size_t cell_order)
 {
   // Get cell type
   mesh::CellType cell_type = mesh::cell_entity_type(mesh.cell_type(), cell_dim);
 
   // Determine VTK cell type
-  std::uint8_t vtk_cell_type = 0;
-  if (cell_type == mesh::CellType::tetrahedron)
-    vtk_cell_type = 10;
-  else if (cell_type == mesh::CellType::hexahedron)
-    vtk_cell_type = 12;
-  else if (cell_type == mesh::CellType::quadrilateral)
-    vtk_cell_type = 9;
-  else if (cell_type == mesh::CellType::triangle)
-    vtk_cell_type = 5;
-  else if (cell_type == mesh::CellType::interval)
-    vtk_cell_type = 3;
-  else if (cell_type == mesh::CellType::point)
-    vtk_cell_type = 1;
-  else
+  switch (cell_type)
   {
+  case mesh::CellType::tetrahedron:
+    return 10;
+  case mesh::CellType::hexahedron:
+    return 12;
+  case mesh::CellType::quadrilateral:
+  {
+    switch (cell_order)
+    {
+    case 1:
+      return 9;
+    default:
+      return 70;
+    }
+  }
+  case mesh::CellType::triangle:
+  {
+    switch (cell_order)
+    {
+    case 1:
+      return 5;
+    default:
+      return 69;
+    }
+  }
+  case mesh::CellType::interval:
+    return 3;
+  case mesh::CellType::point:
+    return 1;
+  default:
     throw std::runtime_error("Unknown cell type");
+    return 0;
   }
-
-  return vtk_cell_type;
 }
 //----------------------------------------------------------------------------
 // Write cell data (ascii)
@@ -107,11 +125,11 @@ void write_ascii_mesh(const mesh::Mesh& mesh, std::size_t cell_dim,
                       std::string filename)
 {
   const std::size_t num_cells = mesh.topology().ghost_offset(cell_dim);
-  const std::size_t num_cell_vertices = mesh::num_cell_vertices(
-      mesh::cell_entity_type(mesh.cell_type(), cell_dim));
+  const int element_degree = mesh.degree();
 
   // Get VTK cell type
-  const std::size_t _vtk_cell_type = vtk_cell_type(mesh, cell_dim);
+  const std::size_t _vtk_cell_type
+      = vtk_cell_type(mesh, cell_dim, element_degree);
 
   // Open file
   std::ofstream file(filename.c_str(), std::ios::app);
@@ -126,11 +144,10 @@ void write_ascii_mesh(const mesh::Mesh& mesh, std::size_t cell_dim,
   file << "<DataArray  type=\"Float64\"  NumberOfComponents=\"3\"  format=\""
        << "ascii"
        << "\">";
-  for (auto& v : mesh::MeshRange(mesh, 0))
-  {
-    Eigen::Vector3d p = mesh.geometry().x(v.index());
-    file << p[0] << " " << p[1] << " " << p[2] << "  ";
-  }
+  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor> points
+      = mesh.geometry().points();
+  for (int i = 0; i < points.rows(); ++i)
+    file << points(i, 0) << " " << points(i, 1) << " " << points(i, 2) << "  ";
   file << "</DataArray>" << std::endl << "</Points>" << std::endl;
 
   // Write cell connectivity
@@ -139,13 +156,20 @@ void write_ascii_mesh(const mesh::Mesh& mesh, std::size_t cell_dim,
        << "ascii"
        << "\">";
 
-  mesh::CellType celltype = mesh::cell_entity_type(mesh.cell_type(), cell_dim);
-  const std::vector<std::int8_t> perm = mesh::vtk_mapping(celltype);
-  const int num_vertices = mesh::cell_num_entities(celltype, 0);
-  for (auto& c : mesh::MeshRange(mesh, cell_dim))
+  const mesh::Connectivity& connectivity_g
+      = mesh.coordinate_dofs().entity_points();
+  Eigen::Ref<const Eigen::Array<std::int32_t, Eigen::Dynamic, 1>>
+      cell_connections = connectivity_g.connections();
+  const Eigen::Ref<const Eigen::Array<std::int32_t, Eigen::Dynamic, 1>> pos_g
+      = connectivity_g.entity_positions();
+  const std::vector<std::uint8_t> perm
+      = mesh.coordinate_dofs().cell_permutation();
+  int num_nodes = perm.size();
+
+  for (int j = 0; j < mesh.num_entities(mesh.topology().dim()); ++j)
   {
-    for (int i = 0; i < num_vertices; ++i)
-      file << c.entities(0)[perm[i]] << " ";
+    for (int i = 0; i < num_nodes; ++i)
+      file << cell_connections(pos_g(j) + perm[i]) << " ";
     file << " ";
   }
   file << "</DataArray>" << std::endl;
@@ -155,7 +179,7 @@ void write_ascii_mesh(const mesh::Mesh& mesh, std::size_t cell_dim,
        << "ascii"
        << "\">";
   for (std::size_t offsets = 1; offsets <= num_cells; offsets++)
-    file << offsets * num_cell_vertices << " ";
+    file << offsets * num_nodes << " ";
   file << "</DataArray>" << std::endl;
 
   // Write cell type

cpp/dolfin/io/XDMFFile.h

@@ -61,6 +61,8 @@ class HDF5File;
 /// data. Output of data in parallel is supported.
 ///
 /// XDMF is not suitable for checkpointing as it may decimate some data.
+/// XDMF is not suitable for higher order geometries, as their currently
+/// only supports 1st and 2nd order geometries.
 
 // FIXME: Set read mode when creating file object?
 

cpp/dolfin/io/xdmf_write.cpp

@@ -210,7 +210,8 @@ std::vector<std::int64_t> compute_topology_data(const mesh::Mesh& mesh,
   // Get mesh communicator
   MPI_Comm comm = mesh.mpi_comm();
 
-  const std::vector<std::int8_t> perm = mesh::vtk_mapping(mesh.cell_type());
+  const std::vector<std::uint8_t> perm
+      = mesh.coordinate_dofs().cell_permutation();
   const int tdim = mesh.topology().dim();
   const auto& global_vertices = mesh.topology().global_indices(0);
   if (dolfin::MPI::size(comm) == 1 or cell_dim == tdim)

cpp/dolfin/mesh/Mesh.cpp

@@ -18,6 +18,7 @@
 #include <dolfin/common/MPI.h>
 #include <dolfin/common/Timer.h>
 #include <dolfin/common/utils.h>
+#include <dolfin/mesh/cell_types.h>
 
 using namespace dolfin;
 using namespace dolfin::mesh;
@@ -140,30 +141,12 @@ Mesh::Mesh(MPI_Comm comm, mesh::CellType type,
   }
 
   // Permutation from VTK to DOLFIN order for cell geometric nodes
-  // FIXME: should do this also for quad/hex
-  // FIXME: remove duplication in mesh::vtk_mapping()
-  std::vector<std::uint8_t> cell_permutation = {0, 1, 2, 3, 4, 5, 6, 7};
+  std::vector<std::uint8_t> cell_permutation;
+  cell_permutation = mesh::vtk_mapping(type, cells.cols());
 
-  // Infer if the mesh has P2 geometry (P1 has num_vertices_per_cell ==
-  // cells.cols())
-  if (num_vertices_per_cell != cells.cols())
-  {
-    if (type == mesh::CellType::triangle and cells.cols() == 6)
-    {
-      _degree = 2;
-      cell_permutation = {0, 1, 2, 5, 3, 4};
-    }
-    else if (type == mesh::CellType::tetrahedron and cells.cols() == 10)
-    {
-      _degree = 2;
-      cell_permutation = {0, 1, 2, 3, 9, 6, 8, 7, 5, 4};
-    }
-    else
-    {
-      throw std::runtime_error(
-          "Mismatch between cell type and number of vertices per cell");
-    }
-  }
+  // Find degree of mesh
+  // FIXME: degree should probably be in MeshGeometry
+  _degree = mesh::cell_degree(type, cells.cols());
 
   // Get number of nodes (global)
   const std::uint64_t num_points_global = MPI::sum(comm, points.rows());

cpp/dolfin/mesh/cell_types.cpp

@@ -482,7 +482,7 @@ mesh::cell_entity_closure(mesh::CellType cell_type)
   return entity_closure;
 }
 //-----------------------------------------------------------------------------
-std::vector<std::int8_t> mesh::vtk_mapping(mesh::CellType type)
+std::vector<std::uint8_t> mesh::vtk_mapping(mesh::CellType type, int num_nodes)
 {
   switch (type)
   {
@@ -491,17 +491,86 @@ std::vector<std::int8_t> mesh::vtk_mapping(mesh::CellType type)
   case mesh::CellType::interval:
     return {0, 1};
   case mesh::CellType::triangle:
-    return {0, 1, 2};
+    switch (num_nodes)
+    {
+    case 3:
+      return {0, 1, 2};
+    case 6:
+      return {0, 1, 2, 5, 3, 4};
+    case 10:
+      return {0, 1, 2, 7, 8, 3, 4, 6, 5, 9};
+    default:
+      throw std::runtime_error("Unknown cell type.");
+    }
   case mesh::CellType::tetrahedron:
-    return {0, 1, 2, 3};
+    switch (num_nodes)
+    {
+    case 4:
+      return {0, 1, 2, 3};
+    default:
+      throw std::runtime_error("Higher order tetrahedron not supported");
+    }
   case mesh::CellType::quadrilateral:
-    return {0, 1, 3, 2};
+    switch (num_nodes)
+    {
+    case 4:
+    {
+      // FIXME: Note that this is not counter clockwise ordering (CC),
+      // as performed by vtk (used by gmsh and other mesh generators),
+      // but lexicographic (LG). A convert function from (CC) to (LG) will
+      // be created in a future PR.
+      return {0, 1, 3, 2};
+    }
+    default:
+      throw std::runtime_error("Higher order quadrilateral not supported");
+    }
   case mesh::CellType::hexahedron:
-    return {0, 1, 3, 2, 4, 5, 7, 6};
+    switch (num_nodes)
+    {
+    case 8:
+      return {0, 1, 3, 2, 4, 5, 7, 6};
+    default:
+      throw std::runtime_error("Higher order hexahedron not supported");
+    }
+  default:
+    throw std::runtime_error("Unknown cell type.");
+  }
+}
+//-----------------------------------------------------------------------------
+int mesh::cell_degree(mesh::CellType type, int num_nodes)
+{
+  switch (type)
+  {
+  case mesh::CellType::point:
+    return 1;
+  case mesh::CellType::interval:
+    return 1;
+  case mesh::CellType::triangle:
+    if (num_nodes == 3)
+      return 1;
+    else if (num_nodes == 6)
+      return 2;
+    else if (num_nodes == 10)
+      return 3;
+    else
+      throw std::runtime_error("Unknown cell type.");
+  case mesh::CellType::tetrahedron:
+    if (num_nodes == 4)
+      return 1;
+    else
+      throw std::runtime_error("Higher order tetrahedron not supported");
+  case mesh::CellType::quadrilateral:
+    if (num_nodes == 4)
+      return 1;
+    else
+      throw std::runtime_error("Higher order quadrilateral not supported");
+  case mesh::CellType::hexahedron:
+    if (num_nodes == 8)
+      return 1;
+    else
+      throw std::runtime_error("Higher order hexahedron not supported");
   default:
     throw std::runtime_error("Unknown cell type.");
   }
-
-  return std::vector<std::int8_t>();
 }
 //-----------------------------------------------------------------------------