diff --git a/__init__.py b/__init__.py
index 222a8f73c40a23f4d98e50fcb7e42f302b820a07..da8a2bd877937598ee4b2f16c302ea29c469346f 100644
--- a/__init__.py
+++ b/__init__.py
@@ -1,31 +1,34 @@
#this file was generated by __init__.py.py
-from compute_displacement_error import *
-from compute_displacements import *
-from compute_quadrature_degree import *
-from compute_strains import *
+from fedic import *
+from plot_binned_strains_vs_radius import *
+from NonlinearSolver import *
+from MeshSeries import *
+from generated_image_expressions_cpp import *
+from plot_strains_vs_radius import *
from compute_unwarped_images import *
+from Problem_ImageRegistration import *
+from compute_quadrature_degree import *
+from Energy_Regularization import *
from compute_warped_images import *
-from compute_warped_mesh import *
+from plot_strains import *
+from image_expressions_cpp import *
+from plot_twist_vs_height import *
from Energy import *
-from Energy_Regularization import *
+from compute_warped_mesh import *
from Energy_WarpedImage import *
-from fedic import *
-from fedic2 import *
-from generate_images import *
+from write_VTU_file import *
from generate_undersampled_images import *
-from generated_image_expressions_cpp import *
-from image_expressions_cpp import *
+from plot_regional_strains import *
from image_expressions_py import *
from ImageIterator import *
-from ImageSeries import *
-from MeshSeries import *
-from NonlinearSolver import *
-from plot_binned_strains_vs_radius import *
+from fedic2 import *
+from generated_grad_expressions_cpp import *
from plot_displacement_error import *
-from plot_regional_strains import *
-from plot_strains import *
-from plot_strains_vs_radius import *
from Problem import *
-from Problem_ImageRegistration import *
-from write_VTU_file import *
+from compute_displacements import *
+from generate_images import *
+from ImageSeries import *
+from compute_displacement_error import *
+from project import *
+from compute_strains import *
diff --git a/compute_strains.py b/compute_strains.py
index 8425d7d0ed1017014ee622239b55633b8e9d06cf..a9d82a97ca8667c996e18fab938b856494feb89f 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
@@ -30,48 +32,66 @@ def compute_strains(
ref_mesh_folder=None, # MG20190612: Mesh with sectors/parts/etc.
ref_mesh_basename=None,
ref_mesh_ext="vtk",
- remove_boundary_layer=True,
CYL_or_PPS="PPS",
+ remove_boundary_layer=True,
+ in_place=True,
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_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_part_id = ref_mesh.GetCellData().GetArray("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_part_id.GetTuple1(k_cell))
+ 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_part_id = None
+ 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
else:
@@ -93,12 +113,16 @@ 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 zz beta\n")
if (ref_frame is not None):
mesh0_filename = working_folder+"/"+working_basename+"_"+str(ref_frame).zfill(working_zfill)+"."+working_ext
@@ -118,16 +142,24 @@ 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 (iarray_part_id is not None):
- mesh.GetCellData().AddArray(iarray_part_id)
- if (iarray_sector_id is not None):
- mesh.GetCellData().AddArray(iarray_sector_id)
- if (write_strains_vs_radius) or (write_binned_strains_vs_radius):
+ 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("ll"))
myvtk.addDeformationGradients(
mesh=mesh,
disp_array_name=disp_array_name,
@@ -146,7 +178,9 @@ def compute_strains(
strain_array_name=strain_array_name,
mesh_w_local_basis=ref_mesh,
verbose=verbose)
- if (remove_boundary_layer):
+ 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",
@@ -155,18 +189,22 @@ def compute_strains(
threshold_by_upper_or_lower="lower")
n_points = mesh.GetNumberOfPoints()
n_cells = mesh.GetNumberOfCells()
- if (iarray_part_id is not None):
- iarray_part_id = mesh.GetCellData().GetArray("part_id")
- n_part_ids = 0
- if (iarray_sector_id is not None):
+ 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)
@@ -176,7 +214,7 @@ def compute_strains(
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_part_id.GetTuple1(k_cell) > 0]
+ 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)]
@@ -187,14 +225,15 @@ def compute_strains(
strains_per_sector[sector_id].append(farray_strain.GetTuple(k_cell))
else:
for k_cell in xrange(n_cells):
- part_id = int(iarray_part_id.GetTuple1(k_cell))
+ 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))
@@ -231,9 +270,60 @@ 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_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):
+ ll = farray_ll.GetTuple1(k_point)
+ if (ll < 1./3): continue
+ 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)
+ 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()
@@ -253,7 +343,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_basenames=[working_basename],
+ suffix=None,
+ verbose=verbose)
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..bc3d2da14b65ae999cc1663e874b7de73fd41548
--- /dev/null
+++ b/plot_twist_vs_height.py
@@ -0,0 +1,83 @@
+#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
+
+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
+
+# a = 1.
+# b = 0.5
+
+c = 1.
+d = 0.5
+
+# c = 1.
+# d = -0.5*c
+
+# FIT_MAXITER = 10
+
+datafile = "'''+working_folder+'''/'''+working_basename+'''-twist_vs_height.dat"
+
+''')
+ for k_frame in xrange(n_frames):
+ plotfile.write('''\
+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")