From 2a65da728392e701bf0b8e91a13c89b477def0f1 Mon Sep 17 00:00:00 2001
From: Martin Genet <martin.genet@polytechnique.edu>
Date: Wed, 16 Jan 2019 21:49:48 +0100
Subject: [PATCH] Using MeshSeries in initialization of displacement

---
 ImageIterator.py | 32 ++++++++++++++++++++------------
 MeshSeries.py    |  9 +++++++++
 __init__.py      | 41 +++++++++++++++++++++--------------------
 fedic2.py        | 12 ++++++++++--
 4 files changed, 60 insertions(+), 34 deletions(-)

diff --git a/ImageIterator.py b/ImageIterator.py
index f88cc8d..3e6fae3 100644
--- a/ImageIterator.py
+++ b/ImageIterator.py
@@ -8,19 +8,17 @@
 ###                                                                          ###
 ################################################################################
 
+import dolfin
 import glob
+import numpy
 import os
 import time
+import vtk
 
 import myPythonLibrary as mypy
 import myVTKPythonLibrary as myvtk
 
 import dolfin_dic as ddic
-import dolfin
-
-from vtk.util import numpy_support as ns
-import numpy as np
-import vtk
 
 ################################################################################
 
@@ -40,6 +38,10 @@ class ImageIterator():
         self.working_folder           = parameters["working_folder"]           if ("working_folder"           in parameters) else "."
         self.working_basename         = parameters["working_basename"]         if ("working_basename"         in parameters) else "sol"
         self.initialize_U_from_file   = parameters["initialize_U_from_file"]   if ("initialize_U_from_file"   in parameters) else False
+        self.initialize_U_folder      = parameters["initialize_U_folder"]      if ("initialize_U_folder"      in parameters) else "."
+        self.initialize_U_basename    = parameters["initialize_U_basename"]    if ("initialize_U_basename"    in parameters) else None
+        self.initialize_U_ext         = parameters["initialize_U_ext"]         if ("initialize_U_ext"         in parameters) else "vtu"
+        self.initialize_U_array_name  = parameters["initialize_U_array_name"]  if ("initialize_U_array_name"  in parameters) else "displacement"
         self.initialize_DU_with_DUold = parameters["initialize_DU_with_DUold"] if ("initialize_DU_with_DUold" in parameters) else False
 
 
@@ -75,6 +77,13 @@ class ImageIterator():
         self.printer.dec()
         self.printer.print_str("Looping over frames…")
 
+        if (self.initialize_U_from_file):
+            mesh_series = ddic.MeshSeries(
+                problem=self.problem,
+                folder=self.initialize_U_folder,
+                basename=self.initialize_U_basename,
+                ext=self.initialize_U_ext)
+
         n_iter_tot = 0
         global_success = True
         for forward_or_backward in ["forward","backward"]:
@@ -101,13 +110,12 @@ class ImageIterator():
                 #self.printer.print_var("k_frame_old",k_frame_old,-1)
 
                 if (self.initialize_U_from_file):
-                    filename = self.initialize_U_from_file+"_"+str(k_frame).zfill(2)+".vtu"
-                    ugrid = myvtk.readUGrid(filename)
-                    vtk_u = ugrid.GetPointData().GetArray("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)]
+                    mesh = mesh_series.get_mesh(k_frame)
+                    array_U = mesh.GetPointData().GetArray(self.initialize_U_array_name)
+                    array_U = vtk.util.numpy_support.vtk_to_numpy(array_U)
+                    array_U = array_U.astype(float)
+                    array_U = numpy.reshape(array_U, array_U.size)
+                    self.problem.U.vector()[:] = array_U[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())
diff --git a/MeshSeries.py b/MeshSeries.py
index 185d909..51a1206 100644
--- a/MeshSeries.py
+++ b/MeshSeries.py
@@ -62,3 +62,12 @@ class MeshSeries():
             k_frame):
 
         return self.folder+"/"+self.basename+"_"+str(k_frame).zfill(self.zfill)+"."+self.ext
+
+
+
+    def get_mesh(self,
+            k_frame):
+
+        return myvtk.readDataSet(
+            filename=self.get_mesh_filename(
+                k_frame=k_frame))
diff --git a/__init__.py b/__init__.py
index 9deddd4..222a8f7 100644
--- a/__init__.py
+++ b/__init__.py
@@ -1,30 +1,31 @@
 #this file was generated by __init__.py.py
 
-from fedic2 import *
-from NonlinearSolver import *
-from compute_warped_mesh import *
-from image_expressions_py import *
-from compute_warped_images import *
+from compute_displacement_error import *
+from compute_displacements import *
+from compute_quadrature_degree import *
+from compute_strains import *
 from compute_unwarped_images import *
+from compute_warped_images import *
+from compute_warped_mesh import *
+from Energy import *
+from Energy_Regularization import *
+from Energy_WarpedImage import *
 from fedic import *
+from fedic2 import *
 from generate_images import *
-from plot_binned_strains_vs_radius import *
-from plot_regional_strains import *
-from ImageIterator import *
-from Problem import *
-from Energy_WarpedImage import *
 from generate_undersampled_images import *
-from Energy_Regularization import *
-from compute_strains import *
-from compute_displacements import *
-from compute_displacement_error import *
-from ImageSeries import *
-from compute_quadrature_degree import *
 from generated_image_expressions_cpp import *
-from write_VTU_file import *
-from Energy import *
 from image_expressions_cpp import *
-from Problem_ImageRegistration import *
-from plot_strains import *
+from image_expressions_py import *
+from ImageIterator import *
+from ImageSeries import *
+from MeshSeries import *
+from NonlinearSolver import *
+from plot_binned_strains_vs_radius import *
 from plot_displacement_error import *
+from plot_regional_strains import *
+from plot_strains import *
 from plot_strains_vs_radius import *
+from Problem import *
+from Problem_ImageRegistration import *
+from write_VTU_file import *
diff --git a/fedic2.py b/fedic2.py
index cc79220..af560c5 100644
--- a/fedic2.py
+++ b/fedic2.py
@@ -42,6 +42,10 @@ def fedic2(
         relax_type="gss", # constant, aitken, gss
         relax_init=1.0,
         initialize_U_from_file=0,
+        initialize_U_folder=None,
+        initialize_U_basename=None,
+        initialize_U_ext="vtu",
+        initialize_U_array_name="displacement",
         initialize_DU_with_DUold=0,
         tol_res=None,
         tol_res_rel=None,
@@ -138,7 +142,11 @@ def fedic2(
         parameters={
             "working_folder":working_folder,
             "working_basename":working_basename,
-            "initialize_DU_with_DUold":initialize_DU_with_DUold,
-            "initialize_U_from_file":initialize_U_from_file})
+            "initialize_U_from_file":initialize_U_from_file,
+            "initialize_U_folder":initialize_U_folder,
+            "initialize_U_basename":initialize_U_basename,
+            "initialize_U_ext":initialize_U_ext,
+            "initialize_U_array_name":initialize_U_array_name,
+            "initialize_DU_with_DUold":initialize_DU_with_DUold})
 
     image_iterator.iterate()
-- 
GitLab