diff --git a/Energy.py b/Energy.py
new file mode 100644
index 0000000000000000000000000000000000000000..a645ef42bb0639043bb6ac83b2253ad86360a0c6
--- /dev/null
+++ b/Energy.py
@@ -0,0 +1,59 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2018 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import dolfin
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+class Energy():
+
+
+
+ def reinit(self,
+ *args,
+ **kwargs):
+
+ pass
+
+
+
+ def call_before_solve(self,
+ *args,
+ **kwargs):
+
+ pass
+
+
+
+ def call_after_solve(self,
+ *args,
+ **kwargs):
+
+ pass
+
+
+
+ def get_qoi_names(self):
+
+ return [self.name+"_ener"]
+
+
+
+ def get_qoi_values(self):
+
+ self.ener = (dolfin.assemble(self.ener_form)/self.problem.mesh_V0)**(1./2)
+ self.printer.print_sci(self.name+"_ener",self.ener)
+
+ return [self.ener]
diff --git a/Energy_Regularization.py b/Energy_Regularization.py
new file mode 100644
index 0000000000000000000000000000000000000000..651d89d813210528acfebed45ec15370954f6b26
--- /dev/null
+++ b/Energy_Regularization.py
@@ -0,0 +1,142 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2018 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import dolfin
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_cm as dcm
+
+import dolfin_dic as ddic
+from Energy import Energy
+
+################################################################################
+
+class RegularizationEnergy(Energy):
+
+
+
+ def __init__(self,
+ problem,
+ name="reg",
+ w=1.,
+ type="equilibrated",
+ model="neohookean",
+ young=1.,
+ poisson=0.,
+ quadrature_degree=None):
+
+ self.problem = problem
+ self.printer = problem.printer
+ self.name = name
+ self.w = w
+ self.type = type
+ self.model = model
+ self.young = young
+ self.poisson = poisson
+ self.quadrature_degree = quadrature_degree
+
+ self.printer.print_str("Defining regularization energy…")
+ self.printer.inc()
+
+ self.printer.print_str("Defining measures…")
+
+ self.form_compiler_parameters = {
+ "representation":"uflacs", # MG20180327: Is that needed?
+ "quadrature_degree":self.quadrature_degree}
+ self.dV = dolfin.Measure(
+ "dx",
+ domain=self.problem.mesh,
+ metadata=self.form_compiler_parameters)
+ self.dF = dolfin.Measure(
+ "dS",
+ domain=self.problem.mesh,
+ metadata=self.form_compiler_parameters)
+ self.dS = dolfin.Measure(
+ "ds",
+ domain=self.problem.mesh,
+ metadata=self.form_compiler_parameters)
+
+ self.printer.print_str("Defining mechanical model…")
+
+ self.E = dolfin.Constant(self.young)
+ self.nu = dolfin.Constant(self.poisson)
+ self.material_parameters = {
+ "E":self.E,
+ "nu":self.nu}
+
+ if (self.model == "linear"): # <- super bad
+ self.e = dolfin.sym(dolfin.grad(self.problem.U))
+ self.material = dcm.HookeElasticMaterial(
+ parameters=self.material_parameters)
+ self.Psi_m, self.S_m = self.material.get_free_energy(
+ epsilon=self.e)
+ self.P_m = self.S_m
+ elif (self.model in ("kirchhoff", "neohookean")):
+ self.dim = self.problem.U.ufl_shape[0]
+ self.I = dolfin.Identity(self.dim)
+ self.F = self.I + dolfin.grad(self.problem.U)
+ self.C = self.F.T * self.F
+ if (self.model == "kirchhoff"): # <- pretty bad too
+ self.E = (self.C - self.I)/2
+ self.material = dcm.KirchhoffElasticMaterial(
+ parameters=self.material_parameters)
+ self.Psi_m, self.S_m = self.material.get_free_energy(
+ E=self.E)
+ elif (self.model == "neohookean"):
+ self.material = dcm.CiarletGeymonatNeoHookeanElasticMaterial(
+ parameters=self.material_parameters)
+ self.Psi_m, self.S_m = self.material.get_free_energy(
+ C=self.C)
+ self.P_m = self.F * self.S_m
+ else:
+ assert (0), "\"model\" ("+str(self.model)+") must be \"linear\", \"kirchhoff\" or \"neohookean\". Aborting."
+
+ self.printer.print_str("Defining regularization energy…")
+
+ if (self.type == "hyperelastic"):
+ self.Psi_m_V = self.Psi_m
+ self.Psi_m_F = dolfin.Constant(0)
+ self.Psi_m_S = dolfin.Constant(0)
+ elif (self.type == "equilibrated"):
+ self.Div_P = dolfin.div(self.P_m)
+ self.Psi_m_V = dolfin.inner(self.Div_P,
+ self.Div_P)
+ self.N = dolfin.FacetNormal(self.problem.mesh)
+ self.Jump_P_N = dolfin.jump(self.P_m,
+ self.N)
+ self.cell_h = dolfin.Constant(self.problem.mesh.hmin())
+ self.Psi_m_F = dolfin.inner(self.Jump_P_N,
+ self.Jump_P_N)/self.cell_h
+ # self.P_N = self.P_m * self.N
+ # self.P_N_N = dolfin.dot(self.N,
+ # self.P_N)
+ # self.P_N_T = self.P_N - self.P_N_N * self.N
+ # self.Psi_m_S = dolfin.inner(self.P_N_T,
+ # self.P_N_T)/self.cell_h
+ # self.Psi_m_S = dolfin.inner(self.P_N,
+ # self.P_N)/self.cell_h
+ self.Psi_m_S = dolfin.Constant(0)
+ else:
+ assert (0), "\"type\" ("+str(self.type)+") must be \"hyperelastic\" or \"equilibrated\". Aborting."
+
+ self.DPsi_m_V = dolfin.derivative( self.Psi_m_V, self.problem.U, self.problem.dU_test )
+ self.DPsi_m_F = dolfin.derivative( self.Psi_m_F, self.problem.U, self.problem.dU_test )
+ self.DPsi_m_S = dolfin.derivative( self.Psi_m_S, self.problem.U, self.problem.dU_test )
+ self.DDPsi_m_V = dolfin.derivative(self.DPsi_m_V, self.problem.U, self.problem.dU_trial)
+ self.DDPsi_m_F = dolfin.derivative(self.DPsi_m_F, self.problem.U, self.problem.dU_trial)
+ self.DDPsi_m_S = dolfin.derivative(self.DPsi_m_S, self.problem.U, self.problem.dU_trial)
+
+ self.ener_form = self.Psi_m_V * self.dV + self.Psi_m_F * self.dF + self.Psi_m_S * self.dS
+ self.res_form = self.DPsi_m_V * self.dV + self.DPsi_m_F * self.dF + self.DPsi_m_S * self.dS
+ self.jac_form = self.DDPsi_m_V * self.dV + self.DDPsi_m_F * self.dF + self.DDPsi_m_S * self.dS
+
+ self.printer.dec()
diff --git a/Energy_WarpedImage.py b/Energy_WarpedImage.py
new file mode 100644
index 0000000000000000000000000000000000000000..454865a43e676fe073247e1a9de1c554daeead63
--- /dev/null
+++ b/Energy_WarpedImage.py
@@ -0,0 +1,283 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2018 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import dolfin
+import numpy
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+from Energy import Energy
+
+################################################################################
+
+class WarpedImageEnergy(Energy):
+
+
+
+ def __init__(self,
+ problem,
+ image_series,
+ quadrature_degree,
+ name="im",
+ w=1.,
+ ref_frame=0,
+ dynamic_scaling=False):
+
+ self.problem = problem
+ self.printer = self.problem.printer
+ self.image_series = image_series
+ self.quadrature_degree = quadrature_degree
+ self.name = name
+ self.w = w
+ self.ref_frame = ref_frame
+ self.dynamic_scaling = dynamic_scaling
+
+ self.printer.print_str("Defining warped image correlation energy…")
+ self.printer.inc()
+
+ self.printer.print_str("Defining quadrature finite elements…")
+
+ # fe
+ self.fe = dolfin.FiniteElement(
+ family="Quadrature",
+ cell=self.problem.mesh.ufl_cell(),
+ degree=self.quadrature_degree,
+ quad_scheme="default")
+ self.fe._quad_scheme = "default" # should not be needed
+ for sub_element in self.fe.sub_elements(): # should not be needed
+ sub_element._quad_scheme = "default" # should not be needed
+
+ # ve
+ self.ve = dolfin.VectorElement(
+ family="Quadrature",
+ cell=self.problem.mesh.ufl_cell(),
+ degree=self.quadrature_degree,
+ quad_scheme="default")
+ self.ve._quad_scheme = "default" # should not be needed
+ for sub_element in self.ve.sub_elements(): # should not be needed
+ sub_element._quad_scheme = "default" # should not be needed
+
+ # te
+ self.te = dolfin.TensorElement(
+ family="Quadrature",
+ cell=self.problem.mesh.ufl_cell(),
+ degree=self.quadrature_degree,
+ quad_scheme="default")
+ self.te._quad_scheme = "default" # should not be needed
+ for sub_element in self.te.sub_elements(): # should not be needed
+ sub_element._quad_scheme = "default" # should not be needed
+
+ self.printer.print_str("Defining measure…")
+
+ self.form_compiler_parameters = {
+ "quadrature_degree":self.quadrature_degree,
+ "quadrature_scheme":"default"}
+ self.dV = dolfin.Measure(
+ "dx",
+ domain=self.problem.mesh,
+ metadata=self.form_compiler_parameters)
+
+ self.printer.print_str("Loading reference image…")
+ self.printer.inc()
+
+ # ref_frame
+ assert (abs(self.ref_frame) < self.image_series.n_frames),\
+ "abs(ref_frame) = "+str(abs(self.ref_frame))+" >= "+str(self.image_series.n_frames)+" = image_series.n_frames. Aborting."
+ self.ref_frame = self.ref_frame%self.image_series.n_frames
+ self.printer.print_var("ref_frame",self.ref_frame)
+
+ # Iref
+ self.Iref = dolfin.Expression(
+ cppcode=ddic.get_ExprIm_cpp(
+ im_dim=self.image_series.dimension,
+ im_type="im",
+ im_is_def=0),
+ element=self.fe)
+ self.ref_image_filename = self.image_series.get_image_filename(self.ref_frame)
+ self.Iref.init_image(self.ref_image_filename)
+
+ self.Iref_int = dolfin.assemble(self.Iref * self.dV)/self.problem.mesh_V0
+ self.printer.print_var("Iref_int",self.Iref_int)
+
+ self.Iref_norm = (dolfin.assemble(self.Iref**2 * self.dV)/self.problem.mesh_V0)**(1./2)
+ assert (self.Iref_norm > 0.),\
+ "Iref_norm = "+str(self.Iref_norm)+" <= 0. Aborting."
+ self.printer.print_var("Iref_norm",self.Iref_norm)
+
+ # DIref
+ self.DIref = dolfin.Expression(
+ cppcode=ddic.get_ExprIm_cpp(
+ im_dim=self.image_series.dimension,
+ im_type="grad" if (self.image_series.grad_basename is None) else "grad_no_deriv",
+ im_is_def=0),
+ element=self.ve)
+ self.ref_image_grad_filename = self.image_series.get_image_grad_filename(self.ref_frame)
+ self.DIref.init_image(self.ref_image_grad_filename)
+
+ self.printer.dec()
+ self.printer.print_str("Defining deformed image…")
+
+ self.scaling = numpy.array([1.,0.])
+ if (self.dynamic_scaling):
+ self.p = numpy.empty((2,2))
+ self.q = numpy.empty(2)
+
+ # Idef
+ self.Idef = dolfin.Expression(
+ cppcode=ddic.get_ExprIm_cpp(
+ im_dim=self.image_series.dimension,
+ im_type="im",
+ im_is_def=1),
+ element=self.fe)
+ self.Idef.init_image(self.ref_image_filename)
+ self.Idef.init_disp(self.problem.U)
+ self.Idef.init_dynamic_scaling(self.scaling)
+
+ # DIdef
+ self.DIdef = dolfin.Expression(
+ cppcode=ddic.get_ExprIm_cpp(
+ im_dim=self.image_series.dimension,
+ im_type="grad" if (self.image_series.grad_basename is None) else "grad_no_deriv",
+ im_is_def=1),
+ element=self.ve)
+ self.DIdef.init_image(self.ref_image_filename)
+ self.DIdef.init_disp(self.problem.U)
+ self.DIdef.init_dynamic_scaling(self.scaling)
+
+ self.printer.print_str("Defining previous image…")
+
+ # Iold
+ self.Iold = dolfin.Expression(
+ cppcode=ddic.get_ExprIm_cpp(
+ im_dim=self.image_series.dimension,
+ im_type="im",
+ im_is_def=1),
+ element=self.fe)
+ self.Iold.init_image(self.ref_image_filename)
+ self.Iold.init_disp(self.problem.Uold)
+ self.Iold.init_dynamic_scaling(self.scaling) # 2016/07/25: ok, same scaling must apply to Idef & Iold…
+
+ # DIold
+ self.DIold = dolfin.Expression(
+ cppcode=ddic.get_ExprIm_cpp(
+ im_dim=self.image_series.dimension,
+ im_type="grad" if (self.image_series.grad_basename is None) else "grad_no_deriv",
+ im_is_def=1),
+ element=self.ve)
+ self.DIold.init_image(self.ref_image_filename)
+ self.DIold.init_disp(self.problem.Uold)
+ self.DIold.init_dynamic_scaling(self.scaling) # 2016/07/25: ok, same scaling must apply to Idef & Iold…
+
+ self.printer.print_str("Defining correlation energy…")
+
+ # Phi_ref
+ self.Phi_Iref = dolfin.Expression(
+ cppcode=ddic.get_ExprCharFuncIm_cpp(
+ im_dim=self.image_series.dimension),
+ element=self.fe)
+ self.Phi_Iref.init_image(self.ref_image_filename)
+
+ # Psi_c
+ self.Psi_c = self.Phi_Iref * (self.Idef - self.Iref)**2/2
+ self.DPsi_c = self.Phi_Iref * (self.Idef - self.Iref) * dolfin.dot(self.DIdef, self.problem.dU_test)
+
+ self.DDPsi_c = self.Phi_Iref * dolfin.dot(self.DIdef, self.problem.dU_trial) * dolfin.dot(self.DIdef, self.problem.dU_test)
+ self.DDPsi_c_old = self.Phi_Iref * dolfin.dot(self.DIold, self.problem.dU_trial) * dolfin.dot(self.DIold, self.problem.dU_test)
+ self.DDPsi_c_ref = self.Phi_Iref * dolfin.dot(self.DIref, self.problem.dU_trial) * dolfin.dot(self.DIref, self.problem.dU_test)
+
+ self.Psi_c_old = self.Phi_Iref * (self.Idef - self.Iold)**2/2
+ self.DPsi_c_old = self.Phi_Iref * (self.Idef - self.Iold) * dolfin.dot(self.DIdef, self.problem.dU_test)
+
+ # forms
+ self.ener_form = self.Psi_c * self.dV
+ self.res_form = self.DPsi_c * self.dV
+ self.jac_form = self.DDPsi_c * self.dV
+
+ self.printer.dec()
+
+
+
+ def reinit(self):
+
+ self.scaling[:] = [1.,0.]
+
+
+
+ def call_before_solve(self,
+ k_frame,
+ k_frame_old):
+
+ self.printer.print_str("Loading deformed image for correlation energy…")
+
+ # Idef
+ self.def_image_filename = self.image_series.get_image_filename(k_frame)
+ self.Idef.init_image(self.def_image_filename)
+
+ # DIdef
+ self.def_grad_image_filename = self.image_series.get_image_grad_filename(k_frame)
+ self.DIdef.init_image(self.def_grad_image_filename)
+
+ self.printer.print_str("Loading previous image for correlation energy…")
+
+ # Iold
+ self.old_image_filename = self.image_series.get_image_filename(k_frame_old)
+ self.Iold.init_image(self.old_image_filename)
+
+ # DIold
+ self.old_grad_image_filename = self.image_series.get_image_grad_filename(k_frame_old)
+ self.DIold.init_image(self.old_grad_image_filename)
+
+
+
+ def call_after_solve(self):
+
+ if (self.dynamic_scaling):
+ self.printer.print_str("Updating dynamic scaling…")
+ self.printer.inc()
+
+ self.get_qoi_values()
+
+ self.p[0,0] = dolfin.assemble(self.Idef**2 * self.dV)
+ self.p[0,1] = dolfin.assemble(self.Idef * self.dV)
+ self.p[1,0] = self.p[0,1]
+ self.p[1,1] = 1.
+ self.q[0] = dolfin.assemble(self.Idef*self.Iref * self.dV)
+ self.q[1] = dolfin.assemble(self.Iref * self.dV)
+ self.scaling[:] = numpy.linalg.solve(self.p, self.q)
+ self.printer.print_var("scaling",self.scaling)
+
+ self.Idef.update_dynamic_scaling(self.scaling) # should not be needed
+ self.DIdef.update_dynamic_scaling(self.scaling) # should not be needed
+
+ self.Iold.update_dynamic_scaling(self.scaling) # should not be needed
+ self.DIold.update_dynamic_scaling(self.scaling) # should not be needed
+
+ self.get_qoi_values()
+
+ self.printer.dec()
+
+
+
+ def get_qoi_names(self):
+
+ return [self.name+"_ener", self.name+"_err"]
+
+
+
+ def get_qoi_values(self):
+
+ self.ener = (dolfin.assemble(self.ener_form)/self.problem.mesh_V0)**(1./2)
+ self.printer.print_sci(self.name+"_ener",self.ener)
+ self.err = self.ener/self.Iref_norm
+ self.printer.print_sci(self.name+"_err",self.err)
+
+ return [self.ener, self.err]
diff --git a/ImageIterator.py b/ImageIterator.py
new file mode 100644
index 0000000000000000000000000000000000000000..08dd004efe5cf9bd61341f3ed2eeb04d4c4193c2
--- /dev/null
+++ b/ImageIterator.py
@@ -0,0 +1,155 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2018 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import glob
+import os
+import time
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+class ImageIterator():
+
+
+
+ def __init__(self,
+ problem,
+ solver,
+ parameters={}):
+
+ self.problem = problem
+ self.printer = self.problem.printer
+ self.solver = solver
+
+ 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_DU_with_DUold = parameters["initialize_DU_with_DUold"] if ("initialize_DU_with_DUold" in parameters) else False
+
+
+
+ def iterate(self):
+
+ self.printer.print_str("Writing initial solution…")
+ self.printer.inc()
+
+ if not os.path.exists(self.working_folder):
+ os.mkdir(self.working_folder)
+ pvd_basename = self.working_folder+"/"+self.working_basename
+ for vtu_filename in glob.glob(pvd_basename+"_[0-9]*.vtu"):
+ os.remove(vtu_filename)
+ ddic.write_VTU_file(
+ filebasename=pvd_basename,
+ function=self.problem.U,
+ time=self.problem.images_ref_frame)
+
+ self.printer.dec()
+ self.printer.print_str("Initializing QOI file…")
+ self.printer.inc()
+
+ qoi_names = ["k_frame"]+self.problem.get_qoi_names()
+ qoi_filebasename = self.working_folder+"/"+self.working_basename+"-qoi"
+ qoi_printer = mypy.DataPrinter(
+ names=qoi_names,
+ filename=qoi_filebasename+".dat")
+ qoi_values = [self.problem.images_ref_frame]+self.problem.get_qoi_values()
+ qoi_printer.write_line(
+ values=qoi_values)
+
+ self.printer.dec()
+ self.printer.print_str("Looping over frames…")
+
+ n_iter_tot = 0
+ global_success = True
+ for forward_or_backward in ["forward","backward"]:
+ self.printer.print_var("forward_or_backward",forward_or_backward)
+
+ if (forward_or_backward == "forward"):
+ k_frames = range(self.problem.images_ref_frame+1, self.problem.images_n_frames, +1)
+ elif (forward_or_backward == "backward"):
+ k_frames = range(self.problem.images_ref_frame-1, -1, -1)
+ #self.printer.print_var("k_frames",k_frames)
+
+ if (forward_or_backward == "backward"):
+ self.problem.reinit()
+
+ self.printer.inc()
+ success = True
+ for k_frame in k_frames:
+ self.printer.print_var("k_frame",k_frame,-1)
+
+ if (forward_or_backward == "forward"):
+ k_frame_old = k_frame-1
+ elif (forward_or_backward == "backward"):
+ k_frame_old = k_frame+1
+ #self.printer.print_var("k_frame_old",k_frame_old,-1)
+
+ if (self.initialize_DU_with_DUold):
+ self.problem.U.vector().axpy(1., self.problem.DUold.vector())
+
+ self.problem.call_before_solve(
+ k_frame=k_frame,
+ k_frame_old=k_frame_old)
+
+ self.printer.print_str("Running registration…")
+
+ success, n_iter = self.solver.solve(
+ k_frame=k_frame)
+ n_iter_tot += n_iter
+
+ if not (success):
+ global_success = False
+ break
+
+ self.problem.call_after_solve()
+
+ self.printer.print_str("Writing solution…")
+ self.printer.inc()
+
+ ddic.write_VTU_file(
+ filebasename=pvd_basename,
+ function=self.problem.U,
+ time=k_frame)
+
+ self.printer.dec()
+ self.printer.print_str("Writing QOI file…")
+ self.printer.inc()
+
+ qoi_printer.write_line(
+ [k_frame]+self.problem.get_qoi_values())
+
+ self.printer.dec()
+
+ self.printer.dec()
+
+ if not (global_success):
+ break
+
+ self.printer.print_str("Image iterator finished…")
+ self.printer.inc()
+
+ self.printer.print_var("n_iter_tot",n_iter_tot)
+
+ self.printer.dec()
+ self.printer.print_str("Plotting QOI…")
+
+ qoi_printer.close()
+ commandline = "gnuplot -e \"set terminal pdf;"
+ commandline += " set output '"+qoi_filebasename+".pdf';"
+ commandline += " set grid;"
+ for k_qoi in xrange(1,len(qoi_names)):
+ commandline += " plot '"+qoi_filebasename+".dat' u 1:"+str(1+k_qoi)+" lw 3 title '"+qoi_names[k_qoi]+"';"
+ commandline += "\""
+ os.system(commandline)
+
+ return global_success
diff --git a/ImageSeries.py b/ImageSeries.py
new file mode 100644
index 0000000000000000000000000000000000000000..89eea3b992c6fd2b86b054d77f3fe6051043f2ef
--- /dev/null
+++ b/ImageSeries.py
@@ -0,0 +1,93 @@
+#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 ImageSeries():
+
+
+
+ def __init__(self,
+ problem,
+ folder,
+ basename,
+ grad_folder=None,
+ grad_basename=None,
+ n_frames=None,
+ ext="vti"):
+
+ self.problem = problem
+ self.printer = self.problem.printer
+ self.folder = folder
+ self.basename = basename
+ self.grad_folder = grad_folder
+ self.grad_basename = grad_basename
+ self.n_frames = n_frames
+ self.ext = ext
+
+ self.printer.print_str("Reading image series…")
+ self.printer.inc()
+
+ self.filenames = glob.glob(self.folder+"/"+self.basename+"_[0-9]*"+"."+self.ext)
+ assert (len(self.filenames) >= 2),\
+ "Not enough images ("+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 >= 2),\
+ "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])
+
+ if (self.grad_basename is not None):
+ if (self.grad_folder is None):
+ self.grad_folder = self.folder
+ self.grad_filenames = glob.glob(self.grad_folder+"/"+self.grad_basename+"_[0-9]*"+"."+self.ext)
+ assert (len(self.grad_filenames) >= self.n_frames)
+
+ image = myvtk.readImage(
+ filename=self.get_image_filename(
+ k_frame=0),
+ verbose=0)
+ self.dimension = myvtk.getImageDimensionality(
+ image=image,
+ verbose=0)
+ self.printer.print_var("dimension",self.dimension)
+
+ self.printer.dec()
+
+
+
+ def get_image_filename(self,
+ k_frame):
+
+ return self.folder+"/"+self.basename+"_"+str(k_frame).zfill(self.zfill)+"."+self.ext
+
+
+
+ def get_image_grad_filename(self,
+ k_frame):
+
+ if (self.grad_basename is None):
+ return self.get_image_filename(k_frame)
+ else:
+ return self.grad_folder+"/"+self.grad_basename+"_"+str(k_frame).zfill(self.zfill)+"."+self.ext
diff --git a/NonlinearSolver.py b/NonlinearSolver.py
new file mode 100644
index 0000000000000000000000000000000000000000..7bba64cf5312416c168934068df05ea97f256378
--- /dev/null
+++ b/NonlinearSolver.py
@@ -0,0 +1,358 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2018 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import dolfin
+import glob
+import math
+import numpy
+import os
+import time
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+class NonlinearSolver():
+
+
+
+ def __init__(self,
+ problem,
+ parameters={}):
+
+ self.problem = problem
+ self.printer = self.problem.printer
+
+ # linear solver
+ self.linear_solver_name = parameters["linear_solver_name"] if ("linear_solver_name" in parameters) else "mumps"
+
+ self.res_vec = dolfin.Vector()
+ self.jac_mat = dolfin.Matrix()
+
+ self.linear_solver = dolfin.LUSolver(
+ self.jac_mat,
+ self.linear_solver_name)
+ self.linear_solver.parameters['report'] = bool(0)
+ self.linear_solver.parameters['reuse_factorization'] = bool(0)
+ self.linear_solver.parameters['same_nonzero_pattern'] = bool(1)
+ self.linear_solver.parameters['symmetric'] = bool(1)
+ self.linear_solver.parameters['verbose'] = bool(0)
+
+ # relaxation
+ self.relax_type = parameters["relax_type"] if ("relax_type" in parameters) else "gss"
+
+ if (self.relax_type == "aitken"):
+ self.compute_relax = self.compute_relax_aitken
+ elif (self.relax_type == "constant"):
+ self.compute_relax = self.compute_relax_constant
+ self.relax_val = parameters["relax"] if ("relax" in parameters) else 1.
+ elif (self.relax_type == "gss"):
+ self.compute_relax = self.compute_relax_gss
+ self.relax_tol = parameters["relax_tol"] if ("relax_tol" in parameters) else 0
+ self.relax_n_iter_max = parameters["relax_n_iter_max"] if ("relax_n_iter_max" in parameters) else 5
+
+ # iterations control
+ self.tol_dU = parameters["tol_dU"] if ("tol_dU" in parameters) else None
+ self.tol_res_rel = parameters["tol_res_rel"] if ("tol_res_rel" in parameters) else None
+ self.n_iter_max = parameters["n_iter_max"] if ("n_iter_max" in parameters) else 32
+
+ # write iterations
+ self.write_iterations = parameters["write_iterations"] if ("write_iterations" in parameters) else False
+
+ if (self.write_iterations):
+ self.working_folder = parameters["working_folder"]
+ self.working_basename = parameters["working_basename"]
+
+ for filename in glob.glob(self.working_folder+"/"+self.working_basename+"-frame=[0-9]*.*"):
+ os.remove(filename)
+
+
+
+ def solve(self,
+ k_frame=None):
+
+ if (self.write_iterations):
+ self.frame_filebasename = self.working_folder+"/"+self.working_basename+"-frame="+str(k_frame).zfill(len(str(self.problem.images_n_frames)))
+
+ self.frame_printer = mypy.DataPrinter(
+ names=["k_iter", "res_norm", "res_err_rel", "relax", "dU_norm", "U_norm", "dU_err"],
+ filename=self.frame_filebasename+".dat")
+
+ ddic.write_VTU_file(
+ filebasename=self.frame_filebasename,
+ function=self.problem.U,
+ time=0)
+
+ self.k_iter = 0
+ self.success = False
+ self.printer.inc()
+ while (True):
+ self.k_iter += 1
+ self.printer.print_var("k_iter",self.k_iter,-1)
+
+ # linear problem
+ self.linear_success = self.linear_solve()
+ if not (self.linear_success):
+ break
+
+ # relaxation
+ self.compute_relax()
+
+ # solution update
+ self.problem.U.vector().axpy(self.relax, self.problem.dU.vector())
+ self.problem.U_norm = self.problem.U.vector().norm("l2")
+ #self.printer.print_sci("U_norm",self.problem.U_norm)
+
+ if (self.write_iterations):
+ ddic.write_VTU_file(
+ filebasename=self.frame_filebasename,
+ function=self.problem.U,
+ time=self.k_iter)
+
+ # displacement error
+ self.problem.dU_norm *= abs(self.relax)
+ #self.printer.print_sci("dU_norm",self.problem.dU_norm)
+ if (self.problem.dU_norm == 0.) and (self.problem.Uold_norm == 0.) and (self.problem.U_norm == 0.):
+ self.problem.dU_err = 0.
+ elif (self.problem.Uold_norm == 0.):
+ self.problem.dU_err = self.problem.dU_norm/self.problem.U_norm
+ else:
+ self.problem.dU_err = self.problem.dU_norm/self.problem.Uold_norm
+ self.printer.print_sci("dU_err",self.problem.dU_err)
+
+ if (self.write_iterations):
+ self.frame_printer.write_line([self.k_iter, self.res_norm, self.res_err_rel, self.relax, self.problem.dU_norm, self.problem.U_norm, self.problem.dU_err])
+
+ # exit test
+ self.success = True
+ if (self.tol_res_rel is not None) and (self.res_err_rel > self.tol_res_rel):
+ self.success = False
+ if (self.tol_dU is not None) and (self.problem.dU_err > self.tol_dU ):
+ self.success = False
+
+ # exit
+ if (self.success):
+ self.printer.print_str("Nonlinear solver converged…")
+ break
+
+ if (self.k_iter == self.n_iter_max):
+ self.printer.print_str("Warning! Nonlinear solver failed to converge… (k_frame = "+str(k_frame)+")")
+ break
+
+ self.printer.dec()
+
+ if (self.write_iterations):
+ self.frame_printer.close()
+ commandline = "gnuplot -e \"set terminal pdf noenhanced;"
+ commandline += " set output '"+self.frame_filebasename+".pdf';"
+ commandline += " set key box textcolor variable;"
+ commandline += " set grid;"
+ commandline += " set logscale y;"
+ commandline += " set yrange [1e-3:1e0];"
+ commandline += " plot '"+self.frame_filebasename+".dat' u 1:7 pt 1 lw 3 title 'dU_err', "+str(self.tol_dU)+" lt -1 notitle;"
+ commandline += " unset logscale y;"
+ commandline += " set yrange [*:*];"
+ commandline += " plot '' u 1:4 pt 1 lw 3 title 'relax'\""
+ os.system(commandline)
+
+ return self.success, self.k_iter
+
+
+
+ def linear_solve(self):
+
+ # res_old
+ if (self.k_iter > 1):
+ if (hasattr(self, "res_old_vec")):
+ self.res_old_vec[:] = self.res_vec[:]
+ else:
+ self.res_old_vec = self.res_vec.copy()
+ self.res_old_norm = self.res_norm
+
+ # linear system: residual assembly
+ self.printer.print_str("Residual assembly…",newline=False)
+ timer = time.time()
+ self.problem.assemble_res(
+ res_vec=self.res_vec)
+ timer = time.time() - timer
+ self.printer.print_str(" "+str(timer)+" s",tab=False)
+
+ self.printer.inc()
+
+ # res_norm
+ self.res_norm = self.res_vec.norm("l2")
+ #self.printer.print_sci("res_norm",self.res_norm)
+
+ # dres
+ if (self.k_iter > 1):
+ if (hasattr(self, "dres_vec")):
+ self.dres_vec[:] = self.res_vec[:] - self.res_old_vec[:]
+ else:
+ self.dres_vec = self.res_vec - self.res_old_vec
+ self.dres_norm = self.dres_vec.norm("l2")
+ self.printer.print_sci("dres_norm",self.dres_norm)
+
+ # res_err_rel
+ if (self.k_iter == 1):
+ self.res_err_rel = 1.
+ else:
+ self.res_err_rel = self.dres_norm / self.res_old_norm
+ self.printer.print_sci("res_err_rel",self.res_err_rel)
+
+ self.printer.dec()
+
+ # linear system: matrix assembly
+ self.printer.print_str("Jacobian assembly…",newline=False)
+ timer = time.time()
+ self.problem.assemble_jac(
+ jac_mat=self.jac_mat)
+ timer = time.time() - timer
+ self.printer.print_str(" "+str(timer)+" s",tab=False)
+
+ # linear system: solve
+ try:
+ self.printer.print_str("Solve…",newline=False)
+ timer = time.time()
+ self.linear_solver.solve(
+ self.problem.dU.vector(),
+ self.res_vec)
+ timer = time.time() - timer
+ self.printer.print_str(" "+str(timer)+" s",tab=False)
+ except:
+ self.printer.print_str("Warning! Linear solver failed!",tab=False)
+ return False
+ #self.printer.print_var("dU",dU.vector().array())
+
+ self.printer.inc()
+
+ # dU_norm
+ self.problem.dU_norm = self.problem.dU.vector().norm("l2")
+ #self.printer.print_sci("dU_norm",self.problem.dU_norm)
+ if not (numpy.isfinite(self.problem.dU_norm)):
+ self.printer.print_str("Warning! Solution increment is NaN! Setting it to 0.",tab=False)
+ self.problem.dU.vector().zero()
+
+ self.printer.dec()
+
+ return True
+
+
+
+ def compute_relax_aitken(self):
+
+ if (self.k_iter == 1):
+ self.relax = 1.
+ else:
+ self.relax *= (-1.) * self.res_old_vec.inner(self.dres_vec) / self.dres_norm**2
+ self.printer.print_sci("relax",self.relax)
+
+
+
+ def compute_relax_constant(self):
+
+ self.relax = self.relax_val
+ self.printer.print_sci("relax",self.relax)
+
+
+
+ def compute_relax_gss(self):
+
+ phi = (1 + math.sqrt(5)) / 2
+ relax_a = (1-phi)/(2-phi)
+ relax_b = 1/(2-phi)
+ need_update_c = True
+ need_update_d = True
+ relax_cur = 0.
+ relax_list = []
+ ener_list = []
+ self.printer.inc()
+ relax_k = 1
+ while (True):
+ self.printer.print_var("relax_k",relax_k,-1)
+ # self.printer.print_sci("relax_a",relax_a)
+ # self.printer.print_sci("relax_b",relax_b)
+ if (need_update_c):
+ relax_c = relax_b - (relax_b - relax_a) / phi
+ relax_list.append(relax_c)
+ self.printer.print_sci("relax_c",relax_c)
+ self.problem.U.vector().axpy(relax_c-relax_cur, self.problem.dU.vector())
+ relax_cur = relax_c
+ relax_fc = self.problem.assemble_ener()
+ #self.printer.print_sci("relax_fc",relax_fc)
+ if (numpy.isnan(relax_fc)):
+ relax_fc = float('+inf')
+ #self.printer.print_sci("relax_fc",relax_fc)
+ self.printer.print_sci("relax_fc",relax_fc)
+ ener_list.append(relax_fc)
+ if (need_update_d):
+ relax_d = relax_a + (relax_b - relax_a) / phi
+ relax_list.append(relax_d)
+ self.printer.print_sci("relax_d",relax_d)
+ self.problem.U.vector().axpy(relax_d-relax_cur, self.problem.dU.vector())
+ relax_cur = relax_d
+ relax_fd = self.problem.assemble_ener()
+ if (numpy.isnan(relax_fd)):
+ relax_fd = float('+inf')
+ #self.printer.print_sci("relax_fd",relax_fd)
+ self.printer.print_sci("relax_fd",relax_fd)
+ ener_list.append(relax_fd)
+ # self.printer.print_var("relax_list",relax_list)
+ # self.printer.print_var("ener_list",ener_list)
+ if (relax_k > 1):
+ ener_min_old = ener_min
+ relax_min = relax_list[numpy.argmin(ener_list)]
+ # self.printer.print_var("relax_min",relax_min)
+ ener_min = min(ener_list)
+ # self.printer.print_var("ener_min",ener_min)
+ if (relax_k > 1):
+ dener_min = ener_min-ener_min_old
+ self.printer.print_var("dener_min",dener_min)
+ relax_err = dener_min/ener_list[0]
+ self.printer.print_var("relax_err",relax_err)
+ if (relax_min != 0.) and (relax_err < self.relax_tol):
+ break
+ if (relax_k == self.relax_n_iter_max):
+ break
+ if (relax_fc < relax_fd):
+ relax_b = relax_d
+ relax_d = relax_c
+ relax_fd = relax_fc
+ need_update_c = True
+ need_update_d = False
+ elif (relax_fc >= relax_fd):
+ relax_a = relax_c
+ relax_c = relax_d
+ relax_fc = relax_fd
+ need_update_c = False
+ need_update_d = True
+ else: assert(0)
+ relax_k += 1
+ self.printer.dec()
+ self.problem.U.vector().axpy(-relax_cur, self.problem.dU.vector())
+ #self.printer.print_var("ener_list",ener_list)
+
+ if (self.write_iterations):
+ self.iter_filebasename = self.frame_filebasename+"-iter="+str(self.k_iter).zfill(3)
+ file_dat_iter = open(self.iter_filebasename+".dat","w")
+ file_dat_iter.write("\n".join([" ".join([str(val) for val in [relax_list[relax_k], ener_list[relax_k]]]) for relax_k in xrange(len(relax_list))]))
+ file_dat_iter.close()
+ commandline = "gnuplot -e \"set terminal pdf;"
+ commandline += " set output '"+self.iter_filebasename+".pdf';"
+ commandline += " plot '"+self.iter_filebasename+".dat' using 1:2 with points title 'psi_int';"
+ commandline += " plot '"+self.iter_filebasename+".dat' u 1:2 with points title 'psi_int', '"+self.iter_filebasename+".dat' using (\$2=='inf'?\$1:1/0):(GPVAL_Y_MIN+(0.8)*(GPVAL_Y_MAX-GPVAL_Y_MIN)):(0):((0.2)*(GPVAL_Y_MAX-GPVAL_Y_MIN)) with vectors notitle\""
+ os.system(commandline)
+
+ self.relax = relax_list[numpy.argmin(ener_list)]
+ self.printer.print_sci("relax",self.relax)
+ if (self.relax == 0.):
+ self.printer.print_str("Warning! Optimal relaxation is null…")
diff --git a/Problem.py b/Problem.py
new file mode 100644
index 0000000000000000000000000000000000000000..cdbba9bc1076127ebeb4b5fd8467c4274c38d03a
--- /dev/null
+++ b/Problem.py
@@ -0,0 +1,20 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2018 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+class Problem():
+
+ pass
diff --git a/Problem_ImageRegistration.py b/Problem_ImageRegistration.py
new file mode 100644
index 0000000000000000000000000000000000000000..fdf6396a8dc1e490f6abd6d46089f2d902bfd691
--- /dev/null
+++ b/Problem_ImageRegistration.py
@@ -0,0 +1,277 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2018 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import dolfin
+import os
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+from Problem import Problem
+
+################################################################################
+
+class ImageRegistrationProblem(Problem):
+
+
+
+ def __init__(self,
+ mesh=None,
+ mesh_folder=None,
+ mesh_basename=None,
+ U_family="Lagrange",
+ U_degree=1):
+
+ self.printer = mypy.Printer()
+
+ self.set_mesh(
+ mesh=mesh,
+ mesh_folder=mesh_folder,
+ mesh_basename=mesh_basename)
+
+ self.set_displacement(
+ U_family=U_family,
+ U_degree=U_degree)
+
+ self.energies = []
+
+
+
+ def __del__(self):
+
+ self.printer.close()
+
+
+
+ def set_mesh(self,
+ mesh=None,
+ mesh_folder=None,
+ mesh_basename=None):
+
+ self.printer.print_str("Loading mesh…")
+ self.printer.inc()
+
+ assert ((mesh is not None) or ((mesh_folder is not None) and (mesh_basename is not None))),\
+ "Must provide a mesh (mesh = "+str(mesh)+") or a mesh file (mesh_folder = "+str(mesh_folder)+", mesh_basename = "+str(mesh_basename)+"). Aborting."
+
+ if (mesh is None):
+ self.mesh_folder = mesh_folder
+ self.mesh_basename = mesh_basename
+ self.mesh_filebasename = self.mesh_folder+"/"+self.mesh_basename
+ self.mesh_filename = self.mesh_filebasename+"."+"xml"
+ assert (os.path.exists(self.mesh_filename)),\
+ "No mesh in "+mesh_filename+". Aborting."
+ self.mesh = dolfin.Mesh(self.mesh_filename)
+ else:
+ self.mesh = mesh
+
+ self.mesh_dimension = self.mesh.ufl_domain().geometric_dimension()
+ assert (self.mesh_dimension in (2,3)),\
+ "mesh_dimension ("+str(self.mesh_dimension)+") must be 2 or 3. Aborting."
+ self.printer.print_var("mesh_dimension",self.mesh_dimension)
+
+ self.printer.print_var("mesh_n_vertices",self.mesh.num_vertices())
+ self.printer.print_var("mesh_n_cells",self.mesh.num_cells())
+
+ self.dV = dolfin.Measure(
+ "dx",
+ domain=self.mesh)
+ self.dS = dolfin.Measure(
+ "ds",
+ domain=self.mesh)
+ self.dF = dolfin.Measure(
+ "dS",
+ domain=self.mesh)
+
+ self.mesh_V0 = dolfin.assemble(dolfin.Constant(1) * self.dV)
+ self.printer.print_sci("mesh_V0",self.mesh_V0)
+ self.mesh_S0 = dolfin.assemble(dolfin.Constant(1) * self.dS)
+ self.printer.print_sci("mesh_S0",self.mesh_S0)
+ self.mesh_h0 = self.mesh_V0**(1./self.mesh_dimension)
+ self.printer.print_sci("mesh_h0",self.mesh_h0)
+ self.mesh_h0 = dolfin.Constant(self.mesh_h0)
+
+ self.printer.dec()
+
+
+
+ def set_displacement(self,
+ U_family="Lagrange",
+ U_degree=1):
+
+ self.printer.print_str("Defining functions…")
+
+ self.U_family = U_family
+ self.U_degree = U_degree
+ self.U_fe = dolfin.VectorElement(
+ family=self.U_family,
+ cell=self.mesh.ufl_cell(),
+ degree=self.U_degree)
+ self.U_fs = dolfin.FunctionSpace(
+ self.mesh,
+ self.U_fe)
+ self.U = dolfin.Function(
+ self.U_fs,
+ name="displacement")
+ self.U.vector().zero()
+ self.U_norm = 0.
+ self.Uold = dolfin.Function(
+ self.U_fs,
+ name="previous displacement")
+ self.Uold.vector().zero()
+ self.Uold_norm = 0.
+ self.DUold = dolfin.Function(
+ self.U_fs,
+ name="previous displacement increment")
+ self.dU = dolfin.Function(
+ self.U_fs,
+ name="displacement correction")
+ self.dU_trial = dolfin.TrialFunction(self.U_fs)
+ self.dU_test = dolfin.TestFunction(self.U_fs)
+
+ # for mesh volume computation
+ self.I = dolfin.Identity(self.mesh_dimension)
+ self.F = self.I + dolfin.grad(self.U)
+ self.J = dolfin.det(self.F)
+
+
+
+ def reinit(self):
+
+ self.U.vector().zero()
+ self.U_norm = 0.
+ self.Uold.vector().zero()
+ self.Uold_norm = 0.
+ self.DUold.vector().zero()
+
+ for energy in self.energies:
+ energy.reinit()
+
+
+
+ def add_image_energy(self,
+ energy):
+
+ if (hasattr(self, "images_n_frames")
+ and hasattr(self, "images_ref_frame")):
+ assert (energy.image_series.n_frames == self.images_n_frames)
+ assert (energy.ref_frame == self.images_ref_frame)
+ else:
+ self.images_n_frames = energy.image_series.n_frames
+ self.images_ref_frame = energy.ref_frame
+
+ self.energies += [energy]
+
+
+
+ def add_regul_energy(self,
+ energy):
+
+ self.energies += [energy]
+
+
+
+ def assemble_ener(self):
+
+ ener_form = 0.
+ for energy in self.energies:
+ ener_form += dolfin.Constant(energy.w) * energy.ener_form
+
+ ener = dolfin.assemble(
+ ener_form)
+ #self.printer.print_var("ener",ener)
+
+ return ener
+
+
+
+ def assemble_res(self,
+ res_vec):
+
+ res_form = 0.
+ for energy in self.energies:
+ res_form -= dolfin.Constant(energy.w) * energy.res_form
+
+ res_vec = dolfin.assemble(
+ res_form,
+ tensor=res_vec)
+ #self.printer.print_var("res_vec",res_vec.array())
+
+
+
+ def assemble_jac(self,
+ jac_mat):
+
+ jac_form = 0.
+ for energy in self.energies:
+ jac_form += dolfin.Constant(energy.w) * energy.jac_form
+
+ jac_mat = dolfin.assemble(
+ jac_form,
+ tensor=jac_mat)
+ #self.printer.print_var("jac_mat",jac_mat.array())
+
+
+
+ def call_before_solve(self,
+ *kargs,
+ **kwargs):
+
+ for energy in self.energies:
+ energy.call_before_solve(
+ *kargs,
+ **kwargs)
+
+
+
+ def call_after_solve(self,
+ *kargs,
+ **kwargs):
+
+ self.DUold.vector()[:] = self.U.vector()[:] - self.Uold.vector()[:]
+ self.Uold.vector()[:] = self.U.vector()[:]
+ self.Uold_norm = self.U_norm
+
+ for energy in self.energies:
+ energy.call_after_solve(
+ *kargs,
+ **kwargs)
+
+
+
+ def get_qoi_names(self):
+
+ names = ["mesh_V"]
+
+ for energy in self.energies:
+ names += energy.get_qoi_names()
+
+ return names
+
+
+
+ def get_qoi_values(self):
+
+ self.compute_mesh_volume()
+ values = [self.mesh_V]
+
+ for energy in self.energies:
+ values += energy.get_qoi_values()
+
+ return values
+
+
+
+ def compute_mesh_volume(self):
+
+ self.mesh_V = dolfin.assemble(self.J * self.dV)
+ self.printer.print_sci("mesh_V",self.mesh_V)
+ return self.mesh_V
diff --git a/__init__.py b/__init__.py
index 4e6273a053f9a6e870475f83f30f4aad186b02db..e67262458e7995999f7626e3578ad9fffdadb3e3 100644
--- a/__init__.py
+++ b/__init__.py
@@ -1,19 +1,29 @@
#this file was generated by __init__.py.py
-from compute_warped_mesh import *
-from image_expressions_py import *
-from compute_warped_images import *
-from compute_unwarped_images import *
from fedic import *
-from generate_images import *
from plot_binned_strains_vs_radius import *
-from plot_regional_strains import *
-from generate_undersampled_images import *
-from compute_strains import *
-from compute_displacements import *
-from compute_displacement_error import *
+from NonlinearSolver import *
+from plot_strains_vs_radius import *
+from compute_unwarped_images import *
+from Problem_ImageRegistration import *
from compute_quadrature_degree import *
-from image_expressions_cpp import *
+from Energy_Regularization import *
+from compute_warped_images import *
from plot_strains import *
+from image_expressions_cpp import *
+from Energy import *
+from compute_warped_mesh import *
+from Energy_WarpedImage import *
+from write_VTU_file import *
+from generate_undersampled_images import *
+from plot_regional_strains import *
+from image_expressions_py import *
+from ImageIterator import *
+from fedic2 import *
from plot_displacement_error import *
-from plot_strains_vs_radius import *
+from Problem import *
+from compute_displacements import *
+from generate_images import *
+from ImageSeries import *
+from compute_displacement_error import *
+from compute_strains import *
diff --git a/compute_quadrature_degree.py b/compute_quadrature_degree.py
index 7a2c2312d658e3734d25176232058aa706314ac8..78755d55ebbca9ecfc145024a437b94374b02cc4 100644
--- a/compute_quadrature_degree.py
+++ b/compute_quadrature_degree.py
@@ -96,7 +96,8 @@ def compute_quadrature_degree_from_points_count(
def compute_quadrature_degree_from_integral(
image_filename,
- mesh,
+ mesh=None,
+ mesh_filebasename=None,
deg_min=1,
deg_max=10,
tol=1e-2,
@@ -111,6 +112,8 @@ def compute_quadrature_degree_from_integral(
verbose=verbose-1)
if (verbose): print "image_dimension = " + str(image_dimension)
+ if (mesh is None):
+ mesh = dolfin.Mesh(mesh_filebasename+"."+"xml")
dX = dolfin.dx(mesh)
first_time = True
diff --git a/fedic.py b/fedic.py
index 833191e77649ef629fb25fd2180b4bbc0cb556d0..ee505f8193732d74952ae7c3d9f1207472a357af 100644
--- a/fedic.py
+++ b/fedic.py
@@ -100,7 +100,7 @@ def fedic(
"images_n_frames = "+str(images_n_frames)+" <= 1. Aborting."
mypy.print_var("images_n_frames",images_n_frames,tab+1)
assert (abs(images_ref_frame) < images_n_frames),\
- "abs(images_ref_frame) = "+str(images_ref_frame)+" >= images_n_frames. Aborting."
+ "abs(images_ref_frame) = "+str(abs(images_ref_frame))+" >= images_n_frames. Aborting."
images_ref_frame = images_ref_frame%images_n_frames
mypy.print_var("images_ref_frame",images_ref_frame,tab+1)
diff --git a/fedic2.py b/fedic2.py
new file mode 100644
index 0000000000000000000000000000000000000000..ffe4b66b185f2c02de545ae6f8c3ee97d58f8f1e
--- /dev/null
+++ b/fedic2.py
@@ -0,0 +1,136 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2018 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+def fedic2(
+ working_folder,
+ working_basename,
+ images_folder,
+ images_basename,
+ images_grad_basename=None,
+ images_ext="vti", # vti, vtk
+ images_n_frames=None,
+ images_ref_frame=0,
+ images_quadrature=None,
+ images_quadrature_from="points_count", # points_count, integral
+ images_expressions_type="cpp", # cpp, py
+ images_dynamic_scaling=1,
+ mesh=None,
+ mesh_folder=None,
+ mesh_basename=None,
+ mesh_degree=1,
+ regul_type="equilibrated", # hyperelastic, equilibrated
+ regul_model="neohookean", # linear, kirchhoff, neohookean, mooneyrivlin
+ regul_quadrature=None,
+ regul_level=0.1,
+ regul_poisson=0.0,
+ tangent_type="Idef", # Idef, Idef-wHess, Iold, Iref
+ residual_type="Iref", # Iref, Iold, Iref-then-Iold
+ relax_type="gss", # constant, aitken, gss
+ relax_init=1.0,
+ initialize_DU_with_DUold=0,
+ tol_res=None,
+ tol_res_rel=None,
+ tol_dU=None,
+ tol_im=None,
+ n_iter_max=100,
+ continue_after_fail=0,
+ print_refined_mesh=0,
+ print_iterations=0):
+
+ assert (images_expressions_type == "cpp"),\
+ "Python image expression are deprecated. Aborting."
+ assert (tangent_type == "Idef"),\
+ "tangent_type must be \"Idef\". Aborting."
+ assert (residual_type == "Iref"),\
+ "residual_type must be \"Iref\". Aborting."
+ assert (relax_init == 1.),\
+ "relax_init must be 1. Aborting."
+ assert (tol_res is None),\
+ "tol_res is deprecated. Aborting."
+ assert (tol_im is None),\
+ "tol_im is deprecated. Aborting."
+ assert (continue_after_fail == 0),\
+ "continue_after_fail is deprecated. Aborting."
+ assert (print_refined_mesh == 0),\
+ "print_refined_mesh is deprecated. Aborting."
+
+ problem = ddic.ImageRegistrationProblem(
+ mesh=mesh,
+ mesh_folder=mesh_folder,
+ mesh_basename=mesh_basename,
+ U_degree=mesh_degree)
+
+ image_series = ddic.ImageSeries(
+ problem=problem,
+ folder=images_folder,
+ basename=images_basename,
+ grad_basename=images_grad_basename,
+ n_frames=images_n_frames,
+ ext=images_ext)
+
+ if (images_quadrature is None):
+ if (images_quadrature_from == "points_count"):
+ images_quadrature = ddic.compute_quadrature_degree_from_points_count(
+ image_filename=image_series.get_image_filename(images_ref_frame),
+ mesh_filebasename=mesh_folder+"/"+mesh_basename,
+ verbose=1)
+ elif (method == "integral"):
+ images_quadrature = ddic.compute_quadrature_degree_from_integral(
+ image_filename=self.get_image_filename(images_ref_frame),
+ mesh_filebasename=mesh_folder+"/"+mesh_basename,
+ verbose=1)
+ else:
+ assert (0), "\"images_quadrature_from\" (="+str(images_quadrature_from)+") must be \"points_count\" or \"integral\". Aborting."
+
+ warped_image_energy = ddic.WarpedImageEnergy(
+ problem=problem,
+ image_series=image_series,
+ quadrature_degree=images_quadrature,
+ w=1.-regul_level,
+ ref_frame=images_ref_frame,
+ dynamic_scaling=images_dynamic_scaling)
+ problem.add_image_energy(warped_image_energy)
+
+ regularization_energy = ddic.RegularizationEnergy(
+ problem=problem,
+ w=regul_level,
+ type=regul_type,
+ model=regul_model,
+ poisson=regul_poisson,
+ quadrature_degree=regul_quadrature)
+ problem.add_regul_energy(regularization_energy)
+
+ solver = ddic.NonlinearSolver(
+ problem=problem,
+ parameters={
+ "working_folder":working_folder,
+ "working_basename":working_basename,
+ "relax_type":relax_type,
+ "tol_res_rel":tol_res_rel,
+ "tol_dU":tol_dU,
+ "n_iter_max":n_iter_max,
+ "write_iterations":print_iterations})
+
+ image_iterator = ddic.ImageIterator(
+ problem=problem,
+ solver=solver,
+ parameters={
+ "working_folder":working_folder,
+ "working_basename":working_basename,
+ "initialize_DU_with_DUold":initialize_DU_with_DUold})
+
+ image_iterator.iterate()
diff --git a/plot_binned_strains_vs_radius.py b/plot_binned_strains_vs_radius.py
index 3b3d6e8240a6d879c2a9249a38305244c32f5cf4..25bfc6fce9bc2c39451dd575ef8f0e7315fc5edd 100644
--- a/plot_binned_strains_vs_radius.py
+++ b/plot_binned_strains_vs_radius.py
@@ -8,8 +8,6 @@
### ###
################################################################################
-import math
-import numpy
import os
import sys
diff --git a/plot_displacement_error.py b/plot_displacement_error.py
index b0007aa42eb874a2bee6d281f1a9e0e42dd63dfe..727e80dc5f294f19f764cf04604f0af4a9d4d017 100644
--- a/plot_displacement_error.py
+++ b/plot_displacement_error.py
@@ -8,8 +8,6 @@
### ###
################################################################################
-import math
-import numpy
import os
import sys
diff --git a/write_VTU_file.py b/write_VTU_file.py
new file mode 100644
index 0000000000000000000000000000000000000000..9a71f67414a5b5c34e6861c7171a0fc721998709
--- /dev/null
+++ b/write_VTU_file.py
@@ -0,0 +1,34 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2018 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import dolfin
+import os
+import shutil
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+def write_VTU_file(
+ filebasename,
+ function,
+ time,
+ zfill=6):
+
+ file_pvd = dolfin.File(filebasename+"__.pvd")
+ file_pvd << (function, float(time))
+ os.remove(
+ filebasename+"__.pvd")
+ shutil.move(
+ filebasename+"__"+"".zfill(zfill)+".vtu",
+ filebasename+"_"+str(time).zfill(zfill)+".vtu")