diff --git a/__init__.py b/__init__.py
index 2424f480631bf5b7dd19c8d062907adbebb1426f..d1e6d57228be1c419dd367bd1d073258a27838a6 100644
--- a/__init__.py
+++ b/__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 *
diff --git a/compute_strains.py b/compute_strains.py
index a1f37c4c2e7d3290662e14e7378cc12414c41c98..9b74a884fea3a70db591594e2d6f8289fae50a0a 100644
--- a/compute_strains.py
+++ b/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)
diff --git a/image_expressions_cpp.py b/image_expressions_cpp.py
index 72b8586fb99585d5d1092551d037e8b1e2a6e5f2..754567bf6f9b0db8e313753f4fa1d6ce2f2fce71 100644
--- a/image_expressions_cpp.py
+++ b/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
diff --git a/plot_binned_strains_vs_radius.py b/plot_binned_strains_vs_radius.py
index 8d0033e959bb4783d2e2e2216af647a8e770917b..e0689da6f31743b85fa8c8046901bd939005b015 100644
--- a/plot_binned_strains_vs_radius.py
+++ b/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"):
diff --git a/plot_twist_vs_height.py b/plot_twist_vs_height.py
new file mode 100644
index 0000000000000000000000000000000000000000..1a41a5772e7b9dff30ab273b50a6e1bd38af7c4a
--- /dev/null
+++ b/plot_twist_vs_height.py
@@ -0,0 +1,84 @@
+#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")