initial commit

__init__.py

+#this file was generated by __init__.py.py
+
+from fedic import *
+from plot_strains_vs_radius import *
+from plot_strains import *
+from compute_strains import *
+from plot_displacement_error import *
+from compute_displacement_error import *
+from compute_quadrature_degree import *
+from image_expressions import *

__init__.py.py

+#!/usr/bin/python
+#coding=utf8
+
+########################################################################
+###                                                                  ###
+### Created by Martin Genet, 2016                                    ###
+###                                                                  ###
+### École Polytechnique, Palaiseau, France                           ###
+###                                                                  ###
+########################################################################
+
+import os
+
+########################################################################
+
+os.system("rm -rf *.pyc")
+
+init_file = open("__init__.py", "w")
+init_file.write("#this file was generated by __init__.py.py\n\n")
+
+for (dirpath, dirnames, filenames) in os.walk("."):
+    for filename in filenames:
+        if filename.endswith(".py") and not filename.startswith("__init__"):
+            init_file.write("from " + filename[:-3] + " import *\n")
+
+init_file.close()

compute_displacement_error.py

+#coding=utf8
+
+########################################################################
+###                                                                  ###
+### Created by Martin Genet, 2016                                    ###
+###                                                                  ###
+### École Polytechnique, Palaiseau, France                           ###
+###                                                                  ###
+########################################################################
+
+import glob
+import numpy
+
+import myVTKPythonLibrary as myVTK
+
+########################################################################
+
+def compute_displacement_error(
+        sol_folder,
+        sol_basename,
+        ref_folder,
+        ref_basename,
+        sol_ext="vtu",
+        ref_ext="vtu",
+        sol_disp_array_name="displacement",
+        ref_disp_array_name="displacement",
+        verbose=1):
+
+    n_frames = len(glob.glob(ref_folder+"/"+ref_basename+"_*."+ref_ext))
+    assert (len(glob.glob(sol_folder+"/"+sol_basename+"_*."+sol_ext)) == n_frames)
+    if (verbose): print "n_frames = " + str(n_frames)
+
+    ref_zfill = len(glob.glob(ref_folder+"/"+ref_basename+"_*."+ref_ext)[0].rsplit("_")[-1].split(".")[0])
+    sol_zfill = len(glob.glob(sol_folder+"/"+sol_basename+"_*."+sol_ext)[0].rsplit("_")[-1].split(".")[0])
+    if (verbose): print "ref_zfill = " + str(ref_zfill)
+    if (verbose): print "sol_zfill = " + str(sol_zfill)
+
+    error_file = open(sol_folder+"/"+sol_basename+"-error.dat", "w")
+    error_file.write("#t e\n")
+
+    err_int = numpy.empty(n_frames)
+    ref_int = numpy.empty(n_frames)
+    for k_frame in xrange(n_frames):
+        ref = myVTK.readUGrid(
+            filename=ref_folder+"/"+ref_basename+"_"+str(k_frame).zfill(ref_zfill)+"."+ref_ext,
+            verbose=0)
+        n_points = ref.GetNumberOfPoints()
+        n_cells = ref.GetNumberOfCells()
+        sol = myVTK.readUGrid(
+            filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(sol_zfill)+"."+sol_ext,
+            verbose=0)
+        assert (sol.GetNumberOfPoints() == n_points)
+        assert (sol.GetNumberOfCells() == n_cells)
+
+        ref_disp = ref.GetPointData().GetArray(ref_disp_array_name)
+        sol_disp = sol.GetPointData().GetArray(sol_disp_array_name)
+
+        err_int[k_frame] = numpy.sqrt(numpy.mean([numpy.sum([(sol_disp.GetTuple(k_point)[k_dim]-ref_disp.GetTuple(k_point)[k_dim])**2 for k_dim in xrange(3)]) for k_point in xrange(n_points)]))
+        ref_int[k_frame] = numpy.sqrt(numpy.mean([numpy.sum([(ref_disp.GetTuple(k_point)[k_dim])**2 for k_dim in xrange(3)]) for k_point in xrange(n_points)]))
+
+    ref_int_int = numpy.mean(ref_int)
+    err_int_rel = err_int/ref_int_int
+
+    error_file.write("\n".join([" ".join([str(val) for val in [k_frame, err_int[k_frame], ref_int[k_frame], err_int_rel[k_frame]]]) for k_frame in xrange(n_frames)]))
+
+    error_file.close()

compute_quadrature_degree.py

+#coding=utf8
+
+########################################################################
+###                                                                  ###
+### Created by Martin Genet, 2016                                    ###
+###                                                                  ###
+### École Polytechnique, Palaiseau, France                           ###
+###                                                                  ###
+########################################################################
+
+import dolfin
+
+import myFEniCSPythonLibrary as myFEniCS
+
+########################################################################
+
+def compute_quadrature_degree(
+        image_filename,
+        dX,
+        image_dimension=3,
+        deg_min=1,
+        deg_max=10,
+        tol=1e-2,
+        n_under_tol=1,
+        verbose=1):
+
+    first_time = True
+    k_under_tol = 0
+    for degree in xrange(deg_min,deg_max+1):
+        if (verbose): print "degree = " + str(degree)
+        fe = dolfin.FiniteElement(
+            family="Quadrature",
+            degree=degree)
+        if (image_dimension == 2):
+            I0 = myFEniCS.ExprIm2(
+                filename=image_filename,
+                element=fe)
+        elif (image_dimension == 3):
+            I0 = myFEniCS.ExprIm3(
+                filename=image_filename,
+                element=fe)
+        else:
+            assert (0), "image_dimension must be 2 or 3. Aborting."
+        I0_norm = dolfin.assemble(I0**2 * dX)**(1./2)
+        if (verbose): print "I0_norm = " + str(I0_norm)
+        if (first_time):
+            first_time = False
+        else:
+            I0_norm_err = abs(I0_norm-I0_norm_old)/I0_norm_old
+            if (verbose): print "I0_norm_err = " + str(I0_norm_err)
+            if (I0_norm_err < tol):
+                k_under_tol += 1
+            else:
+                k_under_tol = 0
+            if (k_under_tol >= n_under_tol):
+                break
+        I0_norm_old = I0_norm
+    return degree

compute_strains.py

+#coding=utf8
+
+########################################################################
+###                                                                  ###
+### Created by Martin Genet, 2016                                    ###
+###                                                                  ###
+### École Polytechnique, Palaiseau, France                           ###
+###                                                                  ###
+########################################################################
+
+import glob
+import numpy
+
+import myVTKPythonLibrary as myVTK
+
+########################################################################
+
+def compute_strains(
+        sol_folder,
+        sol_basename,
+        sol_ext="vtu",
+        disp_array_name="displacement",
+        ref_folder=None,
+        ref_basename=None,
+        CYL_or_PPS="PPS",
+        write_strain_vs_radius=0,
+        verbose=1):
+
+    if (ref_folder is not None) and (ref_basename is not None):
+        ref_mesh = myVTK.readUGrid(
+            filename=ref_folder+"/"+ref_basename+"-WithLocalBasis.vtk",
+            verbose=0)
+        ref_n_cells = ref_mesh.GetNumberOfCells()
+
+        if (ref_mesh.GetCellData().HasArray("sector_id")):
+            iarray_sector_id = ref_mesh.GetCellData().GetArray("sector_id")
+            n_sector_ids = 0
+            for k_cell in xrange(ref_n_cells):
+                sector_id = int(iarray_sector_id.GetTuple1(k_cell))
+                if (sector_id < 0): continue
+                if (sector_id > n_sector_ids-1):
+                    n_sector_ids = sector_id+1
+            if (verbose): print "n_sector_ids = " + str(n_sector_ids)
+        else:
+            n_sector_ids = 0
+    else:
+        ref_mesh = None
+        n_sector_ids = 0
+
+    n_frames = len(glob.glob(sol_folder+"/"+sol_basename+"_*."+sol_ext))
+    if (verbose): print "n_frames = " + str(n_frames)
+
+    n_zfill = len(glob.glob(sol_folder+"/"+sol_basename+"_*."+sol_ext)[0].rsplit("_")[-1].split(".")[0])
+    if (verbose): print "n_zfill = " + str(n_zfill)
+
+    strain_file = open(sol_folder+"/"+sol_basename+"-strains.dat", "w")
+    strain_file.write("#t Err_avg Err_std Ecc_avg Ecc_std Ell_avg Ell_std Erc_avg Erc_std Erl_avg Erl_std Ecl_avg Ecl_std\n")
+
+    if (write_strain_vs_radius):
+        strain_vs_radius_file = open(sol_folder+"/"+sol_basename+"-strains_vs_radius.dat", "w")
+        strain_vs_radius_file.write("#t r Err Ecc Ell Erc Erl Ecl\n")
+
+    for k_frame in xrange(n_frames):
+        mesh = myVTK.readUGrid(
+            filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(n_zfill)+"."+sol_ext,
+            verbose=0)
+        n_cells = mesh.GetNumberOfCells()
+        if (ref_mesh is not None):
+            assert (n_cells == ref_n_cells)
+            mesh.GetCellData().AddArray(ref_mesh.GetCellData().GetArray("part_id"))
+            mesh.GetCellData().AddArray(ref_mesh.GetCellData().GetArray("sector_id"))
+        myVTK.computeStrainsFromDisplacements(
+            mesh=mesh,
+            disp_array_name=disp_array_name,
+            ref_mesh=ref_mesh,
+            verbose=0)
+        myVTK.writeUGrid(
+            ugrid=mesh,
+            filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(n_zfill)+"."+sol_ext,
+            verbose=0)
+
+        if (ref_mesh is not None):
+            assert (mesh.GetCellData().HasArray("Strain_"+CYL_or_PPS))
+            farray_strain = mesh.GetCellData().GetArray("Strain_"+CYL_or_PPS)
+        else:
+            farray_strain = mesh.GetCellData().GetArray("Strain")
+
+        if (n_sector_ids == 0):
+            strains_all = []
+            for k_cell in xrange(n_cells):
+                strains_all.append(farray_strain.GetTuple(k_cell))
+        elif (n_sector_ids == 1):
+            strains_all = []
+            for k_cell in xrange(n_cells):
+                sector_id = int(iarray_sector_id.GetTuple(k_cell)[0])
+                if (sector_id < 0): continue
+                strains_all.append(farray_strain.GetTuple(k_cell))
+        elif (n_sector_ids > 1):
+            strains_all = []
+            strains_per_sector = [[] for sector_id in xrange(n_sector_ids)]
+            for k_cell in xrange(n_cells):
+                sector_id = int(iarray_sector_id.GetTuple(k_cell)[0])
+                if (sector_id < 0): continue
+                strains_all.append(farray_strain.GetTuple(k_cell))
+                strains_per_sector[sector_id].append(farray_strain.GetTuple(k_cell))
+
+        strain_file.write(str(k_frame))
+        strains_all_avg = numpy.mean(strains_all, 0)
+        strains_all_std = numpy.std(strains_all, 0)
+        strain_file.write("".join([" " + str(strains_all_avg[k_comp]) + " " + str(strains_all_std[k_comp]) for k_comp in xrange(6)]))
+        if (n_sector_ids > 1):
+            for sector_id in xrange(n_sector_ids):
+                strains_per_sector_avg = numpy.mean(strains_per_sector[sector_id], 0)
+                strains_per_sector_std = numpy.std(strains_per_sector[sector_id], 0)
+                strain_file.write("".join([" " + str(strains_per_sector_avg[k_comp]) + " " + str(strains_per_sector_std[k_comp]) for k_comp in xrange(6)]))
+        strain_file.write("\n")
+
+        if (write_strain_vs_radius):
+            assert (ref_mesh.GetCellData().HasArray("r"))
+            farray_r = ref_mesh.GetCellData().GetArray("r")
+            for k_cell in xrange(n_cells):
+                strain_vs_radius_file.write(" ".join([str(val) for val in [k_frame, farray_r.GetTuple1(k_cell)]+list(farray_strain.GetTuple(k_cell))]) + "\n")
+            strain_vs_radius_file.write("\n")
+            strain_vs_radius_file.write("\n")
+
+    strain_file.close()
+    strain_vs_radius_file.close()

fedic.py

+#coding=utf8
+
+########################################################################
+###                                                                  ###
+### Created by Martin Genet, 2016                                    ###
+###                                                                  ###
+### École Polytechnique, Palaiseau, France                           ###
+###                                                                  ###
+########################################################################
+
+import dolfin
+import glob
+import numpy
+import os
+
+import myFEniCSPythonLibrary as myFEniCS
+import myVTKPythonLibrary    as myVTK
+
+########################################################################
+
+def fedic(
+        working_folder,
+        working_basename,
+        images_folder,
+        images_basename,
+        mesh_folder,
+        mesh_basename,
+        images_dimension=3,
+        k_ref=0,
+        penalty=0.9,
+        use_I0_tangent=0,
+        relax_type="const",
+        relax_init=0.9,
+        tol_disp=1e-3,
+        n_iter_max=100,
+        print_iterations=0,
+        continue_after_fail=0,
+        verbose=1):
+
+    if not os.path.exists(working_folder):
+        os.mkdir(working_folder)
+
+    print "Loading mesh…"
+    mesh = dolfin.Mesh(mesh_folder+"/"+mesh_basename+".xml")
+
+    fs = dolfin.VectorFunctionSpace(
+        mesh=mesh,
+        family="Lagrange",
+        degree=1)
+    dX = dolfin.dx(mesh)
+
+    print "Computing quadrature degree…"
+    degree = myFEniCS.compute_quadrature_degree(
+        image_filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
+        dX=dX,
+        image_dimension=images_dimension)
+    print "degree = " + str(degree)
+
+    fe = dolfin.FiniteElement(
+        family="Quadrature",
+        degree=degree)
+
+    print "Loading reference image…"
+    if (images_dimension == 2):
+        I0 = myFEniCS.ExprIm2(
+            filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
+            element=fe)
+        if (use_I0_tangent):
+            DXI0 = myFEniCS.ExprGradXIm2(
+                filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
+                element=fe)
+            DYI0 = myFEniCS.ExprGradYIm2(
+                filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
+                element=fe)
+    elif (images_dimension == 3):
+        I0 = myFEniCS.ExprIm3(
+            filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
+            element=fe)
+        if (use_I0_tangent):
+            DXI0 = myFEniCS.ExprGradXIm3(
+                filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
+                element=fe)
+            DYI0 = myFEniCS.ExprGradYIm3(
+                filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
+                element=fe)
+            DZI0 = myFEniCS.ExprGradZIm3(
+                filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
+                element=fe)
+    else:
+        assert (0), "images_dimension must be 2 or 3. Aborting."
+    I0_norm = dolfin.assemble(I0**2 * dX)**(1./2)
+
+    print "Defining functions…"
+    penalty = dolfin.Constant(penalty)
+    U = dolfin.Function(
+        fs,
+        name="displacement")
+    dU = dolfin.Function(
+        fs,
+        name="ddisplacement")
+    ddU = dolfin.TrialFunction(fs)
+    V = dolfin.TestFunction(fs)
+
+    print "Printing initial solution…"
+    file_pvd = dolfin.File(working_folder+"/"+working_basename+"_.pvd")
+    for vtu in glob.glob(working_folder+"/"+working_basename+"_*.vtu"):
+        os.remove(vtu)
+    file_pvd << (U, float(k_ref))
+
+    print "Defining variational forms…"
+    if (images_dimension == 2):
+        I1 = myFEniCS.ExprDefIm2(
+            U=U,
+            element=fe)
+        DXI1 = myFEniCS.ExprGradXDefIm2(
+            U=U,
+            element=fe)
+        DYI1 = myFEniCS.ExprGradYDefIm2(
+            U=U,
+            element=fe)
+        if (use_I0_tangent):
+            a =     penalty  * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI0, DYI0)), ddU),
+                                            dolfin.dot(dolfin.as_vector((DXI0, DYI0)),   V))*dX\
+              + (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
+                                            dolfin.grad(  V))*dX
+        else:
+            a =     penalty  * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI1, DYI1)), ddU),
+                                            dolfin.dot(dolfin.as_vector((DXI1, DYI1)),   V))*dX\
+              + (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
+                                            dolfin.grad(  V))*dX
+        b =     penalty  * dolfin.inner(I0-I1,
+                                        dolfin.dot(dolfin.as_vector((DXI1, DYI1)), V))*dX\
+          - (1.-penalty) * dolfin.inner(dolfin.grad(U),
+                                        dolfin.grad(V))*dX
+    elif (images_dimension == 3):
+        I1 = myFEniCS.ExprDefIm3(
+            U=U,
+            element=fe)
+        DXI1 = myFEniCS.ExprGradXDefIm3(
+            U=U,
+            element=fe)
+        DYI1 = myFEniCS.ExprGradYDefIm3(
+            U=U,
+            element=fe)
+        DZI1 = myFEniCS.ExprGradZDefIm3(
+            U=U,
+            element=fe)
+        if (use_I0_tangent):
+            a =     penalty  * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI0, DYI0, DZI0)), ddU),
+                                            dolfin.dot(dolfin.as_vector((DXI0, DYI0, DZI0)),   V))*dX\
+              + (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
+                                            dolfin.grad(  V))*dX
+        else:
+            a =     penalty  * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI1, DYI1, DZI1)), ddU),
+                                            dolfin.dot(dolfin.as_vector((DXI1, DYI1, DZI1)),   V))*dX\
+              + (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
+                                            dolfin.grad(  V))*dX
+        b =     penalty  * dolfin.inner(I0-I1,
+                                        dolfin.dot(dolfin.as_vector((DXI1, DYI1, DZI1)), V))*dX\
+          - (1.-penalty) * dolfin.inner(dolfin.grad(U),
+                                        dolfin.grad(V))*dX
+    else:
+        assert (0), "images_dimension must be 2 or 3. Aborting."
+
+    # linear system
+    A = None
+    B = None
+    if (use_I0_tangent):
+        A = dolfin.assemble(a, tensor=A)
+
+    print "Checking number of frames…"
+    n_frames = len(glob.glob(images_folder+"/"+images_basename+"_*.vti"))
+    print "n_frames = " + str(n_frames)
+
+    assert (abs(k_ref) < n_frames)
+    k_ref = k_ref%n_frames
+
+    n_iter_tot = 0
+    for forward_or_backward in ["forward", "backward"]: # not sure if this still works…
+        print "forward_or_backward = " + forward_or_backward
+
+        if (forward_or_backward == "forward"):
+            k_next_frames = range(k_ref+1, n_frames, +1)
+        elif (forward_or_backward == "backward"):
+            k_next_frames = range(k_ref-1, -1, -1)
+        print "k_next_frames = " + str(k_next_frames)
+
+        for k_next_frame in k_next_frames:
+            print "k_next_frame = " + str(k_next_frame)
+
+            if (print_iterations):
+                if (os.path.exists(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+".pdf")):
+                    os.remove(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+".pdf")
+                file_dat_iter = open(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+".dat", "w")
+
+                file_pvd_iter = dolfin.File(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+"_.pvd")
+                for vtu in glob.glob(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+"_*.vtu"):
+                    os.remove(vtu)
+                file_pvd_iter << (U, 0.)
+
+            print "Loading image and image gradient…"
+            if (images_dimension == 2):
+                I1.init_image(  filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
+                DXI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
+                DYI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
+            elif (images_dimension == 3):
+                I1.init_image(  filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
+                DXI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
+                DYI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
+                DZI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
+            else:
+                assert (0), "images_dimension must be 2 or 3. Aborting."
+
+            print "Running registration…"
+            k_iter = 0
+            relax = relax_init
+            while (True):
+                print "k_iter = " + str(k_iter)
+                n_iter_tot += 1
+
+                # linear system
+                if not (use_I0_tangent):
+                    A = dolfin.assemble(a, tensor=A)
+                #print "A = " + str(A.array())
+                if (k_iter > 0):
+                    if (k_iter == 1):
+                        B_old = B.copy()
+                    elif (k_iter > 1):
+                        B_old[:] = B[:]
+                    #print "B_old = " + str(B_old[0])
+                B = dolfin.assemble(b, tensor=B)
+                #print "B = " + str(B[0])
+                dolfin.solve(A, dU.vector(), B)
+                #print "dU = " + str(dU.vector().array())
+
+                # relaxation
+                if (relax_type == "aitken") and (k_iter > 0):
+                    if (k_iter == 1):
+                        dB = B - B_old
+                    elif (k_iter > 1):
+                        dB[:] = B[:] - B_old[:]
+                    relax *= (-1.) * B_old.inner(dB) / dB.inner(dB)
+                print "relax = " + str(relax)
+
+                # solution update
+                U.vector()[:] += relax * dU.vector()[:]
+                #U.vector().axpy(relax, dU.vector())
+
+                if (print_iterations):
+                    #print "U = " + str(U.vector().array())
+                    file_pvd_iter << (U, float(k_iter+1))
+
+                # displacement error
+                dU_max = numpy.linalg.norm(dU.vector().array(), ord=numpy.Inf)
+                U_max = numpy.linalg.norm(U.vector().array(), ord=numpy.Inf)
+                err_disp_abs = max(dU_max, U_max)
+                U_norm = numpy.linalg.norm(U.vector().array())
+                err_disp_rel = dU_max/U_norm
+                print "err_disp_abs = " + str(err_disp_abs)
+                print "err_disp_rel = " + str(err_disp_rel)
+
+                # image error
+                I1I0_norm = dolfin.assemble((I1-I0)**2 * dX)**(1./2)
+                err_im = I1I0_norm/I0_norm
+                print "err_im = " + str(err_im)
+
+                if (print_iterations):
+                    file_dat_iter.write(" ".join([str(val) for val in [k_iter, err_disp_abs, err_disp_rel, err_im]])+"\n")
+
+                # exit test
+                if (err_disp_rel < tol_disp) or (err_disp_abs < tol_disp):
+                    print "Nonlinear solver converged…"
+                    success = True
+                    break
+
+                if (k_iter >= n_iter_max-1):
+                    print "Warning! Nonlinear solver failed to converge…"
+                    success = False
+                    break
+
+                # increment counter
+                k_iter += 1
+
+            if (print_iterations):
+                os.remove(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+"_.pvd")
+                file_dat_iter.close()
+                os.system("gnuplot -e \"set terminal pdf; set output '"+working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+".pdf'; set logscale y; plot '"+working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+".dat' using 1:3\"")
+
+            if not (success) and not (continue_after_fail):
+                break
+
+            print "Printing solution…"
+            file_pvd << (U, float(k_next_frame))
+
+        if not (success) and not (continue_after_fail):
+            break
+
+    print "n_iter_tot = " + str(n_iter_tot)
+
+    os.remove(working_folder+"/"+working_basename+"_.pvd")
+
+    return success

plot_displacement_error.py

+#coding=utf8
+
+########################################################################
+###                                                                  ###
+### Created by Martin Genet, 2016                                    ###
+###                                                                  ###
+### École Polytechnique, Palaiseau, France                           ###
+###                                                                  ###
+########################################################################
+
+import math
+import numpy
+import os
+import sys
+
+########################################################################
+
+def plot_displacement_error(
+        working_folder,
+        working_basenames,
+        suffix="",
+        verbose=1):
+
+    n_frames = len(open(working_folder+"/"+working_basenames[0]+"-error.dat").readlines())-1
+    #print "n_frames = " + str(n_frames)
+
+    plotfile = open("plot_displacement_error"+("-"+suffix)*(suffix!="")+".plt", "w")
+    plotfile.write('''\
+set terminal pdf enhanced size 5,3
+
+set output "plot_displacement_error'''+('''-'''+suffix)*(suffix!="")+'''.pdf"
+
+load "Set1.plt"
+set linestyle 1 pointtype 1
+set linestyle 2 pointtype 1
+set linestyle 3 pointtype 1
+set linestyle 4 pointtype 1
+set linestyle 5 pointtype 1
+set linestyle 6 pointtype 1
+set linestyle 7 pointtype 1
+set linestyle 8 pointtype 1
+set linestyle 9 pointtype 1
+
+set key box textcolor variable width +0
+
+set grid
+
+''')
+
+    plotfile.write('''\
+set xlabel "frame number"
+set xrange [0:'''+str(n_frames-1)+''']
+
+set ylabel "displacement error (%)"
+set yrange [0:*]
+
+''')
+
+    for k_basename in range(len(working_basenames)):
+        working_basename = working_basenames[k_basename]
+        plotfile.write('''\
+'''+(k_basename==0)*('''plot ''')+(k_basename>0)*('''     ''')+'''"'''+working_folder+'''/'''+working_basename+'''-error.dat" using ($1):(100*$4) with lines linestyle '''+str(k_basename+1)+''' linewidth 5 title "'''+working_basename+'''"'''+(k_basename<len(working_basenames)-1)*(''',\\
+''')+(k_basename==len(working_basenames)-1)*('''
+
+'''))
+
+    plotfile.close()
+
+    os.system("gnuplot plot_displacement_error"+("-"+suffix)*(suffix!="")+".plt")

plot_strains.py

+#coding=utf8
+
+########################################################################
+###                                                                  ###
+### Created by Martin Genet, 2016                                    ###
+###                                                                  ###
+### École Polytechnique, Palaiseau, France                           ###
+###                                                                  ###
+########################################################################
+
+import math
+import numpy
+import os
+import sys
+
+########################################################################
+
+def plot_strains(
+        working_folder,
+        working_basenames,
+        ref_folder=None,
+        ref_basename=None,
+        components="all",
+        suffix="",
+        verbose=1):
+
+    assert (components in ("all", "circ-long", "rad-circ"))
+
+    if (ref_folder is not None) and (ref_basename is not None):
+        lines = open(ref_folder+"/"+ref_basename+"-strains.dat").readlines()
+    else:
+        lines = open(working_folder+"/"+working_basenames[0]+"-strains.dat").readlines()
+    n_frames = len(lines)-1
+    n_sectors = (len(lines[-1].split(" "))-1)/12
+    #print "n_frames = " + str(n_frames)
+    #print "n_sectors = " + str(n_sectors)
+
+    plotfile = open("plot_strains"+("-"+suffix)*(suffix!="")+".plt", "w")
+    plotfile.write('''\
+set terminal pdf enhanced size 15,'''+('''6''')*(components=="all")+('''3''')*(components in ("circ-long", "rad-circ"))+'''
+
+set output "plot_strains'''+('''-'''+suffix)*(suffix!="")+'''.pdf"
+
+load "Set1.plt"
+set linestyle 1 pointtype 1
+set linestyle 2 pointtype 1
+set linestyle 3 pointtype 1
+set linestyle 4 pointtype 1
+set linestyle 5 pointtype 1
+set linestyle 6 pointtype 1
+set linestyle 7 pointtype 1
+set linestyle 8 pointtype 1
+set linestyle 9 pointtype 1
+
+set key box textcolor variable width +0
+
+set grid
+
+''')
+
+    for k_sector in range(n_sectors):
+        plotfile.write('''\
+set multiplot layout '''+('''2''')*(components=="all")+('''1''')*(components in ("circ-long", "rad-circ"))+''',3
+
+''')
+        if (components == "all"):
+            comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
+        elif (components == "circ-long"):
+            comp_names = ["circumferential","longitudinal","circumferential-longitudinal"]
+        elif (components == "rad-circ"):
+            comp_names = ["radial", "circumferential", "radial-circumferential"]
+        else:
+            assert (0), "components must be \"all\", \"circ-long\" or \"rad-circ\". Aborting."
+        for comp_name in comp_names:
+            if   (comp_name == "radial"                      ): k_comp = 0
+            elif (comp_name == "circumferential"             ): k_comp = 1
+            elif (comp_name == "longitudinal"                ): k_comp = 2
+            elif (comp_name == "radial-circumferential"      ): k_comp = 3
+            elif (comp_name == "radial-longitudinal"         ): k_comp = 4
+            elif (comp_name == "circumferential-longitudinal"): k_comp = 5
+            plotfile.write('''\
+set xrange [0:'''+str(n_frames)+''']
+
+set ylabel "'''+comp_name+''' strain (%)"
+set yrange [-30:30]
+
+plot 0 linecolor rgb "black" notitle,\\
+''')
+            if (ref_folder is not None) and (ref_basename is not None):
+                plotfile.write('''\
+     "'''+ref_folder+'''/'''+ref_basename+'''-strains.dat" using ($1):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linecolor "black" linewidth 5 notitle,\
+     "'''+ref_folder+'''/'''+ref_basename+'''-strains.dat" using ($1):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linecolor "black" pointtype 1 linewidth 1 notitle'''+(len(working_basenames)>0)*(''',\\
+''')+(len(working_basenames)==0)*('''
+
+'''))
+            for k_basename in range(len(working_basenames)):
+                working_basename = working_basenames[k_basename]
+                plotfile.write('''\
+     "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ($1+'''+str(k_basename)+'''./100):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linestyle '''+str(k_basename+1)+''' linewidth 5 title "'''+working_basename+'''",\
+     "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ($1+'''+str(k_basename)+'''./100):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linestyle '''+str(k_basename+1)+''' linewidth 1 notitle'''+(k_basename<len(working_basenames)-1)*(''',\\
+''')+(k_basename==len(working_basenames)-1)*('''
+
+'''))
+
+    plotfile.close()
+
+    os.system("gnuplot plot_strains"+("-"+suffix)*(suffix!="")+".plt")

plot_strains_vs_radius.py

+#coding=utf8
+
+########################################################################
+###                                                                  ###
+### Created by Martin Genet, 2016                                    ###
+###                                                                  ###
+### École Polytechnique, Palaiseau, France                           ###
+###                                                                  ###
+########################################################################
+
+import math
+import numpy
+import os
+import sys
+
+########################################################################
+
+def plot_strains_vs_radius(
+        working_folder,
+        working_basenames,
+        index=10,
+        ref_folder=None,
+        ref_basename=None,
+        components="all",
+        suffix="",
+        verbose=1):
+
+    assert (components in ("all", "circ-long", "rad-circ"))
+
+    plotfile = open("plot_strains_vs_radius"+("-"+suffix)*(suffix!="")+".plt", "w")
+    plotfile.write('''\
+set terminal pdf enhanced size 15,'''+('''6''')*(components=="all")+('''3''')*(components in ("circ-long", "rad-circ"))+'''
+
+set output "plot_strains_vs_radius'''+('''-'''+suffix)*(suffix!="")+'''.pdf"
+
+load "Set1.plt"
+set linestyle 1 pointtype 1
+set linestyle 2 pointtype 1
+set linestyle 3 pointtype 1
+set linestyle 4 pointtype 1
+set linestyle 5 pointtype 1
+set linestyle 6 pointtype 1
+set linestyle 7 pointtype 1
+set linestyle 8 pointtype 1
+set linestyle 9 pointtype 1
+
+set key box textcolor variable width +0
+
+set grid
+
+set multiplot layout '''+('''2''')*(components=="all")+('''1''')*(components in ("circ-long", "rad-circ"))+''',3
+
+''')
+    if (components == "all"):
+        comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
+    elif (components == "circ-long"):
+        comp_names = ["circumferential","longitudinal","circumferential-longitudinal"]
+    elif (components == "rad-circ"):
+        comp_names = ["radial", "circumferential", "radial-circumferential"]
+    else:
+        assert (0), "components must be \"all\", \"circ-long\" or \"rad-circ\". Aborting."
+    for comp_name in comp_names:
+        if   (comp_name == "radial"                      ): k_comp = 0
+        elif (comp_name == "circumferential"             ): k_comp = 1
+        elif (comp_name == "longitudinal"                ): k_comp = 2
+        elif (comp_name == "radial-circumferential"      ): k_comp = 3
+        elif (comp_name == "radial-longitudinal"         ): k_comp = 4
+        elif (comp_name == "circumferential-longitudinal"): k_comp = 5
+        plotfile.write('''\
+
+set xlabel "r (mm)"
+set xrange [0.2:0.4]
+
+set ylabel "'''+comp_name+''' strain (%)"
+set yrange [-40:40]
+
+plot 0 linecolor rgb "black" notitle,\\
+''')
+        if (ref_folder is not None) and (ref_basename is not None):
+            plotfile.write('''\
+     "'''+ref_folder+'''/'''+ref_basename+'''-strains_vs_radius.dat" using ($2):(100*$'''+str(3+k_comp)+''') index '''+str(index)+''' with points linecolor "black" pointtype 1 pointsize 1 linewidth 3 notitle'''+(len(working_basenames)>0)*(''',\\
+''')+(len(working_basenames)==0)*('''
+
+'''))
+        for k_basename in range(len(working_basenames)):
+            working_basename = working_basenames[k_basename]
+            plotfile.write('''\
+     "'''+working_folder+'''/'''+working_basename+'''-strains_vs_radius.dat" using ($2):(100*$'''+str(3+k_comp)+''') index '''+str(index)+''' with points linestyle '''+str(k_basename+1)+''' pointtype 1 pointsize 1 linewidth 3 title "'''+working_basename+'''"'''+(k_basename<len(working_basenames)-1)*(''',\\
+''')+(k_basename==len(working_basenames)-1)*('''
+
+'''))
+
+    plotfile.close()
+
+    os.system("gnuplot plot_strains_vs_radius"+("-"+suffix)*(suffix!="")+".plt")