diff --git a/Energy.py b/Energy.py
index 064cca612fd9e3ac23e70fb2139937635c6e1ada..d84d3ceae41180d334cef1529735f12f1c899c3e 100644
--- a/Energy.py
+++ b/Energy.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/Energy_Regularization.py b/Energy_Regularization.py
index 651d89d813210528acfebed45ec15370954f6b26..1596d14bba92707c0faa4cca6839ff4c6295b62d 100644
--- a/Energy_Regularization.py
+++ b/Energy_Regularization.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/Energy_WarpedImage.py b/Energy_WarpedImage.py
index 454865a43e676fe073247e1a9de1c554daeead63..6ad2162dcdaa50ce6b97e6bab886b670bfdf7be4 100644
--- a/Energy_WarpedImage.py
+++ b/Energy_WarpedImage.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
@@ -30,6 +30,7 @@ class WarpedImageEnergy(Energy):
name="im",
w=1.,
ref_frame=0,
+ im_is_cone=0,
dynamic_scaling=False):
self.problem = problem
@@ -182,20 +183,32 @@ class WarpedImageEnergy(Energy):
# Phi_ref
self.Phi_Iref = dolfin.Expression(
cppcode=ddic.get_ExprCharFuncIm_cpp(
- im_dim=self.image_series.dimension),
+ im_dim=self.image_series.dimension,
+ im_is_def=0,
+ im_is_cone=im_is_cone),
element=self.fe)
self.Phi_Iref.init_image(self.ref_image_filename)
+ # Phi_def
+ self.Phi_Idef = dolfin.Expression(
+ cppcode=ddic.get_ExprCharFuncIm_cpp(
+ im_dim=self.image_series.dimension,
+ im_is_def=1,
+ im_is_cone=im_is_cone),
+ element=self.fe)
+ self.Phi_Idef.init_image(self.ref_image_filename)
+ self.Phi_Idef.init_disp(self.problem.U)
+
# Psi_c
- self.Psi_c = self.Phi_Iref * (self.Idef - self.Iref)**2/2
- self.DPsi_c = self.Phi_Iref * (self.Idef - self.Iref) * dolfin.dot(self.DIdef, self.problem.dU_test)
+ self.Psi_c = self.Phi_Idef * self.Phi_Iref * (self.Idef - self.Iref)**2/2
+ self.DPsi_c = self.Phi_Idef * self.Phi_Iref * (self.Idef - self.Iref) * dolfin.dot(self.DIdef, self.problem.dU_test)
- self.DDPsi_c = self.Phi_Iref * dolfin.dot(self.DIdef, self.problem.dU_trial) * dolfin.dot(self.DIdef, self.problem.dU_test)
- self.DDPsi_c_old = self.Phi_Iref * dolfin.dot(self.DIold, self.problem.dU_trial) * dolfin.dot(self.DIold, self.problem.dU_test)
- self.DDPsi_c_ref = self.Phi_Iref * dolfin.dot(self.DIref, self.problem.dU_trial) * dolfin.dot(self.DIref, self.problem.dU_test)
+ self.DDPsi_c = self.Phi_Idef * self.Phi_Iref * dolfin.dot(self.DIdef, self.problem.dU_trial) * dolfin.dot(self.DIdef, self.problem.dU_test)
+ self.DDPsi_c_old = self.Phi_Idef * self.Phi_Iref * dolfin.dot(self.DIold, self.problem.dU_trial) * dolfin.dot(self.DIold, self.problem.dU_test)
+ self.DDPsi_c_ref = self.Phi_Idef * self.Phi_Iref * dolfin.dot(self.DIref, self.problem.dU_trial) * dolfin.dot(self.DIref, self.problem.dU_test)
- self.Psi_c_old = self.Phi_Iref * (self.Idef - self.Iold)**2/2
- self.DPsi_c_old = self.Phi_Iref * (self.Idef - self.Iold) * dolfin.dot(self.DIdef, self.problem.dU_test)
+ self.Psi_c_old = self.Phi_Idef * self.Phi_Iref * (self.Idef - self.Iold)**2/2
+ self.DPsi_c_old = self.Phi_Idef * self.Phi_Iref * (self.Idef - self.Iold) * dolfin.dot(self.DIdef, self.problem.dU_test)
# forms
self.ener_form = self.Psi_c * self.dV
diff --git a/ImageIterator.py b/ImageIterator.py
index 07e1072df59bcc6cc974e7a3707a5058ac2a4ee5..378d6161066de73bbc00990ff782a0d6a14a7752 100644
--- a/ImageIterator.py
+++ b/ImageIterator.py
@@ -2,15 +2,18 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
################################################################################
+import dolfin
import glob
+import numpy
import os
import time
+import vtk
import myPythonLibrary as mypy
import myVTKPythonLibrary as myvtk
@@ -35,6 +38,11 @@ class ImageIterator():
self.working_folder = parameters["working_folder"] if ("working_folder" in parameters) else "."
self.working_basename = parameters["working_basename"] if ("working_basename" in parameters) else "sol"
self.register_ref_frame = parameters["register_ref_frame"] if ("register_ref_frame" in parameters) else False
+ self.initialize_U_from_file = parameters["initialize_U_from_file"] if ("initialize_U_from_file" in parameters) else False
+ self.initialize_U_folder = parameters["initialize_U_folder"] if ("initialize_U_folder" in parameters) else "."
+ self.initialize_U_basename = parameters["initialize_U_basename"] if ("initialize_U_basename" in parameters) else None
+ self.initialize_U_ext = parameters["initialize_U_ext"] if ("initialize_U_ext" in parameters) else "vtu"
+ self.initialize_U_array_name = parameters["initialize_U_array_name"] if ("initialize_U_array_name" in parameters) else "displacement"
self.initialize_DU_with_DUold = parameters["initialize_DU_with_DUold"] if ("initialize_DU_with_DUold" in parameters) else False
@@ -77,6 +85,15 @@ class ImageIterator():
self.printer.print_str("Looping over frames…")
+ if (self.initialize_U_from_file):
+ mesh_series = ddic.MeshSeries(
+ problem=self.problem,
+ folder=self.initialize_U_folder,
+ basename=self.initialize_U_basename,
+ ext=self.initialize_U_ext)
+
+ dof_to_vertex_map = dolfin.dof_to_vertex_map(self.problem.U_fs)
+
n_iter_tot = 0
global_success = True
for forward_or_backward in ["forward","backward"]:
@@ -108,7 +125,18 @@ class ImageIterator():
k_frame_old = k_frame+1
#self.printer.print_var("k_frame_old",k_frame_old,-1)
- if (self.initialize_DU_with_DUold):
+ if (self.initialize_U_from_file):
+ mesh = mesh_series.get_mesh(k_frame)
+ array_U = mesh.GetPointData().GetArray(self.initialize_U_array_name)
+ array_U = vtk.util.numpy_support.vtk_to_numpy(array_U)[:,:self.problem.mesh_dimension]
+ # print array_U
+ # array_U = array_U.astype(float)
+ # print array_U
+ array_U = numpy.reshape(array_U, array_U.size)
+ # print array_U
+ self.problem.U.vector()[:] = array_U[dof_to_vertex_map]
+
+ elif (self.initialize_DU_with_DUold):
self.problem.U.vector().axpy(1., self.problem.DUold.vector())
self.problem.call_before_solve(
diff --git a/ImageSeries.py b/ImageSeries.py
index 89eea3b992c6fd2b86b054d77f3fe6051043f2ef..c3c219f2d6ea2720acbdbccc64258493c04e50a7 100644
--- a/ImageSeries.py
+++ b/ImageSeries.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000000000000000000000000000000000000..65c5ca88a67c30becee01c5a8816d964b03862f9
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,165 @@
+ GNU LESSER GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+
+ This version of the GNU Lesser General Public License incorporates
+the terms and conditions of version 3 of the GNU General Public
+License, supplemented by the additional permissions listed below.
+
+ 0. Additional Definitions.
+
+ As used herein, "this License" refers to version 3 of the GNU Lesser
+General Public License, and the "GNU GPL" refers to version 3 of the GNU
+General Public License.
+
+ "The Library" refers to a covered work governed by this License,
+other than an Application or a Combined Work as defined below.
+
+ An "Application" is any work that makes use of an interface provided
+by the Library, but which is not otherwise based on the Library.
+Defining a subclass of a class defined by the Library is deemed a mode
+of using an interface provided by the Library.
+
+ A "Combined Work" is a work produced by combining or linking an
+Application with the Library. The particular version of the Library
+with which the Combined Work was made is also called the "Linked
+Version".
+
+ The "Minimal Corresponding Source" for a Combined Work means the
+Corresponding Source for the Combined Work, excluding any source code
+for portions of the Combined Work that, considered in isolation, are
+based on the Application, and not on the Linked Version.
+
+ The "Corresponding Application Code" for a Combined Work means the
+object code and/or source code for the Application, including any data
+and utility programs needed for reproducing the Combined Work from the
+Application, but excluding the System Libraries of the Combined Work.
+
+ 1. Exception to Section 3 of the GNU GPL.
+
+ You may convey a covered work under sections 3 and 4 of this License
+without being bound by section 3 of the GNU GPL.
+
+ 2. Conveying Modified Versions.
+
+ If you modify a copy of the Library, and, in your modifications, a
+facility refers to a function or data to be supplied by an Application
+that uses the facility (other than as an argument passed when the
+facility is invoked), then you may convey a copy of the modified
+version:
+
+ a) under this License, provided that you make a good faith effort to
+ ensure that, in the event an Application does not supply the
+ function or data, the facility still operates, and performs
+ whatever part of its purpose remains meaningful, or
+
+ b) under the GNU GPL, with none of the additional permissions of
+ this License applicable to that copy.
+
+ 3. Object Code Incorporating Material from Library Header Files.
+
+ The object code form of an Application may incorporate material from
+a header file that is part of the Library. You may convey such object
+code under terms of your choice, provided that, if the incorporated
+material is not limited to numerical parameters, data structure
+layouts and accessors, or small macros, inline functions and templates
+(ten or fewer lines in length), you do both of the following:
+
+ a) Give prominent notice with each copy of the object code that the
+ Library is used in it and that the Library and its use are
+ covered by this License.
+
+ b) Accompany the object code with a copy of the GNU GPL and this license
+ document.
+
+ 4. Combined Works.
+
+ You may convey a Combined Work under terms of your choice that,
+taken together, effectively do not restrict modification of the
+portions of the Library contained in the Combined Work and reverse
+engineering for debugging such modifications, if you also do each of
+the following:
+
+ a) Give prominent notice with each copy of the Combined Work that
+ the Library is used in it and that the Library and its use are
+ covered by this License.
+
+ b) Accompany the Combined Work with a copy of the GNU GPL and this license
+ document.
+
+ c) For a Combined Work that displays copyright notices during
+ execution, include the copyright notice for the Library among
+ these notices, as well as a reference directing the user to the
+ copies of the GNU GPL and this license document.
+
+ d) Do one of the following:
+
+ 0) Convey the Minimal Corresponding Source under the terms of this
+ License, and the Corresponding Application Code in a form
+ suitable for, and under terms that permit, the user to
+ recombine or relink the Application with a modified version of
+ the Linked Version to produce a modified Combined Work, in the
+ manner specified by section 6 of the GNU GPL for conveying
+ Corresponding Source.
+
+ 1) Use a suitable shared library mechanism for linking with the
+ Library. A suitable mechanism is one that (a) uses at run time
+ a copy of the Library already present on the user's computer
+ system, and (b) will operate properly with a modified version
+ of the Library that is interface-compatible with the Linked
+ Version.
+
+ e) Provide Installation Information, but only if you would otherwise
+ be required to provide such information under section 6 of the
+ GNU GPL, and only to the extent that such information is
+ necessary to install and execute a modified version of the
+ Combined Work produced by recombining or relinking the
+ Application with a modified version of the Linked Version. (If
+ you use option 4d0, the Installation Information must accompany
+ the Minimal Corresponding Source and Corresponding Application
+ Code. If you use option 4d1, you must provide the Installation
+ Information in the manner specified by section 6 of the GNU GPL
+ for conveying Corresponding Source.)
+
+ 5. Combined Libraries.
+
+ You may place library facilities that are a work based on the
+Library side by side in a single library together with other library
+facilities that are not Applications and are not covered by this
+License, and convey such a combined library under terms of your
+choice, if you do both of the following:
+
+ a) Accompany the combined library with a copy of the same work based
+ on the Library, uncombined with any other library facilities,
+ conveyed under the terms of this License.
+
+ b) Give prominent notice with the combined library that part of it
+ is a work based on the Library, and explaining where to find the
+ accompanying uncombined form of the same work.
+
+ 6. Revised Versions of the GNU Lesser General Public License.
+
+ The Free Software Foundation may publish revised and/or new versions
+of the GNU Lesser General Public License from time to time. Such new
+versions will be similar in spirit to the present version, but may
+differ in detail to address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Library as you received it specifies that a certain numbered version
+of the GNU Lesser General Public License "or any later version"
+applies to it, you have the option of following the terms and
+conditions either of that published version or of any later version
+published by the Free Software Foundation. If the Library as you
+received it does not specify a version number of the GNU Lesser
+General Public License, you may choose any version of the GNU Lesser
+General Public License ever published by the Free Software Foundation.
+
+ If the Library as you received it specifies that a proxy can decide
+whether future versions of the GNU Lesser General Public License shall
+apply, that proxy's public statement of acceptance of any version is
+permanent authorization for you to choose that version for the
+Library.
diff --git a/MeshSeries.py b/MeshSeries.py
new file mode 100644
index 0000000000000000000000000000000000000000..e61c69241269a811c66fe8ca755fb85b40856b5c
--- /dev/null
+++ b/MeshSeries.py
@@ -0,0 +1,73 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2018 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import glob
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+
+
+class MeshSeries():
+
+
+
+ def __init__(self,
+ problem,
+ folder,
+ basename,
+ n_frames=None,
+ ext="vtu"):
+
+ self.problem = problem
+ self.printer = self.problem.printer
+ self.folder = folder
+ self.basename = basename
+ self.n_frames = n_frames
+ self.ext = ext
+
+ self.printer.print_str("Reading mesh series…")
+ self.printer.inc()
+
+ self.filenames = glob.glob(self.folder+"/"+self.basename+"_[0-9]*"+"."+self.ext)
+ assert (len(self.filenames) >= 1),\
+ "Not enough meshes ("+self.folder+"/"+self.basename+"_[0-9]*"+"."+self.ext+"). Aborting."
+
+ if (self.n_frames is None):
+ self.n_frames = len(self.filenames)
+ else:
+ assert (self.n_frames <= len(self.filenames))
+ assert (self.n_frames >= 1),\
+ "n_frames = "+str(self.n_frames)+" < 2. Aborting."
+ self.printer.print_var("n_frames",self.n_frames)
+
+ self.zfill = len(self.filenames[0].rsplit("_",1)[-1].split(".",1)[0])
+
+ self.printer.dec()
+
+
+
+ def get_mesh_filename(self,
+ k_frame):
+
+ return self.folder+"/"+self.basename+"_"+str(k_frame).zfill(self.zfill)+"."+self.ext
+
+
+
+ def get_mesh(self,
+ k_frame):
+
+ return myvtk.readDataSet(
+ filename=self.get_mesh_filename(
+ k_frame=k_frame))
diff --git a/NonlinearSolver.py b/NonlinearSolver.py
index c983d8dbc25274030d3f1af97b4953c3a2f2ea15..959741fab26b2df58b21cd1ac2c73f3005faa559 100644
--- a/NonlinearSolver.py
+++ b/NonlinearSolver.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
@@ -322,7 +322,7 @@ class NonlinearSolver():
self.printer.print_var("dener_min",dener_min)
relax_err = dener_min/ener_list[0]
self.printer.print_var("relax_err",relax_err)
- if (relax_min != 0.) and (relax_err < self.relax_tol):
+ if (relax_min != 0.) and (abs(relax_err) < self.relax_tol):
break
if (relax_k == self.relax_n_iter_max):
break
diff --git a/Problem.py b/Problem.py
index cdbba9bc1076127ebeb4b5fd8467c4274c38d03a..a4028765883fc4571e4de301aaf64dd3830ac8bb 100644
--- a/Problem.py
+++ b/Problem.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/Problem_ImageRegistration.py b/Problem_ImageRegistration.py
index d0898f9edd6cf8346c5add70472f96fb5b00ebb5..f4aa4f23b9cec6dfb54e2871a9d4181ec557c94b 100644
--- a/Problem_ImageRegistration.py
+++ b/Problem_ImageRegistration.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
@@ -68,7 +68,7 @@ class ImageRegistrationProblem(Problem):
self.mesh_filebasename = self.mesh_folder+"/"+self.mesh_basename
self.mesh_filename = self.mesh_filebasename+"."+"xml"
assert (os.path.exists(self.mesh_filename)),\
- "No mesh in "+mesh_filename+". Aborting."
+ "No mesh in "+self.mesh_filename+". Aborting."
self.mesh = dolfin.Mesh(self.mesh_filename)
else:
self.mesh = mesh
@@ -192,14 +192,22 @@ class ImageRegistrationProblem(Problem):
def assemble_ener(self):
- ener_form = 0.
+ ener = 0.
for energy in self.energies:
- ener_form += dolfin.Constant(energy.w) * energy.ener_form
-
- ener = dolfin.assemble(
- ener_form)
+ ener_ = dolfin.assemble(
+ energy.ener_form)
+ self.printer.print_var("ener_"+energy.name,ener_)
+ ener += energy.w * ener_
#self.printer.print_var("ener",ener)
+ # ener_form = 0.
+ # for energy in self.energies:
+ # ener_form += dolfin.Constant(energy.w) * energy.ener_form
+ #
+ # ener = dolfin.assemble(
+ # ener_form)
+ # #self.printer.print_var("ener",ener)
+
return ener
diff --git a/__init__.py b/__init__.py
index 9deddd46220b2a07376a6436b44c404b8e61fc0f..222a8f73c40a23f4d98e50fcb7e42f302b820a07 100644
--- a/__init__.py
+++ b/__init__.py
@@ -1,30 +1,31 @@
#this file was generated by __init__.py.py
-from fedic2 import *
-from NonlinearSolver import *
-from compute_warped_mesh import *
-from image_expressions_py import *
-from compute_warped_images import *
+from compute_displacement_error import *
+from compute_displacements import *
+from compute_quadrature_degree import *
+from compute_strains import *
from compute_unwarped_images import *
+from compute_warped_images import *
+from compute_warped_mesh import *
+from Energy import *
+from Energy_Regularization import *
+from Energy_WarpedImage import *
from fedic import *
+from fedic2 import *
from generate_images import *
-from plot_binned_strains_vs_radius import *
-from plot_regional_strains import *
-from ImageIterator import *
-from Problem import *
-from Energy_WarpedImage import *
from generate_undersampled_images import *
-from Energy_Regularization import *
-from compute_strains import *
-from compute_displacements import *
-from compute_displacement_error import *
-from ImageSeries import *
-from compute_quadrature_degree import *
from generated_image_expressions_cpp import *
-from write_VTU_file import *
-from Energy import *
from image_expressions_cpp import *
-from Problem_ImageRegistration import *
-from plot_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 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 *
diff --git a/compute_displacement_error.py b/compute_displacement_error.py
index 8e7b6acef8b31dd91f13539e33ee576436b1f34d..74afebb4d0aececa71b7e612315fc99fb533ea2b 100644
--- a/compute_displacement_error.py
+++ b/compute_displacement_error.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/compute_displacements.py b/compute_displacements.py
index 17332be5a0a4d7f1a45f20ffac4f3aa3705687cc..8b539df001f97dbcfcf6d08b6fdcac08998d24d1 100644
--- a/compute_displacements.py
+++ b/compute_displacements.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/compute_quadrature_degree.py b/compute_quadrature_degree.py
index 78755d55ebbca9ecfc145024a437b94374b02cc4..e507fc5e5b3e9dd93471ec1c476ebf4da6564c94 100644
--- a/compute_quadrature_degree.py
+++ b/compute_quadrature_degree.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/compute_strains.py b/compute_strains.py
index 6c69d3c810d36730bcdaad3a2fa9e074a26eea66..a1f37c4c2e7d3290662e14e7378cc12414c41c98 100644
--- a/compute_strains.py
+++ b/compute_strains.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
@@ -46,7 +46,7 @@ def compute_strains(
verbose=verbose)
ref_mesh_n_cells = ref_mesh.GetNumberOfCells()
if (verbose): print "ref_mesh_n_cells = " + str(ref_mesh_n_cells)
-
+
if (ref_mesh.GetCellData().HasArray("sector_id")):
iarray_sector_id = ref_mesh.GetCellData().GetArray("sector_id")
n_sector_ids = 0
@@ -57,7 +57,14 @@ def compute_strains(
n_sector_ids = sector_id+1
if (verbose): print "n_sector_ids = " + str(n_sector_ids)
else:
+ iarray_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_sector_ids = 0
@@ -104,10 +111,10 @@ def compute_strains(
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 (ref_mesh.GetCellData().HasArray("part_id")):
- mesh.GetCellData().AddArray(ref_mesh.GetCellData().GetArray("part_id"))
- if (ref_mesh.GetCellData().HasArray("sector_id")):
- mesh.GetCellData().AddArray(ref_mesh.GetCellData().GetArray("sector_id"))
+ 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)
myvtk.addDeformationGradients(
mesh=mesh,
disp_array_name=disp_array_name,
diff --git a/compute_unwarped_images.py b/compute_unwarped_images.py
index cd73bdfed22fc12d0593de0473617b959d2d8e41..1e4c832b5670df7aaf97bca65072c6f52988345f 100644
--- a/compute_unwarped_images.py
+++ b/compute_unwarped_images.py
@@ -2,10 +2,10 @@
################################################################################
### ###
-### Created by Martin Genet, 2012-2016 ###
+### Created by Martin Genet, 2012-2019 ###
### ###
-### University of California at San Francisco (UCSF), USA ###
-### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
+### University of California at San Francisco (UCSF), USA ###
+### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
### École Polytechnique, Palaiseau, France ###
### ###
################################################################################
diff --git a/compute_warped_images.py b/compute_warped_images.py
index db65cd05f475aae8e3608828cb7ee34315befa6b..8e0d366e6e566c5bdf127f08198c5fa8d27e0e46 100644
--- a/compute_warped_images.py
+++ b/compute_warped_images.py
@@ -2,10 +2,10 @@
################################################################################
### ###
-### Created by Martin Genet, 2012-2016 ###
+### Created by Martin Genet, 2012-2019 ###
### ###
-### University of California at San Francisco (UCSF), USA ###
-### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
+### University of California at San Francisco (UCSF), USA ###
+### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
### École Polytechnique, Palaiseau, France ###
### ###
################################################################################
@@ -27,7 +27,10 @@ def compute_warped_images(
working_folder,
working_basename,
ref_frame=0,
+ ref_image_model=None,
working_ext="vtu",
+ working_displacement_field_name="displacement",
+ print_warped_mesh=0,
verbose=0):
ref_image_zfill = len(glob.glob(ref_image_folder+"/"+ref_image_basename+"_*.vti")[0].rsplit("_")[-1].split(".")[0])
@@ -35,10 +38,11 @@ def compute_warped_images(
ref_image = myvtk.readImage(
filename=ref_image_filename)
- interpolator = myvtk.getImageInterpolator(
- image=ref_image)
- #I = numpy.empty(1)
- #interpolator.Interpolate([0.35, 0.25, 0.], I)
+ if (ref_image_model is None):
+ ref_image_interpolator = myvtk.getImageInterpolator(
+ image=ref_image)
+ #I = numpy.empty(1)
+ #ref_image_interpolator.Interpolate([0.35, 0.25, 0.], I)
image = vtk.vtkImageData()
image.SetOrigin(ref_image.GetOrigin())
@@ -54,7 +58,7 @@ def compute_warped_images(
working_zfill = len(glob.glob(working_folder+"/"+working_basename+"_*."+working_ext)[0].rsplit("_")[-1].split(".")[0])
n_frames = len(glob.glob(working_folder+"/"+working_basename+"_"+"[0-9]"*working_zfill+"."+working_ext))
- #n_frames = 1
+ # n_frames = 1
X = numpy.empty(3)
U = numpy.empty(3)
@@ -66,7 +70,9 @@ def compute_warped_images(
mesh = myvtk.readUGrid(
filename=working_folder+"/"+working_basename+"_"+str(k_frame).zfill(working_zfill)+"."+working_ext)
- #print mesh
+ # print mesh
+ assert (mesh.GetPointData().HasArray(working_displacement_field_name)), "no array '" + working_displacement_field_name + "' in mesh"
+ mesh.GetPointData().SetActiveVectors(working_displacement_field_name)
warp = vtk.vtkWarpVector()
if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
@@ -75,9 +81,10 @@ def compute_warped_images(
warp.SetInput(mesh)
warp.Update()
warped_mesh = warp.GetOutput()
- #myvtk.writeUGrid(
- #ugrid=warped_mesh,
- #filename=working_folder+"/"+working_basename+"-warped_"+str(k_frame).zfill(working_zfill)+"."+working_ext)
+ if print_warped_mesh:
+ myvtk.writeUGrid(
+ ugrid=warped_mesh,
+ filename=working_folder+"/"+working_basename+"-warped_"+str(k_frame).zfill(working_zfill)+"."+working_ext)
probe = vtk.vtkProbeFilter()
if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
@@ -89,7 +96,7 @@ def compute_warped_images(
probe.Update()
probed_image = probe.GetOutput()
scalars_mask = probed_image.GetPointData().GetArray("vtkValidPointMask")
- scalars_U = probed_image.GetPointData().GetArray("displacement")
+ scalars_U = probed_image.GetPointData().GetArray(working_displacement_field_name)
#myvtk.writeImage(
#image=probed_image,
#filename=working_folder+"/"+working_basename+"_"+str(k_frame).zfill(working_zfill)+".vti")
@@ -102,7 +109,10 @@ def compute_warped_images(
image.GetPoint(k_point, x)
scalars_U.GetTuple(k_point, U)
X = x - U
- interpolator.Interpolate(X, I)
+ if (ref_image_model is None):
+ ref_image_interpolator.Interpolate(X, I)
+ else:
+ I[0] = ref_image_model(X)
scalars.SetTuple(k_point, I)
myvtk.writeImage(
diff --git a/compute_warped_mesh.py b/compute_warped_mesh.py
index e1968c9e6512214628f8c8799417e60d6f965c53..cbe4af546cf9929e8d859a2c6f2e833daeafbcd6 100644
--- a/compute_warped_mesh.py
+++ b/compute_warped_mesh.py
@@ -2,10 +2,10 @@
################################################################################
### ###
-### Created by Martin Genet, 2012-2016 ###
+### Created by Martin Genet, 2012-2019 ###
### ###
-### University of California at San Francisco (UCSF), USA ###
-### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
+### University of California at San Francisco (UCSF), USA ###
+### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
### École Polytechnique, Palaiseau, France ###
### ###
################################################################################
@@ -72,7 +72,7 @@ def compute_warped_mesh(
myvtk.addStrainsFromDisplacements(
mesh=mesh,
disp_array_name="displacement",
- ref_mesh=ref_mesh,
+ mesh_w_local_basis=ref_mesh,
verbose=verbose-1)
myvtk.writeUGrid(
diff --git a/fedic.py b/fedic.py
index ee505f8193732d74952ae7c3d9f1207472a357af..d7bfac8586cb43c885905babdc8a649bbbb0c1bf 100644
--- a/fedic.py
+++ b/fedic.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/fedic2.py b/fedic2.py
index ffe4b66b185f2c02de545ae6f8c3ee97d58f8f1e..a3b5f2a00b164cf9045b00e8278441aec7e29086 100644
--- a/fedic2.py
+++ b/fedic2.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
@@ -28,6 +28,7 @@ def fedic2(
images_quadrature_from="points_count", # points_count, integral
images_expressions_type="cpp", # cpp, py
images_dynamic_scaling=1,
+ images_is_cone=0,
mesh=None,
mesh_folder=None,
mesh_basename=None,
@@ -41,6 +42,11 @@ def fedic2(
residual_type="Iref", # Iref, Iold, Iref-then-Iold
relax_type="gss", # constant, aitken, gss
relax_init=1.0,
+ initialize_U_from_file=0,
+ initialize_U_folder=None,
+ initialize_U_basename=None,
+ initialize_U_ext="vtu",
+ initialize_U_array_name="displacement",
initialize_DU_with_DUold=0,
tol_res=None,
tol_res_rel=None,
@@ -59,6 +65,8 @@ def fedic2(
"residual_type must be \"Iref\". Aborting."
assert (relax_init == 1.),\
"relax_init must be 1. Aborting."
+ assert (not ((initialize_U_from_file) and (initialize_DU_with_DUold))),\
+ "Cannot initialize U from file and DU with DUold together. Aborting."
assert (tol_res is None),\
"tol_res is deprecated. Aborting."
assert (tol_im is None),\
@@ -83,6 +91,8 @@ def fedic2(
ext=images_ext)
if (images_quadrature is None):
+ problem.printer.print_str("Computing quadrature degree…")
+ problem.printer.inc()
if (images_quadrature_from == "points_count"):
images_quadrature = ddic.compute_quadrature_degree_from_points_count(
image_filename=image_series.get_image_filename(images_ref_frame),
@@ -95,6 +105,8 @@ def fedic2(
verbose=1)
else:
assert (0), "\"images_quadrature_from\" (="+str(images_quadrature_from)+") must be \"points_count\" or \"integral\". Aborting."
+ problem.printer.print_var("images_quadrature",images_quadrature)
+ problem.printer.dec()
warped_image_energy = ddic.WarpedImageEnergy(
problem=problem,
@@ -102,7 +114,8 @@ def fedic2(
quadrature_degree=images_quadrature,
w=1.-regul_level,
ref_frame=images_ref_frame,
- dynamic_scaling=images_dynamic_scaling)
+ dynamic_scaling=images_dynamic_scaling,
+ im_is_cone=images_is_cone)
problem.add_image_energy(warped_image_energy)
regularization_energy = ddic.RegularizationEnergy(
@@ -131,6 +144,11 @@ def fedic2(
parameters={
"working_folder":working_folder,
"working_basename":working_basename,
+ "initialize_U_from_file":initialize_U_from_file,
+ "initialize_U_folder":initialize_U_folder,
+ "initialize_U_basename":initialize_U_basename,
+ "initialize_U_ext":initialize_U_ext,
+ "initialize_U_array_name":initialize_U_array_name,
"initialize_DU_with_DUold":initialize_DU_with_DUold})
image_iterator.iterate()
diff --git a/generate_images.py b/generate_images.py
index fe2a56f9a64679ed4efea8717598bfa7ce2879cf..ffcf5017505b02462e51095b3336ea617e00016c 100644
--- a/generate_images.py
+++ b/generate_images.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
@@ -103,6 +103,10 @@ class Image():
# structure
if (structure["type"] == "no"):
self.I0_structure = self.I0_structure_no
+ elif (structure["type"] == "box"):
+ self.I0_structure = self.I0_structure_box
+ self.Xmin = structure["Xmin"]
+ self.Xmax = structure["Xmax"]
elif (structure["type"] == "heart"):
if (images["n_dim"] == 2):
self.I0_structure = self.I0_structure_heart_2
@@ -201,6 +205,13 @@ class Image():
i[0] = 1.
if (g is not None): g[:] = 1. # MG 20180806: gradient is given by texture; here it is just indicator function
+ def I0_structure_box(self, X, i, g=None):
+ if all(numpy.greater_equal(X, self.Xmin)) and all(numpy.less_equal(X, self.Xmax)):
+ i[0] = 1.
+ else:
+ i[0] = 0.
+ if (g is not None): g[:] = 0. # MG 20180806: gradient is given by texture; here it is just indicator function
+
def I0_structure_heart_2(self, X, i, g=None):
self.R = ((X[0]-self.L[0]/2)**2 + (X[1]-self.L[1]/2)**2)**(1./2)
if (self.R >= self.Ri) and (self.R <= self.Re):
@@ -437,9 +448,10 @@ class Mapping:
pass
def init_t_translation(self, t):
- self.D[0] = self.deformation["Dx"]*self.phi(t) if ("Dx" in self.deformation.keys()) else 0.
- self.D[1] = self.deformation["Dy"]*self.phi(t) if ("Dy" in self.deformation.keys()) else 0.
- self.D[2] = self.deformation["Dz"]*self.phi(t) if ("Dz" in self.deformation.keys()) else 0.
+ self.D[0] = self.deformation["Dx"] if ("Dx" in self.deformation.keys()) else 0.
+ self.D[1] = self.deformation["Dy"] if ("Dy" in self.deformation.keys()) else 0.
+ self.D[2] = self.deformation["Dz"] if ("Dz" in self.deformation.keys()) else 0.
+ self.D *= self.phi(t)
def init_t_rotation(self, t):
self.C[0] = self.deformation["Cx"] if ("Cx" in self.deformation.keys()) else 0.
@@ -461,23 +473,31 @@ class Mapping:
self.Rinv[:,:] = numpy.linalg.inv(self.R)
def init_t_homogeneous(self, t):
- Exx = self.deformation["Exx"]*self.phi(t) if ("Exx" in self.deformation.keys()) else 0.
- Eyy = self.deformation["Eyy"]*self.phi(t) if ("Eyy" in self.deformation.keys()) else 0.
- Ezz = self.deformation["Ezz"]*self.phi(t) if ("Ezz" in self.deformation.keys()) else 0.
- Exy = self.deformation["Exy"]*self.phi(t) if ("Exy" in self.deformation.keys()) else 0.
- Eyx = self.deformation["Eyx"]*self.phi(t) if ("Eyx" in self.deformation.keys()) else 0.
- Exz = self.deformation["Exz"]*self.phi(t) if ("Exz" in self.deformation.keys()) else 0.
- Ezx = self.deformation["Ezx"]*self.phi(t) if ("Ezx" in self.deformation.keys()) else 0.
- Eyz = self.deformation["Eyz"]*self.phi(t) if ("Eyz" in self.deformation.keys()) else 0.
- Ezy = self.deformation["Ezy"]*self.phi(t) if ("Ezy" in self.deformation.keys()) else 0.
- self.F = numpy.array([[Exx, Exy, Exz],
- [Eyx, Eyy, Eyz],
- [Ezx, Ezy, Ezz]])
- self.F *= 2
- self.F += numpy.eye(3)
- w, v = numpy.linalg.eig(self.F)
- #assert (numpy.diag(numpy.dot(numpy.dot(numpy.transpose(v), self.F), v)) == w).all(), str(numpy.dot(numpy.dot(numpy.transpose(v), self.F), v))+" ≠"+str(numpy.diag(w))+". Aborting."
- self.F = numpy.dot(numpy.dot(v, numpy.diag(numpy.sqrt(w))), numpy.transpose(v))
+ if (any(E in self.deformation.keys() for E in ("Exx", "Eyy", "Ezz"))): # build F from E
+ Exx = self.deformation["Exx"] if ("Exx" in self.deformation.keys()) else 0.
+ Eyy = self.deformation["Eyy"] if ("Eyy" in self.deformation.keys()) else 0.
+ Ezz = self.deformation["Ezz"] if ("Ezz" in self.deformation.keys()) else 0.
+ self.F = numpy.array([[Exx, 0., 0.],
+ [ 0., Eyy, 0.],
+ [ 0., 0., Ezz]])*self.phi(t)
+ self.F *= 2
+ self.F += numpy.eye(3)
+ w, v = numpy.linalg.eig(self.F)
+ # assert (numpy.diag(numpy.dot(numpy.dot(numpy.transpose(v), self.F), v)) == w).all(), str(numpy.dot(numpy.dot(numpy.transpose(v), self.F), v))+" ≠"+str(numpy.diag(w))+". Aborting."
+ self.F = numpy.dot(numpy.dot(v, numpy.diag(numpy.sqrt(w))), numpy.transpose(v))
+ else:
+ Fxx = self.deformation["Fxx"] if ("Fxx" in self.deformation.keys()) else 0.
+ Fyy = self.deformation["Fyy"] if ("Fyy" in self.deformation.keys()) else 0.
+ Fzz = self.deformation["Fzz"] if ("Fzz" in self.deformation.keys()) else 0.
+ Fxy = self.deformation["Fxy"] if ("Fxy" in self.deformation.keys()) else 0.
+ Fyx = self.deformation["Fyx"] if ("Fyx" in self.deformation.keys()) else 0.
+ Fyz = self.deformation["Fyz"] if ("Fyz" in self.deformation.keys()) else 0.
+ Fzy = self.deformation["Fzy"] if ("Fzy" in self.deformation.keys()) else 0.
+ Fzx = self.deformation["Fzx"] if ("Fzx" in self.deformation.keys()) else 0.
+ Fxz = self.deformation["Fxz"] if ("Fxz" in self.deformation.keys()) else 0.
+ self.F = numpy.eye(3) + (numpy.array([[Fxx, Fxy, Fxz],
+ [Fyx, Fyy, Fyz],
+ [Fzx, Fzy, Fzz]])-numpy.eye(3))*self.phi(t)
self.Finv = numpy.linalg.inv(self.F)
def init_t_heart(self, t):
diff --git a/generate_undersampled_images.py b/generate_undersampled_images.py
index 7deda2065669f79fd0344e841c06dfa30800d3c5..bcd23330727d1434c3800bd410a834aa49fa30b8 100644
--- a/generate_undersampled_images.py
+++ b/generate_undersampled_images.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/generated_grad_expressions_cpp.py b/generated_grad_expressions_cpp.py
new file mode 100644
index 0000000000000000000000000000000000000000..fc64124bbb8141e2c149c1d2feb565aa8c94a9f0
--- /dev/null
+++ b/generated_grad_expressions_cpp.py
@@ -0,0 +1,176 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2019 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+################################################################################
+
+def get_ExprGenGrad_cpp(
+ im_dim,
+ im_type="grad",
+ im_is_def=0,
+ disp_type="fenics", # "vtk"
+ verbose=0):
+
+ assert (im_dim in (2,3))
+ assert (im_type in ("grad"))
+
+ ExprGenGrad_cpp = '''\
+#include <string.h>
+
+#include <vtkSmartPointer.h>
+#include <vtkStructuredPointsReader.h>
+#include <vtkXMLImageDataReader.h>
+#include <vtkImageData.h>
+#include <vtkImageGradient.h>
+#include <vtkImageInterpolator.h>
+
+double getStaticScalingFactor(const char* scalar_type_as_string)
+{
+ if (strcmp(scalar_type_as_string, "unsigned char" ) == 0) return pow(2, 8)-1;
+ if (strcmp(scalar_type_as_string, "unsigned short") == 0) return pow(2, 16)-1;
+ if (strcmp(scalar_type_as_string, "unsigned int" ) == 0) return pow(2, 32)-1;
+ if (strcmp(scalar_type_as_string, "unsigned long" ) == 0) return pow(2, 64)-1;
+ if (strcmp(scalar_type_as_string, "float" ) == 0) return 1.;
+ if (strcmp(scalar_type_as_string, "double" ) == 0) return 1.;
+ assert (0);
+}
+
+namespace dolfin
+{
+
+class MyExpr : public Expression
+{
+ vtkSmartPointer<vtkImageData> image;
+ vtkSmartPointer<vtkImageGradient> gradient;
+ vtkSmartPointer<vtkImageInterpolator> interpolator;
+ double static_scaling;
+ std::shared_ptr<dolfin::Mesh> mesh;
+ std::shared_ptr<dolfin::Function> U;
+ mutable Array<double> UX;
+ mutable Array<double> x;
+
+public:
+
+ MyExpr():
+ Expression('''+str(im_dim)*(im_type in ("grad"))+'''),
+ UX('''+str(im_dim)+'''),
+ x(3),
+ interpolator(vtkSmartPointer<vtkImageInterpolator>::New())
+
+ {
+ }
+
+ void init_mesh(
+ const std::shared_ptr<dolfin::Mesh> mesh_)
+ {
+ mesh = mesh_;
+
+ std::cout << mesh->num_vertices() << std::endl;
+ std::cout << mesh->num_cells() << std::endl;
+ }
+
+ void init_disp(
+ const std::shared_ptr<dolfin::Function> U_)
+ {
+ U = U_;
+ }
+
+ void init_image()
+ {
+ image = vtkSmartPointer<vtkImageData>::New(); // MG20180913: another way to instantiate object
+ gradient = vtkSmartPointer<vtkImageGradient>::New();
+ }
+
+
+
+ void eval(Array<double>& expr, const Array<double>& X) const
+ {'''+('''
+
+ std::cout << "X = " << X.str(1) << std::endl;''')*(verbose)+'''
+
+ U->eval(UX, X);'''+('''
+ std::cout << "UX = " << UX.str(1) << std::endl;''')*(verbose)+('''
+ x[0] = X[0] + UX[0];
+ x[1] = X[1] + UX[1];''')*(im_dim==2)+('''
+ x[0] = X[0] + UX[0];
+ x[1] = X[1] + UX[1];
+ x[2] = X[2] + UX[2];''')*(im_dim==3)+('''
+ std::cout << "x = " << x.str(1) << std::endl;''')*(verbose)+'''
+
+ interpolator->Interpolate(x.data(), expr.data());'''+('''
+ std::cout << "expr = " << expr.str(1) << std::endl;''')*(verbose)+('''
+
+ expr[0] /= static_scaling;''')*(im_type=="im")+('''
+
+ expr[0] /= static_scaling;
+ expr[1] /= static_scaling;''')*(im_type=="grad")*(im_dim==2)+('''
+
+ expr[0] /= static_scaling;
+ expr[1] /= static_scaling;
+ expr[2] /= static_scaling;''')*(im_type=="grad")*(im_dim==3)+('''
+ std::cout << "expr = " << expr.str(1) << std::endl;''')*(verbose)+'''
+ }
+
+
+ void prepareImage( int n1, int n2, int n3, float origin1, float origin2, float origin3, float spac1, float spac2, float spac3)
+ {
+ image->SetDimensions(n1,n2,n3);
+ image->AllocateScalars(VTK_FLOAT,1);
+
+ image->SetSpacing(spac1, spac2, spac3);
+ image->SetOrigin(origin1, origin2, origin3);
+
+ static_scaling = getStaticScalingFactor("double");
+ }
+
+
+ void assignValues( int i, int j, int k, float value)
+ {
+ image->SetScalarComponentFromFloat(i,j,k,0,value);
+ }
+
+ void setGradient_Interpol()
+ {
+ interpolator->SetInterpolationModeToLinear();
+ interpolator->SetOutValue(.0);
+
+ gradient->SetInputData(image);
+ gradient->SetDimensionality('''+str(im_dim)+''');
+ gradient->Update();
+ interpolator->Initialize(gradient->GetOutput());
+
+ interpolator->Update();
+ }
+
+
+ void prepareImage_2D( int n1, int n2, float origin1, float origin2, float spac1, float spac2)
+ {
+ const double &Z=0.;
+ image->SetExtent(0, n1-1, 0, n2-1, 0, 0);
+
+ image->AllocateScalars(VTK_FLOAT,1);
+
+ image->SetSpacing(spac1, spac2, 1);
+ image->SetOrigin(origin1, origin2, 0);
+
+ static_scaling = getStaticScalingFactor("double");
+ x[2] = Z;
+
+ }
+
+
+ void assignValues_2D( int i, int j, float value)
+ {
+ image->SetScalarComponentFromFloat(i,j,0,0,value);
+ }
+};
+
+}'''
+ #print ExprIm_cpp
+ return ExprGenGrad_cpp
diff --git a/generated_image_expressions_cpp.py b/generated_image_expressions_cpp.py
index 814ea78cb56374a8d5a186fb30fe60e74fd6d55c..dbbef064dcae8e2dc297b4d3916acae035d97b20 100644
--- a/generated_image_expressions_cpp.py
+++ b/generated_image_expressions_cpp.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
@@ -61,8 +61,8 @@ public:
reader->SetFileName(filename);
reader->Update();
- image = reader->GetOutput();
// image = vtkSmartPointer<vtkImageData>::New(); // MG20180913: another way to instantiate object
+ image = reader->GetOutput();
static_scaling = getStaticScalingFactor(image->GetScalarTypeAsString());
diff --git a/image_expressions_cpp.py b/image_expressions_cpp.py
index 361fe93094eb34320c0ead07f6bf9cec2b04de09..9eb43bb04e367af99e061cac04e3ab890d2191fc 100644
--- a/image_expressions_cpp.py
+++ b/image_expressions_cpp.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
@@ -15,6 +15,7 @@ def get_ExprIm_cpp(
im_type="im",
im_is_def=0,
disp_type="fenics", # "vtk"
+ u_is_vtk=0,
verbose=0):
assert (im_dim in (2,3))
@@ -29,6 +30,12 @@ def get_ExprIm_cpp(
#include <vtkImageData.h>'''+('''
#include <vtkImageGradient.h>''')*(im_type=="grad")+'''
#include <vtkImageInterpolator.h>
+#include <vtkXMLUnstructuredGridReader.h>
+#include <vtkUnstructuredGrid.h>
+#include <vtkProbeFilter.h>
+#include <vtkPointData.h>
+#include <vtkPolyData.h>
+
double getStaticScalingFactor(const char* scalar_type_as_string)
{
@@ -50,8 +57,11 @@ class MyExpr : public Expression
double static_scaling;'''+('''
Array<double>* dynamic_scaling; // does not work
double dynamic_scaling_a; // should not be needed
- double dynamic_scaling_b; // should not be needed
- Function* U;
+ double dynamic_scaling_b; // should not be needed'''+('''
+ Function* U;''')*(not u_is_vtk)+('''
+ vtkSmartPointer<vtkProbeFilter> probe_filter;
+ vtkSmartPointer<vtkPoints> probe_points;
+ vtkSmartPointer<vtkPolyData> probe_polydata;''')*(u_is_vtk)+'''
mutable Array<double> UX;
mutable Array<double> x;''')*(im_is_def)+('''
mutable Array<double> X3D;''')*(not im_is_def)*(im_dim==2)+'''
@@ -65,7 +75,10 @@ public:
UX('''+str(im_dim)+'''),
x(3)''')*(im_is_def)+(''',
X3D(3)''')*(not im_is_def)*(im_dim==2)+'''
- {
+ {'''+('''
+ probe_filter = vtkSmartPointer<vtkProbeFilter>::New();
+ probe_points = vtkSmartPointer<vtkPoints>::New();
+ probe_polydata = vtkSmartPointer<vtkPolyData>::New();''')*(u_is_vtk)+'''
}'''+('''
void init_dynamic_scaling(
@@ -84,11 +97,23 @@ public:
dynamic_scaling_b = scaling[1]; // should not be needed
} // should not be needed
+ '''+('''
void init_disp(
Function* UU)
{
U = UU;
- }''')*(im_is_def)+'''
+ }''')*(not u_is_vtk)+('''
+ void init_disp(
+ const char* mesh_filename)
+ {
+ vtkSmartPointer<vtkXMLUnstructuredGridReader> reader = vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
+ reader->SetFileName(mesh_filename);
+ reader->Update();
+
+ vtkSmartPointer<vtkUnstructuredGrid> mesh = reader->GetOutput();
+
+ probe_filter->SetSourceData(mesh);
+ }''')*(u_is_vtk))*(im_is_def)+'''
void init_image(
const char* filename,
@@ -140,9 +165,18 @@ public:
void eval(Array<double>& expr, const Array<double>& X) const
{'''+('''
- std::cout << "X = " << X.str(1) << std::endl;''')*(verbose)+('''
+ std::cout << "X = " << X.str(1) << std::endl;''')*(verbose)+(('''
+
+ U->eval(UX, X);''')*(not u_is_vtk)+('''
+
+ probe_points->SetNumberOfPoints(1);
+ probe_points->SetPoint(0,X.data());
+ probe_polydata->SetPoints(probe_points);
+ probe_filter->SetInputData(probe_polydata);
+ probe_filter->Update();
+ probe_filter->GetOutput()->GetPointData()->GetArray("U")->GetTuple(0, UX.data());
- U->eval(UX, X);'''+('''
+ ''')*(u_is_vtk)+('''
std::cout << "UX = " << UX.str(1) << std::endl;''')*(verbose)+('''
x[0] = X[0] + UX[0];
x[1] = X[1] + UX[1];''')*(im_dim==2)+('''
@@ -209,13 +243,19 @@ public:
#print ExprIm_cpp
return ExprIm_cpp
+
def get_ExprCharFuncIm_cpp(
im_dim,
+ im_is_def=0,
+ im_is_cone=0,
verbose=0):
assert (im_dim in (2,3))
+ if (im_is_cone):
+ assert (im_dim == 3)
ExprCharFuncIm_cpp = '''\
+#include <math.h>
#include <string.h>
#include <vtkSmartPointer.h>
@@ -228,13 +268,35 @@ namespace dolfin
class MyExpr : public Expression
{
double xmin, xmax, ymin, ymax, zmin, zmax;
+ mutable Array<double> UX;
+ mutable Array<double> x;'''+('''
+ Function* U;''')*(im_is_def)+('''
+ Array<double> O;
+ Array<double> n1;
+ Array<double> n2;
+ Array<double> n3;
+ Array<double> n4;
+ mutable double d1, d2, d3, d4;''')*(im_is_cone)+'''
public:
MyExpr():
- Expression()
+ Expression(),
+ UX('''+str(im_dim)+'''),
+ x('''+str(im_dim)+''')'''+(''',
+ O(3),
+ n1(3),
+ n2(3),
+ n3(3),
+ n4(3)''')*(im_is_cone)+'''
{
- }
+ }'''+('''
+
+ void init_disp(
+ Function* UU)
+ {
+ U = UU;
+ }''')*(im_is_def)+'''
void init_image(
const char* filename)
@@ -257,16 +319,71 @@ public:
std::cout << "ymin = " << ymin << std::endl;
std::cout << "ymax = " << ymax << std::endl;
std::cout << "zmin = " << zmin << std::endl;
- std::cout << "zmax = " << zmax << std::endl;''')*(verbose)+'''
+ std::cout << "zmax = " << zmax << std::endl;''')*(verbose)+('''
+
+ O[0] = (xmin+xmax)/2;
+ O[1] = (ymin+ymax)/2;
+ O[2] = zmax;
+
+ n1[0] = +cos(35. * M_PI/180.);
+ n1[1] = 0.;
+ n1[2] = -sin(35. * M_PI/180.);
+
+ n2[0] = -cos(35. * M_PI/180.);
+ n2[1] = 0.;
+ n2[2] = -sin(35. * M_PI/180.);
+
+ n3[0] = 0.;
+ n3[1] = +cos(40. * M_PI/180.);
+ n3[2] = -sin(40. * M_PI/180.);
+
+ n4[0] = 0.;
+ n4[1] = -cos(40. * M_PI/180.);
+ n4[2] = -sin(40. * M_PI/180.);''')*(im_is_cone)+'''
}
void eval(Array<double>& expr, const Array<double>& X) const
{'''+('''
std::cout << "X = " << X.str(1) << std::endl;''')*(verbose)+('''
- if ((X[0] >= xmin) && (X[0] <= xmax) && (X[1] >= ymin) && (X[1] <= ymax))''')*(im_dim==2)+('''
+ U->eval(UX, X);''')*(im_is_def)+('''
+ std::cout << "UX = " << UX.str(1) << std::endl;''')*(verbose)+('''
- if ((X[0] >= xmin) && (X[0] <= xmax) && (X[1] >= ymin) && (X[1] <= ymax) && (X[2] >= zmin) && (X[2] <= zmax))''')*(im_dim==3)+'''
+ x[0] = X[0] + UX[0];
+ x[1] = X[1] + UX[1];''')*(im_dim==2)+('''
+ x[0] = X[0] + UX[0];
+ x[1] = X[1] + UX[1];
+ x[2] = X[2] + UX[2];''')*(im_dim==3)+('''
+ std::cout << "x = " << x.str(1) << std::endl;''')*(verbose)+(('''
+
+ if ((x[0] >= xmin)
+ && (x[0] <= xmax)
+ && (x[1] >= ymin)
+ && (x[1] <= ymax))''')*(im_dim==2)+('''
+ if ((x[0] >= xmin)
+ && (x[0] <= xmax)
+ && (x[1] >= ymin)
+ && (x[1] <= ymax)
+ && (x[2] >= zmin)
+ && (x[2] <= zmax))''')*(im_dim==3))*(not im_is_cone)+('''
+
+ d1 = (x[0]-O[0])*n1[0]
+ + (x[1]-O[1])*n1[1]
+ + (x[2]-O[2])*n1[2];
+ d2 = (x[0]-O[0])*n2[0]
+ + (x[1]-O[1])*n2[1]
+ + (x[2]-O[2])*n2[2];
+ d3 = (x[0]-O[0])*n3[0]
+ + (x[1]-O[1])*n3[1]
+ + (x[2]-O[2])*n3[2];
+ d4 = (x[0]-O[0])*n4[0]
+ + (x[1]-O[1])*n4[1]
+ + (x[2]-O[2])*n4[2];
+
+ if ((d1 >= 0.)
+ && (d2 >= 0.)
+ && (d3 >= 0.)
+ && (d4 >= 0.))''')*(im_is_cone)+'''
{
expr[0] = 1.;
}
@@ -279,5 +396,5 @@ public:
};
}'''
- #print ExprCharFuncIm_cpp
+ # print ExprCharFuncIm_cpp
return ExprCharFuncIm_cpp
diff --git a/image_expressions_py.py b/image_expressions_py.py
index 7e7585bc81d09ab4128da02d85643cfc37fee449..dabb8f138c3cc5a6ba4b0f1cd839cfaa54c4d89b 100644
--- a/image_expressions_py.py
+++ b/image_expressions_py.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/plot_binned_strains_vs_radius.py b/plot_binned_strains_vs_radius.py
index 25bfc6fce9bc2c39451dd575ef8f0e7315fc5edd..8d0033e959bb4783d2e2e2216af647a8e770917b 100644
--- a/plot_binned_strains_vs_radius.py
+++ b/plot_binned_strains_vs_radius.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/plot_displacement_error.py b/plot_displacement_error.py
index 727e80dc5f294f19f764cf04604f0af4a9d4d017..57a575678a2f47592b05baed7e082325e8bed488 100644
--- a/plot_displacement_error.py
+++ b/plot_displacement_error.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/plot_regional_strains.py b/plot_regional_strains.py
index 4237804bb5dc51375a0de51e06b3db3b86824d54..ab3472bfeb05f4dc27ad81d0e58c0075dffa9368 100644
--- a/plot_regional_strains.py
+++ b/plot_regional_strains.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/plot_strains.py b/plot_strains.py
index 90872397ca9841c7a74ca58e27bd9f777001b90c..ffa967bed80b6d20056eaa742c987f587529a6df 100644
--- a/plot_strains.py
+++ b/plot_strains.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/plot_strains_vs_radius.py b/plot_strains_vs_radius.py
index 73b4e8d8f6d3fc07d5ca45a1ad24d510546b5f6d..8294398de7e47abaf53d647c77817fa6cf015cc6 100644
--- a/plot_strains_vs_radius.py
+++ b/plot_strains_vs_radius.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/project.py b/project.py
new file mode 100644
index 0000000000000000000000000000000000000000..62d518487c7ba588c8d5de39248f19f9a9359ea4
--- /dev/null
+++ b/project.py
@@ -0,0 +1,121 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Ezgi BerberoÄŸlu, 2018-2019 ###
+### ###
+### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
+### ###
+### And Martin Genet, 2016-2019 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import dolfin
+import numpy
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+def project(
+ mesh,
+ image_filename,
+ image_field_name="displacement",
+ image_quadrature=1):
+
+ dV = dolfin.Measure("dx", domain=mesh)
+ form_compiler_parameters_for_images = {}
+ form_compiler_parameters_for_images["quadrature_degree"] = image_quadrature
+
+ vfe = dolfin.VectorElement(
+ family="Quadrature",
+ cell=mesh.ufl_cell(),
+ degree=image_quadrature,
+ quad_scheme="default")
+
+ vfs = dolfin.VectorFunctionSpace(
+ mesh=mesh,
+ family="Lagrange",
+ degree=1)
+
+ projected_func = dolfin.Function(
+ vfs,
+ name=image_field_name)
+ U = dolfin.TrialFunction(
+ vfs)
+ V = dolfin.TestFunction(
+ vfs)
+
+ source_expr = dolfin.Expression(
+ cppcode='''\
+#include <vtkSmartPointer.h>
+#include <vtkStructuredGridReader.h>
+#include <vtkProbeFilter.h>
+#include <vtkStructuredGrid.h>
+#include <vtkPolyData.h>
+#include <vtkPointData.h>
+
+namespace dolfin
+{
+
+class MyExpr : public Expression
+{
+ vtkSmartPointer<vtkStructuredGrid> sgrid;
+ vtkSmartPointer<vtkPoints> probe_points;
+ vtkSmartPointer<vtkPolyData> probe_polydata;
+ vtkSmartPointer<vtkProbeFilter> probe_filter;
+
+public:
+
+ MyExpr():
+ Expression(3)
+ {
+ sgrid = vtkSmartPointer<vtkStructuredGrid>::New();
+ probe_points = vtkSmartPointer<vtkPoints>::New();
+ probe_polydata = vtkSmartPointer<vtkPolyData>::New();
+ probe_filter = vtkSmartPointer<vtkProbeFilter>::New();
+ }
+
+ void init_image(
+ const char* image_filename)
+ {
+ vtkSmartPointer<vtkStructuredGridReader> reader = vtkSmartPointer<vtkStructuredGridReader>::New();
+ reader->SetFileName(image_filename);
+ reader->Update();
+ sgrid = reader->GetOutput();
+ probe_filter->SetSourceData(sgrid);
+}
+
+ void eval(
+ Array<double>& expr,
+ const Array<double>& X) const
+ {
+ probe_points->SetNumberOfPoints(1);
+ probe_points->SetPoint(0, X.data());
+ probe_polydata->SetPoints(probe_points);
+ probe_filter->SetInputData(probe_polydata);
+ probe_filter->Update();
+ probe_filter->GetOutput()->GetPointData()->GetArray("'''+image_field_name+'''")->GetTuple(0, expr.data());
+ }
+};
+
+}''',
+ element=vfe)
+ source_expr.init_image(
+ image_filename)
+
+ M = dolfin.assemble(
+ dolfin.inner(U,V)*dV)
+
+ N = dolfin.assemble(
+ dolfin.inner(source_expr,V)*dV,
+ form_compiler_parameters=form_compiler_parameters_for_images)
+
+ dolfin.solve(M, projected_func.vector(), N)
+
+ return projected_func
diff --git a/write_VTU_file.py b/write_VTU_file.py
index 9a71f67414a5b5c34e6861c7171a0fc721998709..512d1b46e966665464213e2d098a789eb5e30c1e 100644
--- a/write_VTU_file.py
+++ b/write_VTU_file.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###