diff --git a/ImageIterator.py b/ImageIterator.py
index 08dd004efe5cf9bd61341f3ed2eeb04d4c4193c2..e650c94a79d970e8f14af73361e1ed66a9e0ff24 100644
--- a/ImageIterator.py
+++ b/ImageIterator.py
@@ -8,9 +8,12 @@
### ###
################################################################################
+import dolfin
import glob
+import numpy
import os
import time
+import vtk
import myPythonLibrary as mypy
import myVTKPythonLibrary as myvtk
@@ -34,6 +37,11 @@ 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
@@ -69,6 +77,15 @@ 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)
+
+ dof_to_vertex_map = dolfin.dof_to_vertex_map(self.problem.U_fs)
+
n_iter_tot = 0
global_success = True
for forward_or_backward in ["forward","backward"]:
@@ -94,7 +111,18 @@ class ImageIterator():
k_frame_old = k_frame+1
#self.printer.print_var("k_frame_old",k_frame_old,-1)
- if (self.initialize_DU_with_DUold):
+ if (self.initialize_U_from_file):
+ 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)[:,:self.problem.mesh_dimension]
+ # print array_U
+ # array_U = array_U.astype(float)
+ # print array_U
+ array_U = numpy.reshape(array_U, array_U.size)
+ # print array_U
+ self.problem.U.vector()[:] = array_U[dof_to_vertex_map]
+
+ elif (self.initialize_DU_with_DUold):
self.problem.U.vector().axpy(1., self.problem.DUold.vector())
self.problem.call_before_solve(
diff --git a/MeshSeries.py b/MeshSeries.py
new file mode 100644
index 0000000000000000000000000000000000000000..e61c69241269a811c66fe8ca755fb85b40856b5c
--- /dev/null
+++ b/MeshSeries.py
@@ -0,0 +1,73 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2018 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import glob
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+
+
+class MeshSeries():
+
+
+
+ def __init__(self,
+ problem,
+ folder,
+ basename,
+ n_frames=None,
+ ext="vtu"):
+
+ self.problem = problem
+ self.printer = self.problem.printer
+ self.folder = folder
+ self.basename = basename
+ self.n_frames = n_frames
+ self.ext = ext
+
+ self.printer.print_str("Reading mesh series…")
+ self.printer.inc()
+
+ self.filenames = glob.glob(self.folder+"/"+self.basename+"_[0-9]*"+"."+self.ext)
+ assert (len(self.filenames) >= 1),\
+ "Not enough meshes ("+self.folder+"/"+self.basename+"_[0-9]*"+"."+self.ext+"). Aborting."
+
+ if (self.n_frames is None):
+ self.n_frames = len(self.filenames)
+ else:
+ assert (self.n_frames <= len(self.filenames))
+ assert (self.n_frames >= 1),\
+ "n_frames = "+str(self.n_frames)+" < 2. Aborting."
+ self.printer.print_var("n_frames",self.n_frames)
+
+ self.zfill = len(self.filenames[0].rsplit("_",1)[-1].split(".",1)[0])
+
+ self.printer.dec()
+
+
+
+ def get_mesh_filename(self,
+ 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 9deddd46220b2a07376a6436b44c404b8e61fc0f..222a8f73c40a23f4d98e50fcb7e42f302b820a07 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 ac4143cb0e71c72ebc5268ee5f0d62b5ed897595..af560c58da46f45725ea79cafa251276ba141983 100644
--- a/fedic2.py
+++ b/fedic2.py
@@ -41,6 +41,11 @@ def fedic2(
residual_type="Iref", # Iref, Iold, Iref-then-Iold
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,
@@ -59,6 +64,8 @@ def fedic2(
"residual_type must be \"Iref\". Aborting."
assert (relax_init == 1.),\
"relax_init must be 1. Aborting."
+ assert (not ((initialize_U_from_file) and (initialize_DU_with_DUold))),\
+ "Cannot initialize U from file and DU with DUold together. Aborting."
assert (tol_res is None),\
"tol_res is deprecated. Aborting."
assert (tol_im is None),\
@@ -135,6 +142,11 @@ def fedic2(
parameters={
"working_folder":working_folder,
"working_basename":working_basename,
+ "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()