Merge branch 'init_disp' into 'master'

Init disp

See merge request mgenet/dolfin_dic!13

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(

MeshSeries.py

+#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))

__init__.py

 #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 *

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()