initialize U with from file is working

ImageIterator.py

@@ -101,82 +101,13 @@ class ImageIterator():
                 #self.printer.print_var("k_frame_old",k_frame_old,-1)
 
                 if (self.initialize_U_from_file):
-                    print k_frame
-
-                    # xdmf_filename = self.initialize_U_from_file+"_"+str(k_frame).zfill(2)+".xmf"
-                    # print xdmf_filename
-                    # xdmf_file = dolfin.XDMFFile(xdmf_filename)
-                    # xdmf_file.read(self.problem.mesh)   #can not read type of cell
-
-                    # xdmf_filename = self.initialize_U_from_file+"_"+str(k_frame).zfill(2)+".h5"
-                    # print xdmf_filename
-                    # xdmf_file = dolfin.HDF5File(dolfin.mpi_comm_world(),xdmf_filename, 'r')
-                    # xdmf_file.read(self.problem.mesh, "grid",True) #unknown cell type
-
-                    # xdmf_filename = self.initialize_U_from_file+".vtu"
-                    # print xdmf_filename
-                    # xdmf_file = dolfin.XDMFFile(xdmf_filename)
-                    # xdmf_file.read(self.problem.mesh) #can not read vtu file
-
-                    # xdmf_filename = self.initialize_U_from_file+"_"+str(k_frame).zfill(2)+".xmf"
-                    # print xdmf_filename
-                    # xdmf_file = dolfin.XDMFFile(xdmf_filename)
-                    # mesh = dolfin.MeshValueCollectionDouble()
-                    # xdmf_file.read(mesh,'U')   #can not read type of cell
-
-                    # mesh_test = dolfin.MeshValueCollection()
-                    # xdmf_file.read(self.problem.U, "U")
-                    # xdmf_file.read(mesh_test, "U")
-
-
-                    # xdmf_filename = self.initialize_U_from_file+"_"+str(k_frame).zfill(2)+".h5"
-                    # print xdmf_filename
-                    # mesh_test = dolfin.Mesh()
-                    # xdmf_file = dolfin.HDF5File(mesh_test.mpi_comm(),xdmf_filename, 'r')
-                    # xdmf_file.read(self.problem.U, "U", False)
-
                     filename = self.initialize_U_from_file+"_"+str(k_frame).zfill(2)+".vtu"
                     ugrid = myvtk.readUGrid(filename)
                     vtk_u = ugrid.GetPointData().GetArray("U")
-
-                    # print 'get cells', ugrid.GetCells()
-                    # print 'node 0', ugrid.GetCell(0)
-                    # print 'node 1', ugrid.GetCell(1)
-
-                    print 'vtk nb cells = ', ugrid.GetNumberOfCells()
-                    print 'dolfin nb cells = ', self.problem.mesh.num_cells()
-                    print 'vtk nb points = ', ugrid.GetNumberOfPoints()
-                    print 'dolfin nb points = ', self.problem.mesh.num_vertices()
-                    assert ugrid.GetNumberOfCells() == self.problem.mesh.num_cells()
-                    assert ugrid.GetNumberOfPoints() == self.problem.mesh.num_vertices()
-
-                    dofmap = self.problem.U_fs.dofmap()
-                    n_cells = self.problem.mesh.num_cells()
-                    nodes_per_element = dofmap.num_element_dofs(0)
-                    np_connectivity = np.zeros([n_cells, nodes_per_element], dtype=np.int)
-                    for i in range(n_cells):
-                        np_connectivity[i,:] = dofmap.cell_dofs(i)
-                    PERM_DOLFIN_TO_VTK = [0, 1, 2, 3]
-                    np_connectivity = np_connectivity[:, PERM_DOLFIN_TO_VTK]
-                    np_connectivity = np.hstack((np.ones([n_cells, 1], dtype=np.int)*nodes_per_element,
-                                                 np_connectivity)).flatten()
-                    vtk_connectivity = ns.numpy_to_vtkIdTypeArray(np_connectivity)
-                    vtk_cells = vtk.vtkCellArray()
-                    vtk_cells.SetCells(n_cells, vtk_connectivity)
-                    print 'vtk cells', vtk_cells
-                    vtk_grid = vtk.vtkUnstructuredGrid()
-                    vtk_grid.SetCells(vtk.VTK_QUADRATIC_TETRA, vtk_cells)
-                    print 'get cells', vtk_grid.GetCells()
-                    # print 'node 0', vtk_grid.GetCell(0)
-                    print 'node 1', vtk_grid.GetCell(1)
-
-                    # assert vtk_grid.GetCells() == ugrid.GetCells()
-                    # assert vtk_grid.GetCell(0) == ugrid.GetCell(0)
-
-
-
-                    self.problem.U = ns.vtk_to_numpy(vtk_u)
-
+                    np_U = ns.vtk_to_numpy(vtk_u)
+                    np_U_float = np_U.astype(float)
+                    np_U_float_reshape = np.reshape(np_U_float, np_U_float.size)
+                    self.problem.U.vector()[:] = np_U_float_reshape[dolfin.dof_to_vertex_map(self.problem.U_fs)]
 
                 elif (self.initialize_DU_with_DUold):
                     self.problem.U.vector().axpy(1., self.problem.DUold.vector())