From 257a79275c177f80d7a6c8aa1e21cc2fb30f94c9 Mon Sep 17 00:00:00 2001
From: Martin Genet <martin.genet@polytechnique.edu>
Date: Wed, 19 Jun 2019 12:24:19 +0200
Subject: [PATCH] mesh2ugrid now works for scalar & vector functions
---
fedic.py | 4 +-
mesh2ugrid.py | 100 ++++++++++++++++++++++++++++++--------------------
2 files changed, 62 insertions(+), 42 deletions(-)
diff --git a/fedic.py b/fedic.py
index 77a9f8f..dd77ee6 100644
--- a/fedic.py
+++ b/fedic.py
@@ -274,7 +274,7 @@ def fedic(
if (print_refined_mesh):
U.set_allow_extrapolation(True)
U_for_plot = dolfin.interpolate(U, function_space_for_plot)
- U_for_plot.rename("displacement", "a Function")
+ U_for_plot.rename("displacement", "displacement")
file_pvd = dolfin.File(pvd_basename+"-refined__.pvd")
file_pvd << (U_for_plot, float(images_ref_frame))
os.remove(
@@ -932,7 +932,7 @@ def fedic(
if (print_refined_mesh):
U_for_plot = dolfin.interpolate(U, function_space_for_plot)
- U_for_plot.rename("displacement", "a Function")
+ U_for_plot.rename("displacement", "displacement")
file_pvd = dolfin.File(pvd_basename+"-refined__.pvd")
file_pvd << (U_for_plot, float(k_frame))
os.remove(
diff --git a/mesh2ugrid.py b/mesh2ugrid.py
index 4b9a96e..b011352 100644
--- a/mesh2ugrid.py
+++ b/mesh2ugrid.py
@@ -25,11 +25,9 @@ import dolfin_dic as ddic
################################################################################
def mesh2ugrid(
- mesh,
- function_space):
+ mesh):
n_dim = mesh.geometry().dim()
- # n_dim = mesh.ufl_domain().geometric_dimension()
assert (n_dim in (2,3))
# print "n_dim = "+str(n_dim)
@@ -38,56 +36,49 @@ def mesh2ugrid(
n_cells = mesh.num_cells()
# print "n_cells = "+str(n_cells)
+ # Create function space
+ fe = dolfin.FiniteElement(
+ family="CG",
+ cell=mesh.ufl_cell(),
+ degree=1,
+ quad_scheme='default')
+ fs = dolfin.FunctionSpace(
+ mesh,
+ fe)
+
# Store nodes coordinates as numpy array
- n_dofs = function_space.dim()
- np_coordinates = function_space.tabulate_dof_coordinates().reshape([n_dofs, n_dim])
+ n_nodes = fs.dim()
+ assert (n_nodes == n_verts)
+ # print "n_nodes = "+str(n_nodes)
+ np_coordinates = fs.tabulate_dof_coordinates().reshape([n_nodes, n_dim])
# print "np_coordinates = "+str(np_coordinates)
if (n_dim == 2):
- np_coordinates = numpy.hstack((np_coordinates, numpy.zeros([n_dofs, 1])))
+ np_coordinates = numpy.hstack((np_coordinates, numpy.zeros([n_nodes, 1])))
# print "np_coordinates = "+str(np_coordinates)
# Convert nodes coordinates to VTK
- vtk_coordinates = vtk.util.numpy_support.numpy_to_vtk(np_coordinates)
+ vtk_coordinates = vtk.util.numpy_support.numpy_to_vtk(
+ num_array=np_coordinates,
+ deep=1)
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
+ # print "n_points = "+str(vtk_points.GetNumberOfPoints())
# Store connectivity as numpy array
+ n_nodes_per_cell = fs.dofmap().num_element_dofs(0)
+ # print "n_nodes_per_cell = "+str(n_nodes_per_cell)
np_connectivity = numpy.empty(
- [n_cells, n_nodes_per_element],
+ [n_cells, n_nodes_per_cell+1],
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]
+ np_connectivity[i, 0] = n_nodes_per_cell
+ np_connectivity[i,1:] = fs.dofmap().cell_dofs(i)
+ # print "np_connectivity = "+str(np_connectivity)
# 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()
+ np_connectivity = np_connectivity.flatten()
+ # print "np_connectivity = "+str(np_connectivity)
# Convert connectivity to VTK
vtk_connectivity = vtk.util.numpy_support.numpy_to_vtkIdTypeArray(np_connectivity)
@@ -97,6 +88,10 @@ def mesh2ugrid(
vtk_cells.SetCells(n_cells, vtk_connectivity)
# Create unstructured grid and set points and connectivity
+ if (n_dim == 2):
+ vtk_cell_type = vtk.VTK_TRIANGLE
+ elif (n_dim == 3):
+ vtk_cell_type = vtk.VTK_TETRA
ugrid = vtk.vtkUnstructuredGrid()
ugrid.SetPoints(vtk_points)
ugrid.SetCells(vtk_cell_type, vtk_cells)
@@ -109,10 +104,35 @@ def add_function_to_ugrid(
function,
ugrid):
+ # print ugrid.GetPoints()
+
# Convert function values and add as scalar data
+ n_dofs = function.function_space().dim()
+ # print "n_dofs = "+str(n_dofs)
+ n_dim = function.value_size()
+ # print "n_dim = "+str(n_dim)
+ assert (n_dofs/n_dim == ugrid.GetNumberOfPoints()),\
+ "Only CG1 functions can be converted to VTK. Aborting."
np_array = function.vector().get_local()
- vtk_array = vtk.util.numpy_support.numpy_to_vtk(np_array)
+ # print "np_array = "+str(np_array)
+ np_array = np_array.reshape([n_dofs/n_dim, n_dim])
+ # print "np_array = "+str(np_array)
+ vtk_array = vtk.util.numpy_support.numpy_to_vtk(
+ num_array=np_array,
+ deep=1)
vtk_array.SetName(function.name())
- # There are multiple ways of adding this
- ugrid.GetPointData().SetScalars(vtk_array)
+ # print ugrid.GetPoints()
+ ugrid.GetPointData().AddArray(vtk_array)
+ # print ugrid.GetPoints()
+
+################################################################################
+
+def add_functions_to_ugrid(
+ functions,
+ ugrid):
+
+ for function in functions:
+ ddic.add_function_to_ugrid(
+ function=function,
+ ugrid=ugrid)
--
GitLab