Merge branch 'GIMIC'

Energy.py

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

Energy_Regularization.py

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

Energy_WarpedImage.py

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

ImageIterator.py

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

ImageSeries.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 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

NonlinearSolver.py

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

Problem.py

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

Problem_ImageRegistration.py

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

__init__.py

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

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

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)
 

fedic2.py

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

plot_binned_strains_vs_radius.py

@@ -8,8 +8,6 @@
 ###                                                                          ###
 ################################################################################
 
-import math
-import numpy
 import os
 import sys
 

plot_displacement_error.py

@@ -8,8 +8,6 @@
 ###                                                                          ###
 ################################################################################
 
-import math
-import numpy
 import os
 import sys
 

write_VTU_file.py

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