dolfin_dic now works with both fenics 2017 & 2019

Energy_WarpedImage.py

@@ -97,35 +97,59 @@ class WarpedImageEnergy(Energy):
         self.printer.print_var("ref_frame",self.ref_frame)
 
         # Iref
-        self.Iref = dolfin.Expression(
-            ddic.get_ExprIm_cpp(
+        if (int(dolfin.__version__.split('.')[0]) >= 2018):
+            name, cpp = ddic.get_ExprIm_cpp_pybind(
                 im_dim=self.image_series.dimension,
                 im_type="im",
-                im_is_def=0),
-            element=self.fe)
+                im_is_def=0)
+            module = dolfin.compile_cpp_code(cpp)
+            expr = getattr(module, name)
+            self.Iref = dolfin.CompiledExpression(
+                expr(),
+                element=self.fe)
+        else:
+            self.Iref = dolfin.Expression(
+                ddic.get_ExprIm_cpp_swig(
+                    im_dim=self.image_series.dimension,
+                    im_type="im",
+                    im_is_def=0),
+                element=self.fe)
         self.ref_image_filename = self.image_series.get_image_filename(self.ref_frame)
         self.Iref.init_image(self.ref_image_filename)
 
         self.Iref_int = dolfin.assemble(self.Iref * self.dV)/self.problem.mesh_V0
-        self.printer.print_var("Iref_int",self.Iref_int)
+        self.printer.print_sci("Iref_int",self.Iref_int)
 
         self.Iref_norm = (dolfin.assemble(self.Iref**2 * self.dV)/self.problem.mesh_V0)**(1./2)
         assert (self.Iref_norm > 0.),\
             "Iref_norm = "+str(self.Iref_norm)+" <= 0. Aborting."
-        self.printer.print_var("Iref_norm",self.Iref_norm)
+        self.printer.print_sci("Iref_norm",self.Iref_norm)
 
         # DIref
-        self.DIref = dolfin.Expression(
-            ddic.get_ExprIm_cpp(
+        if (int(dolfin.__version__.split('.')[0]) >= 2018):
+            name, cpp = ddic.get_ExprIm_cpp_pybind(
                 im_dim=self.image_series.dimension,
                 im_type="grad" if (self.image_series.grad_basename is None) else "grad_no_deriv",
-                im_is_def=0),
-            element=self.ve)
+                im_is_def=0)
+            module = dolfin.compile_cpp_code(cpp)
+            expr = getattr(module, name)
+            self.DIref = dolfin.CompiledExpression(
+                expr(),
+                element=self.ve)
+        else:
+            cpp = ddic.get_ExprIm_cpp_swig(
+                im_dim=self.image_series.dimension,
+                im_type="grad" if (self.image_series.grad_basename is None) else "grad_no_deriv",
+                im_is_def=0)
+            self.DIref = dolfin.Expression(
+                cppcode=cpp,
+                element=self.ve)
         self.ref_image_grad_filename = self.image_series.get_image_grad_filename(self.ref_frame)
         self.DIref.init_image(self.ref_image_grad_filename)
 
         self.printer.dec()
         self.printer.print_str("Defining deformed image…")
+        self.printer.inc()
 
         self.scaling = numpy.array([1.,0.])
         if (self.dynamic_scaling):
@@ -133,82 +157,172 @@ class WarpedImageEnergy(Energy):
             self.q = numpy.empty(2)
 
         # Idef
-        self.Idef = dolfin.Expression(
-            ddic.get_ExprIm_cpp(
+        if (int(dolfin.__version__.split('.')[0]) >= 2018):
+            name, cpp = ddic.get_ExprIm_cpp_pybind(
                 im_dim=self.image_series.dimension,
                 im_type="im",
-                im_is_def=1),
-            element=self.fe)
+                im_is_def=1)
+            module = dolfin.compile_cpp_code(cpp)
+            expr = getattr(module, name)
+            self.Idef = dolfin.CompiledExpression(
+                expr(),
+                element=self.fe)
+            self.Idef.init_disp(self.problem.U.cpp_object())
+        else:
+            cpp = ddic.get_ExprIm_cpp_swig(
+                im_dim=self.image_series.dimension,
+                im_type="im",
+                im_is_def=1)
+            self.Idef = dolfin.Expression(
+                cppcode=cpp,
+                element=self.fe)
+            self.Idef.init_disp(self.problem.U)
         self.Idef.init_image(self.ref_image_filename)
-        self.Idef.init_disp(self.problem.U)
         self.Idef.init_dynamic_scaling(self.scaling)
 
+        self.Idef_int = dolfin.assemble(self.Idef * self.dV)/self.problem.mesh_V0
+        self.printer.print_sci("Idef_int",self.Idef_int)
+
         # DIdef
-        self.DIdef = dolfin.Expression(
-            ddic.get_ExprIm_cpp(
+        if (int(dolfin.__version__.split('.')[0]) >= 2018):
+            name, cpp = ddic.get_ExprIm_cpp_pybind(
                 im_dim=self.image_series.dimension,
                 im_type="grad" if (self.image_series.grad_basename is None) else "grad_no_deriv",
-                im_is_def=1),
-            element=self.ve)
+                im_is_def=1)
+            module = dolfin.compile_cpp_code(cpp)
+            expr = getattr(module, name)
+            self.DIdef = dolfin.CompiledExpression(
+                expr(),
+                element=self.ve)
+            self.DIdef.init_disp(self.problem.U.cpp_object())
+        else:
+            cpp = ddic.get_ExprIm_cpp_swig(
+                im_dim=self.image_series.dimension,
+                im_type="grad" if (self.image_series.grad_basename is None) else "grad_no_deriv",
+                im_is_def=1)
+            self.DIdef = dolfin.Expression(
+                cppcode=cpp,
+                element=self.ve)
+            self.DIdef.init_disp(self.problem.U)
         self.DIdef.init_image(self.ref_image_filename)
-        self.DIdef.init_disp(self.problem.U)
         self.DIdef.init_dynamic_scaling(self.scaling)
 
+        self.printer.dec()
         self.printer.print_str("Defining previous image…")
+        self.printer.inc()
 
         # Iold
-        self.Iold = dolfin.Expression(
-            ddic.get_ExprIm_cpp(
+        if (int(dolfin.__version__.split('.')[0]) >= 2018):
+            name, cpp = ddic.get_ExprIm_cpp_pybind(
+                im_dim=self.image_series.dimension,
+                im_type="im",
+                im_is_def=1)
+            module = dolfin.compile_cpp_code(cpp)
+            expr = getattr(module, name)
+            self.Iold = dolfin.CompiledExpression(
+                expr(),
+                element=self.fe)
+            self.Iold.init_disp(self.problem.Uold.cpp_object())
+        else:
+            cpp = ddic.get_ExprIm_cpp_swig(
                 im_dim=self.image_series.dimension,
                 im_type="im",
-                im_is_def=1),
-            element=self.fe)
+                im_is_def=1)
+            self.Iold = dolfin.Expression(
+                cppcode=cpp,
+                element=self.fe)
+            self.Iold.init_disp(self.problem.Uold)
         self.Iold.init_image(self.ref_image_filename)
-        self.Iold.init_disp(self.problem.Uold)
         self.Iold.init_dynamic_scaling(self.scaling) # 2016/07/25: ok, same scaling must apply to Idef & Iold…
 
         # DIold
-        self.DIold = dolfin.Expression(
-            ddic.get_ExprIm_cpp(
+        if (int(dolfin.__version__.split('.')[0]) >= 2018):
+            name, cpp = ddic.get_ExprIm_cpp_pybind(
+                im_dim=self.image_series.dimension,
+                im_type="grad" if (self.image_series.grad_basename is None) else "grad_no_deriv",
+                im_is_def=1)
+            module = dolfin.compile_cpp_code(cpp)
+            expr = getattr(module, name)
+            self.DIold = dolfin.CompiledExpression(
+                expr(),
+                element=self.ve)
+            self.DIold.init_disp(self.problem.Uold.cpp_object())
+        else:
+            cpp = ddic.get_ExprIm_cpp_swig(
                 im_dim=self.image_series.dimension,
                 im_type="grad" if (self.image_series.grad_basename is None) else "grad_no_deriv",
-                im_is_def=1),
-            element=self.ve)
+                im_is_def=1)
+            self.DIold = dolfin.Expression(
+                cppcode=cpp,
+                element=self.ve)
+            self.DIold.init_disp(self.problem.Uold)
         self.DIold.init_image(self.ref_image_filename)
-        self.DIold.init_disp(self.problem.Uold)
         self.DIold.init_dynamic_scaling(self.scaling) # 2016/07/25: ok, same scaling must apply to Idef & Iold…
 
+        self.printer.dec()
         self.printer.print_str("Defining correlation energy…")
+        self.printer.inc()
 
         # Phi_ref
-        self.Phi_Iref = dolfin.Expression(
-            ddic.get_ExprCharFuncIm_cpp(
+        if (int(dolfin.__version__.split('.')[0]) >= 2018):
+            name, cpp = ddic.get_ExprCharFuncIm_cpp_pybind(
                 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)
+                im_is_cone=im_is_cone)
+            module = dolfin.compile_cpp_code(cpp)
+            expr = getattr(module, name)
+            self.Phi_ref = dolfin.CompiledExpression(
+                expr(),
+                element=self.fe)
+        else:
+            cpp = ddic.get_ExprCharFuncIm_cpp_swig(
+                im_dim=self.image_series.dimension,
+                im_is_def=0,
+                im_is_cone=im_is_cone)
+            self.Phi_ref = dolfin.Expression(
+                cppcode=cpp,
+                element=self.fe)
+        self.Phi_ref.init_image(self.ref_image_filename)
+
+        self.Phi_ref_int = dolfin.assemble(self.Phi_ref * self.dV)/self.problem.mesh_V0
+        self.printer.print_sci("Phi_ref_int",self.Phi_ref_int)
 
         # Phi_def
-        self.Phi_Idef = dolfin.Expression(
-            ddic.get_ExprCharFuncIm_cpp(
+        if (int(dolfin.__version__.split('.')[0]) >= 2018):
+            name, cpp = ddic.get_ExprCharFuncIm_cpp_pybind(
                 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)
+                im_is_cone=im_is_cone)
+            module = dolfin.compile_cpp_code(cpp)
+            expr = getattr(module, name)
+            self.Phi_def = dolfin.CompiledExpression(
+                expr(),
+                element=self.fe)
+            self.Phi_def.init_disp(self.problem.U.cpp_object())
+        else:
+            cpp = ddic.get_ExprCharFuncIm_cpp_swig(
+                im_dim=self.image_series.dimension,
+                im_is_def=1,
+                im_is_cone=im_is_cone)
+            self.Phi_def = dolfin.Expression(
+                cppcode=cpp,
+                element=self.fe)
+            self.Phi_def.init_disp(self.problem.U)
+        self.Phi_def.init_image(self.ref_image_filename)
+
+        self.Phi_def_int = dolfin.assemble(self.Phi_def * self.dV)/self.problem.mesh_V0
+        self.printer.print_sci("Phi_def_int",self.Phi_def_int)
 
         # Psi_c
-        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.Psi_c  = self.Phi_def * self.Phi_ref * (self.Idef - self.Iref)**2/2
+        self.DPsi_c = self.Phi_def * self.Phi_ref * (self.Idef - self.Iref) * dolfin.dot(self.DIdef, 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.DDPsi_c     = self.Phi_def * self.Phi_ref * dolfin.dot(self.DIdef, self.problem.dU_trial) * dolfin.dot(self.DIdef, self.problem.dU_test)
+        self.DDPsi_c_old = self.Phi_def * self.Phi_ref * dolfin.dot(self.DIold, self.problem.dU_trial) * dolfin.dot(self.DIold, self.problem.dU_test)
+        self.DDPsi_c_ref = self.Phi_def * self.Phi_ref * dolfin.dot(self.DIref, self.problem.dU_trial) * dolfin.dot(self.DIref, 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)
+        self.Psi_c_old   = self.Phi_def * self.Phi_ref * (self.Idef - self.Iold)**2/2
+        self.DPsi_c_old  = self.Phi_def * self.Phi_ref * (self.Idef - self.Iold) * dolfin.dot(self.DIdef, self.problem.dU_test)
 
         # forms
         self.ener_form = self.Psi_c   * self.dV
@@ -268,11 +382,12 @@ class WarpedImageEnergy(Energy):
             self.scaling[:] = numpy.linalg.solve(self.p, self.q)
             self.printer.print_var("scaling",self.scaling)
 
-            self.Idef.update_dynamic_scaling(self.scaling)  # should not be needed
-            self.DIdef.update_dynamic_scaling(self.scaling) # should not be needed
+            if (int(dolfin.__version__.split('.')[0]) <= 2017):
+                self.Idef.update_dynamic_scaling(self.scaling)  # should not be needed
+                self.DIdef.update_dynamic_scaling(self.scaling) # should not be needed
 
-            self.Iold.update_dynamic_scaling(self.scaling)  # should not be needed
-            self.DIold.update_dynamic_scaling(self.scaling) # should not be needed
+                self.Iold.update_dynamic_scaling(self.scaling)  # should not be needed
+                self.DIold.update_dynamic_scaling(self.scaling) # should not be needed
 
             self.get_qoi_values()
 

ImageIterator.py

@@ -67,6 +67,7 @@ class ImageIterator():
             filename=qoi_filebasename+".dat")
 
         if not (self.register_ref_frame):
+            self.printer.dec()
             self.printer.print_str("Writing initial solution…")
             self.printer.inc()
 
@@ -83,8 +84,7 @@ class ImageIterator():
             qoi_printer.write_line(
                 values=qoi_values)
 
-            self.printer.dec()
-
+        self.printer.dec()
         self.printer.print_str("Looping over frames…")
 
         if (self.initialize_U_from_file):

ImageSeries.py

@@ -9,7 +9,6 @@
 ################################################################################
 
 # from builtins     import *
-# from future.utils import native_str
 
 import dolfin
 import glob
@@ -84,7 +83,6 @@ class ImageSeries():
     def get_image_filename(self,
             k_frame):
 
-        # return native_str(self.folder+"/"+self.basename+"_"+str(k_frame).zfill(self.zfill)+"."+self.ext)
         return self.folder+"/"+self.basename+"_"+str(k_frame).zfill(self.zfill)+"."+self.ext
 
 
@@ -92,10 +90,6 @@ class ImageSeries():
     def get_image_grad_filename(self,
             k_frame):
 
-        # if (self.grad_basename is None):
-        #     return native_str(self.get_image_filename(k_frame))
-        # else:
-        #     return native_str(self.grad_folder+"/"+self.grad_basename+"_"+str(k_frame).zfill(self.zfill)+"."+self.ext)
         if (self.grad_basename is None):
             return self.get_image_filename(k_frame)
         else:

NonlinearSolver.py

@@ -45,8 +45,8 @@ class NonlinearSolver():
             self.jac_mat,
             self.linear_solver_name)
         self.linear_solver.parameters['report']               = bool(0)
-        self.linear_solver.parameters['reuse_factorization']  = bool(0)
-        self.linear_solver.parameters['same_nonzero_pattern'] = bool(1)
+        # self.linear_solver.parameters['reuse_factorization']  = bool(0)
+        # self.linear_solver.parameters['same_nonzero_pattern'] = bool(1)
         self.linear_solver.parameters['symmetric']            = bool(1)
         self.linear_solver.parameters['verbose']              = bool(0)
 
@@ -316,14 +316,14 @@ class NonlinearSolver():
             if (relax_k > 1):
                 ener_min_old = ener_min
             relax_min = relax_list[numpy.argmin(ener_list)]
-            # self.printer.print_var("relax_min",relax_min)
+            # self.printer.print_sci("relax_min",relax_min)
             ener_min = min(ener_list)
-            # self.printer.print_var("ener_min",ener_min)
+            # self.printer.print_sci("ener_min",ener_min)
             if (relax_k > 1):
                 dener_min = ener_min-ener_min_old
-                self.printer.print_var("dener_min",dener_min)
+                self.printer.print_sci("dener_min",dener_min)
                 relax_err = dener_min/ener_list[0]
-                self.printer.print_var("relax_err",relax_err)
+                self.printer.print_sci("relax_err",relax_err)
                 if (relax_min != 0.) and (abs(relax_err) < self.relax_tol):
                     break
             if (relax_k == self.relax_n_iter_max):

Problem_ImageRegistration.py

@@ -196,9 +196,9 @@ class ImageRegistrationProblem(Problem):
         for energy in self.energies:
             ener_ = dolfin.assemble(
                 energy.ener_form)
-            self.printer.print_var("ener_"+energy.name,ener_)
+            self.printer.print_sci("ener_"+energy.name,ener_)
             ener += energy.w * ener_
-        #self.printer.print_var("ener",ener)
+        #self.printer.print_sci("ener",ener)
 
         # ener_form = 0.
         # for energy in self.energies:
@@ -206,7 +206,7 @@ class ImageRegistrationProblem(Problem):
         #
         # ener = dolfin.assemble(
         #     ener_form)
-        # #self.printer.print_var("ener",ener)
+        # #self.printer.print_sci("ener",ener)
 
         return ener
 

__init__.py

@@ -30,6 +30,7 @@ from .generate_undersampled_images import *
 from .plot_regional_strains import *
 from .compute_displacements_from_ref import *
 from .ImageIterator import *
+from .expressions_char_func_cpp import *
 from .fedic2 import *
 from .plot_displacement_error import *
 from .expressions_generated_images_cpp import *

expressions_char_func_cpp.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2019                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import dolfin
+
+import dolfin_dic as ddic
+
+################################################################################
+
+def get_ExprCharFuncIm_cpp_pybind(
+        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)
+
+    name  = "Expr"
+    name += str(im_dim)
+    name += "CharFuncIm"
+    if   (im_is_def == 0):
+        name += "Ref"
+    elif (im_is_def == 1):
+        name += "Def"
+    # print(name)
+
+    if   (im_is_def == 0):
+        X_or_x = "X"
+    elif (im_is_def == 1):
+        X_or_x = "x"
+
+    cpp = '''\
+#include <math.h>
+#include <string.h>
+
+#include <dolfin/function/Expression.h>
+#include <dolfin/function/Function.h>
+
+#include <vtkSmartPointer.h>
+#include <vtkXMLImageDataReader.h>
+#include <vtkImageData.h>
+
+#include <pybind11/eigen.h>
+#include <pybind11/pybind11.h>
+
+class '''+name+''' : public dolfin::Expression
+{
+public:
+
+    double xmin, xmax, ymin, ymax, zmin, zmax;'''+('''
+
+    mutable Eigen::Vector'''+str(im_dim)+'''d UX, x;
+    std::shared_ptr<dolfin::Function> U;''')*(im_is_def)+('''
+
+    Eigen::Vector3d O, n1, n2, n3, n4;
+    mutable double d1, d2, d3, d4;''')*(im_is_cone)+'''
+
+    '''+name+'''() : dolfin::Expression() {}
+
+    void init_image(
+        const char* filename) {
+
+        vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New();
+        reader->SetFileName(filename);
+        reader->Update();
+
+        vtkSmartPointer<vtkImageData> image = reader->GetOutput();
+        double* bounds = image->GetBounds();
+        xmin = bounds[0];
+        xmax = bounds[1];
+        ymin = bounds[2];
+        ymax = bounds[3];
+        zmin = bounds[4];
+        zmax = bounds[5];'''+('''
+        std::cout << "bounds = " << bounds[0] << " " << bounds[1] << " " << bounds[2] << " " << bounds[3] << " " << bounds[4] << " " << bounds[5] << std::endl;
+        std::cout << "xmin = " << xmin << std::endl;
+        std::cout << "xmax = " << xmax << std::endl;
+        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)+('''
+
+        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 init_disp(
+        std::shared_ptr<dolfin::Function> U_) {
+
+        U = U_;}''')*(im_is_def)+'''
+
+    void eval(
+        Eigen::Ref<      Eigen::VectorXd> expr,
+        Eigen::Ref<const Eigen::VectorXd> X   ) const {'''+('''
+
+        std::cout << "X = " << X << std::endl;''')*(verbose)+('''
+
+        U->eval(UX, X);'''+('''
+        std::cout << "UX = " << UX << 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 << std::endl;''')*(verbose))*(im_is_def)+('''
+
+        d1 = ('''+X_or_x+'''[0]-O[0])*n1[0]
+           + ('''+X_or_x+'''[1]-O[1])*n1[1]
+           + ('''+X_or_x+'''[2]-O[2])*n1[2];
+        d2 = ('''+X_or_x+'''[0]-O[0])*n2[0]
+           + ('''+X_or_x+'''[1]-O[1])*n2[1]
+           + ('''+X_or_x+'''[2]-O[2])*n2[2];
+        d3 = ('''+X_or_x+'''[0]-O[0])*n3[0]
+           + ('''+X_or_x+'''[1]-O[1])*n3[1]
+           + ('''+X_or_x+'''[2]-O[2])*n3[2];
+        d4 = ('''+X_or_x+'''[0]-O[0])*n4[0]
+           + ('''+X_or_x+'''[1]-O[1])*n4[1]
+           + ('''+X_or_x+'''[2]-O[2])*n4[2];''')*(im_is_cone)+('''
+
+        if (('''+X_or_x+'''[0] >= xmin)
+         && ('''+X_or_x+'''[0] <= xmax)
+         && ('''+X_or_x+'''[1] >= ymin)
+         && ('''+X_or_x+'''[1] <= ymax))''')*(im_dim==2)+('''
+        if (('''+X_or_x+'''[0] >= xmin)
+         && ('''+X_or_x+'''[0] <= xmax)
+         && ('''+X_or_x+'''[1] >= ymin)
+         && ('''+X_or_x+'''[1] <= ymax)
+         && ('''+X_or_x+'''[2] >= zmin)
+         && ('''+X_or_x+'''[2] <= zmax)'''+('''
+         && (d1 >= 0.)
+         && (d2 >= 0.)
+         && (d3 >= 0.)
+         && (d4 >= 0.)''')*(im_is_cone)+''')''')*(im_dim==3)+''' {
+            expr[0] = 1.;}
+        else {
+            expr[0] = 0.;}'''+('''
+        std::cout << "expr = " << expr << std::endl;''')*(verbose)+'''}
+};
+
+PYBIND11_MODULE(SIGNATURE, m)
+{
+    pybind11::class_<'''+name+''', std::shared_ptr<'''+name+'''>, dolfin::Expression>
+    (m, "'''+name+'''")
+    .def(pybind11::init<>())
+    .def("init_image", &'''+name+'''::init_image, pybind11::arg("filename"))'''+('''
+    .def("init_disp", &'''+name+'''::init_disp, pybind11::arg("U_"))''')*(im_is_def)+''';
+}
+'''
+    # print(cpp)
+
+    return name, cpp
+
+################################################################################
+
+def get_ExprCharFuncIm_cpp_swig(
+        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)
+
+    if   (im_is_def == 0):
+        X_or_x = "X"
+    elif (im_is_def == 1):
+        X_or_x = "x"
+
+    ExprCharFuncIm_cpp = '''\
+#include <math.h>
+#include <string.h>
+
+#include <vtkSmartPointer.h>
+#include <vtkXMLImageDataReader.h>
+#include <vtkImageData.h>
+
+namespace dolfin
+{
+
+class MyExpr : public Expression
+{
+
+    double xmin, xmax, ymin, ymax, zmin, zmax;'''+('''
+
+    mutable Array<double> UX, x;
+    std::shared_ptr<Function> U;''')*(im_is_def)+('''
+
+    Array<double> O, n1, n2, n3, n4;
+    mutable double d1, d2, d3, d4;''')*(im_is_cone)+'''
+
+public:
+
+    MyExpr() :
+        Expression()'''+(''',
+        UX('''+str(im_dim)+'''),
+        x('''+str(im_dim)+''')''')*(im_is_def)+(''',
+        O(3),
+        n1(3),
+        n2(3),
+        n3(3),
+        n4(3)''')*(im_is_cone)+''' {}
+
+    void init_image(
+        const char* filename) {
+
+        vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New();
+        reader->SetFileName(filename);
+        reader->Update();
+
+        vtkSmartPointer<vtkImageData> image = reader->GetOutput();
+        double* bounds = image->GetBounds();
+        xmin = bounds[0];
+        xmax = bounds[1];
+        ymin = bounds[2];
+        ymax = bounds[3];
+        zmin = bounds[4];
+        zmax = bounds[5];'''+('''
+        std::cout << "bounds = " << bounds[0] << " " << bounds[1] << " " << bounds[2] << " " << bounds[3] << " " << bounds[4] << " " << bounds[5] << std::endl;
+        std::cout << "xmin = " << xmin << std::endl;
+        std::cout << "xmax = " << xmax << std::endl;
+        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)+('''
+
+        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 init_disp(
+        std::shared_ptr<Function> U_) {
+
+        U = U_;}''')*(im_is_def)+'''
+
+    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))*(im_is_def)+('''
+
+       d1 = ('''+X_or_x+'''[0]-O[0])*n1[0]
+          + ('''+X_or_x+'''[1]-O[1])*n1[1]
+          + ('''+X_or_x+'''[2]-O[2])*n1[2];
+       d2 = ('''+X_or_x+'''[0]-O[0])*n2[0]
+          + ('''+X_or_x+'''[1]-O[1])*n2[1]
+          + ('''+X_or_x+'''[2]-O[2])*n2[2];
+       d3 = ('''+X_or_x+'''[0]-O[0])*n3[0]
+          + ('''+X_or_x+'''[1]-O[1])*n3[1]
+          + ('''+X_or_x+'''[2]-O[2])*n3[2];
+       d4 = ('''+X_or_x+'''[0]-O[0])*n4[0]
+          + ('''+X_or_x+'''[1]-O[1])*n4[1]
+          + ('''+X_or_x+'''[2]-O[2])*n4[2];''')*(im_is_cone)+('''
+
+        if (('''+X_or_x+'''[0] >= xmin)
+         && ('''+X_or_x+'''[0] <= xmax)
+         && ('''+X_or_x+'''[1] >= ymin)
+         && ('''+X_or_x+'''[1] <= ymax))''')*(im_dim==2)+('''
+        if (('''+X_or_x+'''[0] >= xmin)
+         && ('''+X_or_x+'''[0] <= xmax)
+         && ('''+X_or_x+'''[1] >= ymin)
+         && ('''+X_or_x+'''[1] <= ymax)
+         && ('''+X_or_x+'''[2] >= zmin)
+         && ('''+X_or_x+'''[2] <= zmax)'''+('''
+         && (d1 >= 0.)
+         && (d2 >= 0.)
+         && (d3 >= 0.)
+         && (d4 >= 0.)''')*(im_is_cone)+''')''')*(im_dim==3)+''' {
+            expr[0] = 1.;}
+        else {
+            expr[0] = 0.;}'''+('''
+        std::cout << "expr = " << expr.str(1) << std::endl;''')*(verbose)+'''}
+};
+
+}
+'''
+    # print(ExprCharFuncIm_cpp)
+
+    return ExprCharFuncIm_cpp

expressions_static_scaling_cpp.py

@@ -14,14 +14,19 @@ def get_StaticScaling_cpp():
 
     return '''\
 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);
-}
+    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);}
 '''