New mesh2ugrid function (thanks to Miguel A. Rodriguez)

__init__.py

@@ -9,6 +9,7 @@ from generated_image_expressions_cpp import *
 from plot_strains_vs_radius import *
 from Energy_GeneratedImage import *
 from compute_unwarped_images import *
+from mesh2ugrid import *
 from Problem_ImageRegistration import *
 from compute_quadrature_degree import *
 from Energy_Regularization import *

compute_quadrature_degree.py

@@ -20,8 +20,7 @@ import dolfin_dic as ddic
 
 def compute_quadrature_degree_from_points_count(
         image_filename,
-        mesh_filebasename,
-        mesh_ext="vtk",
+        mesh,
         compute_n_quad=False,
         deg_min=1,
         deg_max=20,
@@ -30,15 +29,13 @@ def compute_quadrature_degree_from_points_count(
     image = myvtk.readImage(
         filename=image_filename,
         verbose=verbose-1)
-    n_points = image.GetNumberOfPoints()
+    image_n_points = image.GetNumberOfPoints()
     image_dimension = myvtk.getImageDimensionality(
         image=image,
         verbose=verbose-1)
 
-    mesh = myvtk.readUGrid(
-        filename=mesh_filebasename+"."+mesh_ext,
-        verbose=verbose-1)
-    n_cells = mesh.GetNumberOfCells()
+    ugrid = ddic.mesh2ugrid(mesh)
+    n_cells = ugrid.GetNumberOfCells()
 
     (cell_locator,
      closest_point,
@@ -46,12 +43,12 @@ def compute_quadrature_degree_from_points_count(
      k_cell,
      subId,
      dist) = myvtk.getCellLocator(
-        mesh=mesh,
+        mesh=ugrid,
         verbose=verbose-1)
 
     point = numpy.empty(3)
     n_pixels_per_cell = numpy.zeros(n_cells)
-    for k_point in xrange(n_points):
+    for k_point in xrange(image_n_points):
         image.GetPoint(k_point, point)
 
         k_cell = cell_locator.FindCell(point)
@@ -66,7 +63,6 @@ def compute_quadrature_degree_from_points_count(
     #if (verbose): print "n_pixels_per_cell_avg = "+str(n_pixels_per_cell_avg)
 
     if (compute_n_quad):
-        mesh = dolfin.Mesh(mesh_filebasename+"."+"xml")
         #n_quads = []
         for degree in xrange(deg_min,deg_max+1):
             if (verbose): print "degree = "+str(degree)
@@ -96,8 +92,7 @@ def compute_quadrature_degree_from_points_count(
 
 def compute_quadrature_degree_from_integral(
         image_filename,
-        mesh=None,
-        mesh_filebasename=None,
+        mesh,
         deg_min=1,
         deg_max=10,
         tol=1e-2,
@@ -110,10 +105,10 @@ def compute_quadrature_degree_from_integral(
     image_dimension = myvtk.getImageDimensionality(
         image=image,
         verbose=verbose-1)
+    assert (image_dimension in (2,3)),\
+        "image_dimension must be 2 or 3. Aborting."
     if (verbose): print "image_dimension = " + str(image_dimension)
 
-    if (mesh is None):
-        mesh = dolfin.Mesh(mesh_filebasename+"."+"xml")
     dX = dolfin.dx(mesh)
 
     first_time = True
@@ -125,7 +120,7 @@ def compute_quadrature_degree_from_integral(
             cell=mesh.ufl_cell(),
             degree=degree,
             quad_scheme="default")
-        if (image_dimension == 2):
+        if   (image_dimension == 2):
             I0 = ddic.ExprIm2(
                 filename=image_filename,
                 element=fe)
@@ -133,8 +128,6 @@ def compute_quadrature_degree_from_integral(
             I0 = ddic.ExprIm3(
                 filename=image_filename,
                 element=fe)
-        else:
-            assert (0), "image_dimension must be 2 or 3. Aborting."
         if not (first_time):
             I0_norm_old = I0_norm
         I0_norm = dolfin.assemble(I0**2 * dX, form_compiler_parameters={'quadrature_degree':degree})**(1./2)

fedic.py

@@ -145,7 +145,7 @@ def fedic(
         if (images_quadrature_from == "points_count"):
             images_quadrature = ddic.compute_quadrature_degree_from_points_count(
                 image_filename=ref_image_filename,
-                mesh_filebasename=mesh_filebasename,
+                mesh=mesh,
                 verbose=1)
         elif (images_quadrature_from == "integral"):
             images_quadrature = ddic.compute_quadrature_degree_from_integral(

fedic2.py

@@ -96,12 +96,12 @@ def fedic2(
         if (images_quadrature_from == "points_count"):
             images_quadrature = ddic.compute_quadrature_degree_from_points_count(
                 image_filename=image_series.get_image_filename(images_ref_frame),
-                mesh_filebasename=mesh_folder+"/"+mesh_basename,
+                mesh=mesh,
                 verbose=1)
         elif (method == "integral"):
             images_quadrature = ddic.compute_quadrature_degree_from_integral(
                 image_filename=self.get_image_filename(images_ref_frame),
-                mesh_filebasename=mesh_folder+"/"+mesh_basename,
+                mesh=mesh,
                 verbose=1)
         else:
             assert (0), "\"images_quadrature_from\" (="+str(images_quadrature_from)+") must be \"points_count\" or \"integral\". Aborting."

mesh2ugrid.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2019                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+### This function inspired by Miguel A. Rodriguez                            ###
+###  https://fenicsproject.org/qa/12933/                                     ###
+###            making-vtk-python-object-from-solution-object-the-same-script ###
+###                                                                          ###
+################################################################################
+
+import dolfin
+import glob
+import numpy
+import vtk
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+def mesh2ugrid(
+        mesh,
+        function_space):
+
+    n_dim = mesh.geometry().dim()
+    # n_dim = mesh.ufl_domain().geometric_dimension()
+    assert (n_dim in (2,3))
+    # print "n_dim = "+str(n_dim)
+
+    n_verts = mesh.num_vertices()
+    # print "n_verts = "+str(n_verts)
+    n_cells = mesh.num_cells()
+    # print "n_cells = "+str(n_cells)
+
+    # Store nodes coordinates as numpy array
+    n_dofs = function_space.dim()
+    np_coordinates = function_space.tabulate_dof_coordinates().reshape([n_dofs, n_dim])
+    # print "np_coordinates = "+str(np_coordinates)
+
+    if (n_dim == 2):
+        np_coordinates = numpy.hstack((np_coordinates, numpy.zeros([n_dofs, 1])))
+        # print "np_coordinates = "+str(np_coordinates)
+
+    # Convert nodes coordinates to VTK
+    vtk_coordinates = vtk.util.numpy_support.numpy_to_vtk(np_coordinates)
+    vtk_points = vtk.vtkPoints()
+    vtk_points.SetData(vtk_coordinates)
+
+    # Check element type
+    dofmap = function_space.dofmap()
+    # print "dofmap = "+str(dofmap)
+    n_nodes_per_element = dofmap.num_element_dofs(0)
+    # print "n_nodes_per_element = "+str(n_nodes_per_element)
+    if   (n_dim == 2):
+        assert (n_nodes_per_element in (3,  6))
+        if   (n_nodes_per_element == 3):
+            vtk_cell_type = vtk.VTK_TRIANGLE
+        elif (n_nodes_per_element == 6):
+            vtk_cell_type = vtk.VTK_QUADRATIC_TRIANGLE
+    elif (n_dim == 3):
+        assert (n_nodes_per_element in (4, 10))
+        if   (n_nodes_per_element == 4):
+            vtk_cell_type = vtk.VTK_TETRA
+        elif (n_nodes_per_element == 10):
+            vtk_cell_type = vtk.VTK_QUADRATIC_TETRA
+
+    # Store connectivity as numpy array
+    np_connectivity = numpy.empty(
+        [n_cells, n_nodes_per_element],
+        dtype=numpy.int)
+    for i in range(n_cells):
+        np_connectivity[i,:] = dofmap.cell_dofs(i)
+
+    # Permute connectivity
+    # (Because node numbering scheme in FEniCS is different from VTK for 2nd order.)
+    # (Explains why the connectivity is first filled as an array and then flattened.)
+    if   (vtk_cell_type == vtk.VTK_QUADRATIC_TRIANGLE):
+        assert (0), "ToCheck. Aborting."
+    elif (vtk_cell_type == vtk.VTK_QUADRATIC_TETRA):
+        PERM_DOLFIN_TO_VTK = [0, 1, 2, 3, 9, 6, 8, 7, 5, 4]
+        np_connectivity = np_connectivity[:, PERM_DOLFIN_TO_VTK]
+
+    # Add left column specifying number of nodes per cell and flatten array
+    np_connectivity = numpy.hstack((numpy.ones([n_cells, 1], dtype=numpy.int)*n_nodes_per_element, np_connectivity)).flatten()
+
+    # Convert connectivity to VTK
+    vtk_connectivity = vtk.util.numpy_support.numpy_to_vtkIdTypeArray(np_connectivity)
+
+    # Create cell array
+    vtk_cells = vtk.vtkCellArray()
+    vtk_cells.SetCells(n_cells, vtk_connectivity)
+
+    # Create unstructured grid and set points and connectivity
+    ugrid = vtk.vtkUnstructuredGrid()
+    ugrid.SetPoints(vtk_points)
+    ugrid.SetCells(vtk_cell_type, vtk_cells)
+
+    return ugrid
+
+################################################################################
+
+def add_function_to_ugrid(
+        function,
+        ugrid):
+
+    # Convert function values and add as scalar data
+    np_array = function.vector().get_local()
+    vtk_array = vtk.util.numpy_support.numpy_to_vtk(np_array)
+    vtk_array.SetName(function.name())
+
+    # There are multiple ways of adding this
+    ugrid.GetPointData().SetScalars(vtk_array)