Merge branch 'master' into GIMIC

__init__.py

@@ -16,6 +16,7 @@ from Energy_Regularization import *
 from compute_warped_images import *
 from plot_strains import *
 from image_expressions_cpp import *
+from plot_twist_vs_height import *
 from Energy import *
 from compute_warped_mesh import *
 from Energy_WarpedImage import *

compute_strains.py

@@ -10,10 +10,12 @@
 
 import dolfin
 import glob
+import math
 import numpy
 
 import myPythonLibrary as mypy
 import myVTKPythonLibrary as myvtk
+import vtkpython_cbl as cbl
 
 import dolfin_dic as ddic
 
@@ -23,50 +25,78 @@ def compute_strains(
         working_folder,
         working_basename,
         working_ext="vtu",
-        ref_frame=None,
+        ref_frame=None,                             # MG20190612: Reference configuration.
         disp_array_name="displacement",
         defo_grad_array_name="DeformationGradient",
         strain_array_name="Strain",
-        ref_mesh_folder=None,
+        ref_mesh_folder=None,                       # MG20190612: Mesh with sectors/parts/etc.
         ref_mesh_basename=None,
         ref_mesh_ext="vtk",
         CYL_or_PPS="PPS",
+        remove_boundary_layer=1,
+        in_place=1,
         write_strains=1,
         temporal_offset=None,
         temporal_resolution=None,
         plot_strains=1,
         plot_regional_strains=0,
         write_strains_vs_radius=0,
+        plot_strains_vs_radius=0,
         write_binned_strains_vs_radius=0,
+        plot_binned_strains_vs_radius=0,
+        write_twist_vs_height=0,
+        plot_twist_vs_height=0,
         verbose=1):
 
-    if (ref_mesh_folder is not None) and (ref_mesh_basename is not None):
+    if (write_strains_vs_radius):
+        assert (remove_boundary_layer),\
+            "write_strains_vs_radius only works after removing the boundary layer. Aborting"
+    if (write_binned_strains_vs_radius):
+        assert (remove_boundary_layer),\
+            "write_binned_strains_vs_radius only works after removing the boundary layer. Aborting"
+    if (write_twist_vs_height):
+        assert (remove_boundary_layer),\
+            "write_twist_vs_height only works after removing the boundary layer. Aborting"
+
+    if  (ref_mesh_folder   is not None)\
+    and (ref_mesh_basename is not None):
+        ref_mesh_filename = ref_mesh_folder+"/"+ref_mesh_basename+"."+ref_mesh_ext
         ref_mesh = myvtk.readUGrid(
-            filename=ref_mesh_folder+"/"+ref_mesh_basename+"."+ref_mesh_ext,
+            filename=ref_mesh_filename,
             verbose=verbose)
+        ref_mesh_n_points = ref_mesh.GetNumberOfPoints()
         ref_mesh_n_cells = ref_mesh.GetNumberOfCells()
+        if (verbose): print "ref_mesh_n_points = " + str(ref_mesh_n_points)
         if (verbose): print "ref_mesh_n_cells = " + str(ref_mesh_n_cells)
 
+        if (ref_mesh.GetCellData().HasArray("part_id")):
+            iarray_ref_part_id = ref_mesh.GetCellData().GetArray("part_id")
+            n_part_ids = 0
+            for k_cell in xrange(ref_mesh_n_cells):
+                part_id = int(iarray_ref_part_id.GetTuple1(k_cell))
+                if (part_id > n_part_ids-1):
+                    n_part_ids = part_id+1
+            if (verbose): print "n_part_ids = " + str(n_part_ids)
+        else:
+            iarray_ref_part_id = None
+            n_part_ids = 0
+
         if (ref_mesh.GetCellData().HasArray("sector_id")):
-            iarray_sector_id = ref_mesh.GetCellData().GetArray("sector_id")
+            iarray_ref_sector_id = ref_mesh.GetCellData().GetArray("sector_id")
             n_sector_ids = 0
             for k_cell in xrange(ref_mesh_n_cells):
-                sector_id = int(iarray_sector_id.GetTuple1(k_cell))
+                sector_id = int(iarray_ref_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:
-            iarray_sector_id = None
+            iarray_ref_sector_id = None
             n_sector_ids = 0
 
-        if (ref_mesh.GetCellData().HasArray("part_id")):
-            part_id_array = ref_mesh.GetCellData().GetArray("part_id")
-        else:
-            part_id_array = None
-
     else:
         ref_mesh = None
+        n_part_ids = 0
         n_sector_ids = 0
 
     working_filenames = glob.glob(working_folder+"/"+working_basename+"_[0-9]*."+working_ext)
@@ -83,23 +113,27 @@ def compute_strains(
         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_strains_vs_radius):
-        strain_vs_radius_file = open(working_folder+"/"+working_basename+"-strains_vs_radius.dat", "w")
-        strain_vs_radius_file.write("#t rr Err Ecc Ell Erc Erl Ecl\n")
+        strains_vs_radius_file = open(working_folder+"/"+working_basename+"-strains_vs_radius.dat", "w")
+        strains_vs_radius_file.write("#t rr Err Ecc Ell Erc Erl Ecl\n")
 
     if (write_binned_strains_vs_radius):
-        binned_strain_vs_radius_file = open(working_folder+"/"+working_basename+"-binned_strains_vs_radius.dat", "w")
-        binned_strain_vs_radius_file.write("#t rr Err Ecc Ell Erc Erl Ecl\n")
+        binned_strains_vs_radius_file = open(working_folder+"/"+working_basename+"-binned_strains_vs_radius.dat", "w")
+        binned_strains_vs_radius_file.write("#t rr Err Ecc Ell Erc Erl Ecl\n")
+
+    if (write_twist_vs_height):
+        twist_vs_height_file = open(working_folder+"/"+working_basename+"-twist_vs_height.dat", "w")
+        twist_vs_height_file.write("#t z beta\n")
 
     if (ref_frame is not None):
-        ref_mesh_filename = working_folder+"/"+working_basename+"_"+str(ref_frame).zfill(working_zfill)+"."+working_ext
-        ref_mesh = myvtk.readUGrid(
-            filename=ref_mesh_filename,
+        mesh0_filename = working_folder+"/"+working_basename+"_"+str(ref_frame).zfill(working_zfill)+"."+working_ext
+        mesh0 = myvtk.readUGrid(
+            filename=mesh0_filename,
             verbose=verbose)
         myvtk.addDeformationGradients(
-            mesh=ref_mesh,
+            mesh=mesh0,
             disp_array_name=disp_array_name,
             verbose=verbose)
-        farray_F0 = ref_mesh.GetCellData().GetArray(defo_grad_array_name)
+        farray_F0 = mesh0.GetCellData().GetArray(defo_grad_array_name)
 
     for k_frame in xrange(n_frames):
         print "k_frame = "+str(k_frame)
@@ -108,13 +142,25 @@ def compute_strains(
         mesh = myvtk.readUGrid(
             filename=mesh_filename,
             verbose=verbose)
+        n_points = mesh.GetNumberOfPoints()
         n_cells = mesh.GetNumberOfCells()
         if (ref_mesh is not None):
-            assert (n_cells == ref_mesh_n_cells), "ref_mesh_n_cells ("+str(ref_mesh_n_cells)+") ≠ n_cells ("+str(n_cells)+"). Aborting."
-            if (part_id_array is not None):
-                mesh.GetCellData().AddArray(part_id_array)
-            if (iarray_sector_id is not None):
-                mesh.GetCellData().AddArray(iarray_sector_id)
+            assert (n_points == ref_mesh_n_points),\
+                "ref_mesh_n_points ("+str(ref_mesh_n_points)+") ≠ n_points ("+str(n_points)+"). Aborting."
+            assert (n_cells == ref_mesh_n_cells),\
+                "ref_mesh_n_cells ("+str(ref_mesh_n_cells)+") ≠ n_cells ("+str(n_cells)+"). Aborting."
+            if (iarray_ref_part_id is not None):
+                mesh.GetCellData().AddArray(iarray_ref_part_id)
+            if (iarray_ref_sector_id is not None):
+                mesh.GetCellData().AddArray(iarray_ref_sector_id)
+            if (write_strains_vs_radius       )\
+            or (write_binned_strains_vs_radius):
+                assert (ref_mesh.GetCellData().HasArray("rr"))
+                mesh.GetCellData().AddArray(ref_mesh.GetCellData().GetArray("rr"))
+            if (write_twist_vs_height):
+                assert (ref_mesh.GetPointData().HasArray("ll"))
+                mesh.GetPointData().AddArray(ref_mesh.GetPointData().GetArray("r"))
+                mesh.GetPointData().AddArray(ref_mesh.GetPointData().GetArray("ll"))
         myvtk.addDeformationGradients(
             mesh=mesh,
             disp_array_name=disp_array_name,
@@ -133,38 +179,62 @@ def compute_strains(
             strain_array_name=strain_array_name,
             mesh_w_local_basis=ref_mesh,
             verbose=verbose)
+        if  (ref_mesh           is not None)\
+        and (iarray_ref_part_id is not None)\
+        and (remove_boundary_layer         ):
+            mesh = myvtk.getThresholdedUGrid(
+                ugrid=mesh,
+                field_support="cells",
+                field_name="part_id",
+                threshold_value=0.5,
+                threshold_by_upper_or_lower="lower")
+            n_points = mesh.GetNumberOfPoints()
+            n_cells = mesh.GetNumberOfCells()
+            n_part_ids = 0
+            if (iarray_ref_sector_id is not None):
+                iarray_sector_id = mesh.GetCellData().GetArray("sector_id")
+        else:
+            iarray_sector_id = iarray_ref_sector_id
+        mesh_filename = working_folder+"/"+working_basename+("-wStrains")*(not in_place)+"_"+str(k_frame).zfill(working_zfill)+"."+working_ext
         myvtk.writeUGrid(
             ugrid=mesh,
-            filename=working_folder+"/"+working_basename+"_"+str(k_frame).zfill(working_zfill)+"."+working_ext,
+            filename=mesh_filename,
             verbose=verbose)
 
-        if (write_strains) or (write_strains_vs_radius) or (write_binned_strains_vs_radius):
-            if (ref_mesh is not None) and (mesh.GetCellData().HasArray(strain_array_name+"_"+CYL_or_PPS)):
+        if (write_strains                 )\
+        or (write_strains_vs_radius       )\
+        or (write_binned_strains_vs_radius):
+            if  (ref_mesh is not None)\
+            and (mesh.GetCellData().HasArray(strain_array_name+"_"+CYL_or_PPS)):
                 farray_strain = mesh.GetCellData().GetArray(strain_array_name+"_"+CYL_or_PPS)
             else:
                 farray_strain = mesh.GetCellData().GetArray(strain_array_name)
 
         if (write_strains):
-            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))
+            if (n_sector_ids in (0,1)):
+                if (n_part_ids == 0):
+                    strains_all = [farray_strain.GetTuple(k_cell) for k_cell in xrange(n_cells)]
+                else:
+                    strains_all = [farray_strain.GetTuple(k_cell) for k_cell in xrange(n_cells) if (iarray_ref_part_id.GetTuple1(k_cell) > 0)]
             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))
+                if (n_part_ids == 0):
+                    for k_cell in xrange(n_cells):
+                        strains_all.append(farray_strain.GetTuple(k_cell))
+                        sector_id = int(iarray_sector_id.GetTuple1(k_cell))
+                        strains_per_sector[sector_id].append(farray_strain.GetTuple(k_cell))
+                else:
+                    for k_cell in xrange(n_cells):
+                        part_id = int(iarray_ref_part_id.GetTuple1(k_cell))
+                        if (part_id > 0): continue
+                        strains_all.append(farray_strain.GetTuple(k_cell))
+                        sector_id = int(iarray_sector_id.GetTuple1(k_cell))
+                        if (sector_id < 0): continue
+                        strains_per_sector[sector_id].append(farray_strain.GetTuple(k_cell))
 
-            if (temporal_offset is not None) and (temporal_resolution is not None):
+            if  (temporal_offset     is not None)\
+            and (temporal_resolution is not None):
                 strain_file.write(str(temporal_offset + k_frame*temporal_resolution))
             else:
                 strain_file.write(str(k_frame))
@@ -179,16 +249,14 @@ def compute_strains(
             strain_file.write("\n")
 
         if (write_strains_vs_radius):
-            assert (ref_mesh.GetCellData().HasArray("rr"))
-            farray_rr = ref_mesh.GetCellData().GetArray("rr")
+            farray_rr = mesh.GetCellData().GetArray("rr")
             for k_cell in xrange(n_cells):
-                strain_vs_radius_file.write(" ".join([str(val) for val in [k_frame, farray_rr.GetTuple1(k_cell)]+list(farray_strain.GetTuple(k_cell))]) + "\n")
-            strain_vs_radius_file.write("\n")
-            strain_vs_radius_file.write("\n")
+                strains_vs_radius_file.write(" ".join([str(val) for val in [k_frame, farray_rr.GetTuple1(k_cell)]+list(farray_strain.GetTuple(k_cell))]) + "\n")
+            strains_vs_radius_file.write("\n")
+            strains_vs_radius_file.write("\n")
 
         if (write_binned_strains_vs_radius):
-            assert (ref_mesh.GetCellData().HasArray("rr"))
-            farray_rr = ref_mesh.GetCellData().GetArray("rr")
+            farray_rr = mesh.GetCellData().GetArray("rr")
             n_r = 10
             binned_strains = [[] for k_r in xrange(n_r)]
             for k_cell in xrange(n_cells):
@@ -203,9 +271,63 @@ def compute_strains(
             #print binned_strains_avg
             #print binned_strains_std
             for k_r in xrange(n_r):
-                binned_strain_vs_radius_file.write(" ".join([str(val) for val in [k_frame, (k_r+0.5)/n_r]+[val for k_comp in xrange(6) for val in [binned_strains_avg[k_r][k_comp], binned_strains_std[k_r][k_comp]]]]) + "\n")
-            binned_strain_vs_radius_file.write("\n")
-            binned_strain_vs_radius_file.write("\n")
+                binned_strains_vs_radius_file.write(" ".join([str(val) for val in [k_frame, (k_r+0.5)/n_r]+[val for k_comp in xrange(6) for val in [binned_strains_avg[k_r][k_comp], binned_strains_std[k_r][k_comp]]]]) + "\n")
+            binned_strains_vs_radius_file.write("\n")
+            binned_strains_vs_radius_file.write("\n")
+
+        if (write_twist_vs_height):
+            points_AB = cbl.getABPointsFromBoundsAndCenter(
+                mesh=mesh,
+                AB=[0,0,1],
+                verbose=verbose)
+            C = points_AB.GetPoint(0)
+
+            farray_r  = mesh.GetPointData().GetArray("r")
+            farray_ll = mesh.GetPointData().GetArray("ll")
+            farray_U  = mesh.GetPointData().GetArray(disp_array_name)
+            farray_Theta = myvtk.createFloatArray(
+                name="Theta",
+                n_components=1,
+                n_tuples=n_points)
+            farray_theta = myvtk.createFloatArray(
+                name="theta",
+                n_components=1,
+                n_tuples=n_points)
+            farray_beta = myvtk.createFloatArray(
+                name="beta",
+                n_components=1,
+                n_tuples=n_points)
+            mesh.GetPointData().AddArray(farray_Theta)
+            mesh.GetPointData().AddArray(farray_theta)
+            mesh.GetPointData().AddArray(farray_beta)
+            X = numpy.empty(3)
+            U = numpy.empty(3)
+            x = numpy.empty(3)
+            for k_point in xrange(n_points):
+                r  = farray_r.GetTuple1(k_point)
+                ll = farray_ll.GetTuple1(k_point)
+                mesh.GetPoint(k_point, X)
+                X -= C
+                Theta = math.degrees(math.atan2(X[1], X[0]))
+                farray_U.GetTuple(k_point, U)
+                x[:] = X[:] + U[:]
+                theta = math.degrees(math.atan2(x[1], x[0]))
+                beta = theta - Theta
+                if (beta > +180.): beta -= 360.
+                if (beta < -180.): beta += 360.
+                farray_Theta.SetTuple1(k_point, Theta)
+                farray_theta.SetTuple1(k_point, theta)
+                farray_beta.SetTuple1(k_point, beta)
+                if (r < 15.): continue
+                # if (ll < 1./3): continue
+                twist_vs_height_file.write(" ".join([str(val) for val in [k_frame, ll, beta]]) + "\n")
+            mesh_filename = working_folder+"/"+working_basename+("-wStrains")*(not in_place)+"_"+str(k_frame).zfill(working_zfill)+"."+working_ext
+            myvtk.writeUGrid(
+                ugrid=mesh,
+                filename=mesh_filename,
+                verbose=verbose)
+            twist_vs_height_file.write("\n")
+            twist_vs_height_file.write("\n")
 
     if (write_strains):
         strain_file.close()
@@ -225,7 +347,31 @@ def compute_strains(
                 verbose=verbose)
 
     if (write_strains_vs_radius):
-        strain_vs_radius_file.close()
+        strains_vs_radius_file.close()
+
+        if (plot_strains_vs_radius):
+            ddic.plot_strains_vs_radius(
+                working_folder=working_folder,
+                working_basenames=[working_basename],
+                suffix=None,
+                verbose=verbose)
 
     if (write_binned_strains_vs_radius):
-        binned_strain_vs_radius_file.close()
+        binned_strains_vs_radius_file.close()
+
+        if (plot_binned_strains_vs_radius):
+            ddic.plot_binned_strains_vs_radius(
+                working_folder=working_folder,
+                working_basenames=[working_basename],
+                suffix=None,
+                verbose=verbose)
+
+    if (write_twist_vs_height):
+        twist_vs_height_file.close()
+
+        if (plot_twist_vs_height):
+            ddic.plot_twist_vs_height(
+                working_folder=working_folder,
+                working_basename=working_basename,
+                suffix=None,
+                verbose=verbose)

image_expressions_cpp.py

@@ -21,6 +21,7 @@ def get_ExprIm_cpp(
         im_default_interpol_out_value="0.",
         grad_default_interpol_mode="linear", # linear, nearest
         grad_default_interpol_out_value="0.",
+        static_scaling_factor=0,
         verbose=0):
 
     assert (im_dim in (2,3))
@@ -121,9 +122,10 @@ public:
     {
         vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New();
         reader->SetFileName(filename);
-        reader->Update();
+        reader->Update();'''+('''
 
-        static_scaling = getStaticScalingFactor(reader->GetOutput()->GetScalarTypeAsString());'''+('''
+        static_scaling = getStaticScalingFactor(reader->GetOutput()->GetScalarTypeAsString());''')*(not static_scaling_factor)+('''
+        static_scaling = '''+str(static_scaling_factor)+''';''')*(static_scaling_factor)+('''
 
         vtkSmartPointer<vtkImageGradient> gradient = vtkSmartPointer<vtkImageGradient>::New();
 #if VTK_MAJOR_VERSION >= 6

plot_binned_strains_vs_radius.py

@@ -23,10 +23,10 @@ def plot_binned_strains_vs_radius(
         verbose=1):
 
     if (ref_folder is not None) and (ref_basename is not None):
-        lines = open(ref_folder+"/"+ref_basename+"-strains.dat").readlines()
+        lines = open(ref_folder+"/"+ref_basename+"-strains.dat").readlines()[1:]
     else:
-        lines = open(working_folder+"/"+working_basenames[0]+"-strains.dat").readlines()
-    n_frames = len(lines)-1
+        lines = open(working_folder+"/"+working_basenames[0]+"-strains.dat").readlines()[1:]
+    n_frames = len(lines)
 
     assert (components in ("all", "circ-long", "rad-circ"))
     if (components == "all"):

plot_twist_vs_height.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2019                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import os
+
+################################################################################
+
+def plot_twist_vs_height(
+        working_folder,
+        working_basename,
+        beta_range=15.,
+        suffix="",
+        verbose=1):
+
+    lines = open(working_folder+"/"+working_basename+"-strains.dat").readlines()[1:]
+    n_frames = len(lines)
+
+    if (suffix is None):
+        plotfile_basename = working_folder+"/"+working_basename+"-twist_vs_height"
+    else:
+        plotfile_basename = "plot_twist_vs_height"+("-"+suffix)*(suffix!="")
+    plotfile = open(plotfile_basename+".plt", "w")
+
+    plotfile.write('''\
+set terminal pdf enhanced
+
+load "Set1.plt"
+
+set output "'''+plotfile_basename+'''.pdf"
+
+set key off
+# set key box textcolor variable width +0
+
+set grid
+
+set xlabel "beta (deg)"
+xrange = '''+str(beta_range)+'''
+set xrange [-xrange:+xrange]
+
+# set ylabel "ll ()"
+set yrange [0.:1.]
+set ytics("apex" 0, "base" 1)
+
+# f(x) = a*x + b
+
+# g(x) = x/c - d/c
+# h(x) = c*x + d
+
+g(x) = c*x + d
+h(x) = x/c - d/c
+
+# FIT_MAXITER = 10
+
+datafile = "'''+working_folder+'''/'''+working_basename+'''-twist_vs_height.dat"
+
+''')
+    for k_frame in xrange(n_frames):
+        plotfile.write('''\
+# a = 1.
+# b = 0.5
+
+# c = 1.
+# d = 0.5
+
+c = 1.
+d = -0.5*c
+
+set title "index '''+str(k_frame)+'''"
+# fit f(x) datafile using 3:2 index '''+str(k_frame)+''' via a,b
+'''+('''# ''')*(k_frame==0)+'''fit g(x) datafile using 2:3 index '''+str(k_frame)+''' via c,d
+plot datafile using 3:2 index '''+str(k_frame)+''' notitle, h(x) notitle
+
+''')
+
+    plotfile.close()
+
+    os.system("gnuplot "+plotfile_basename+".plt")