diff --git a/fedic.py b/fedic.py
index c539797c7fd4c8680b35e5a1ce0e2684d0c0c81f..6c46a14c0ee7d3ab3478897bdda9856325d9b333 100644
--- a/fedic.py
+++ b/fedic.py
@@ -330,7 +330,7 @@ def fedic(
Div_P)
N = dolfin.FacetNormal(mesh)
Jump_P_N = dolfin.jump(P_m, N)
- h = mesh.hmin()
+ h = dolfin.Constant(mesh.hmin())
psi_m_F = dolfin.dot(Jump_P_N,
Jump_P_N)/h
#P_N = P_m * N
@@ -456,9 +456,14 @@ def fedic(
assert (0), "\"images_expressions_type\" (="+str(images_expressions_type)+") must be \"cpp\" or \"py\". Aborting."
mypy.print_str("Defining correlation energy…",tab)
- psi_c = (Idef - Iref)**2/2
- Dpsi_c = (Idef - Iref) * dolfin.dot(DIdef, dU_test)
- DDpsi_c = dolfin.dot(DIdef, dU_trial) * dolfin.dot(DIdef, dU_test)
+ phi_Iref = dolfin.Expression(
+ cppcode=ddic.get_ExprCharFuncIm_cpp(
+ im_dim=images_dimension),
+ element=fe)
+ phi_Iref.init_image(ref_image_filename)
+ psi_c = phi_Iref * (Idef - Iref)**2/2
+ Dpsi_c = phi_Iref * (Idef - Iref) * dolfin.dot(DIdef, dU_test)
+ DDpsi_c = phi_Iref * dolfin.dot(DIdef, dU_trial) * dolfin.dot(DIdef, dU_test)
if ("-wHess" in tangent_type):
DDpsi_c += (Idef - Iref) * dolfin.inner(dolfin.dot(DDIdef, dU_trial), dU_test)
@@ -831,7 +836,7 @@ def fedic(
if (print_iterations):
os.remove(frame_basename+"_.pvd")
file_dat_frame.close()
- os.system("gnuplot -e \"set terminal pdf; set output '"+frame_basename+".pdf'; set key box textcolor variable; set grid; set logscale y; set yrange [1e-3:1e0]; plot '"+frame_basename+".dat' u 1:3 pt 1 lw 3 title 'err_res', '' u 1:7 pt 1 lw 3 title 'err_dU', '' using 1:9 pt 1 lw 3 title 'err_im', "+str(tol_res or tol_dU or tol_im)+" lt -1 notitle; unset logscale y; set yrange [*:*]; plot '' u 1:4 pt 1 lw 3 title 'relax'\"")
+ os.system("gnuplot -e \"set terminal pdf; set output '"+frame_basename+".pdf'; set key box textcolor variable; set grid; set logscale y; set yrange [1e-3:1e0]; plot '"+frame_basename+".dat' u 1:3 pt 1 lw 3 title 'err_res', '' u 1:8 pt 1 lw 3 title 'err_dU', '' using 1:10 pt 1 lw 3 title 'err_im', "+str(tol_res or tol_dU or tol_im)+" lt -1 notitle; unset logscale y; set yrange [*:*]; plot '' u 1:4 pt 1 lw 3 title 'relax'\"")
if not (success) and not (continue_after_fail):
break
diff --git a/image_expressions_cpp.py b/image_expressions_cpp.py
index 4a801e4f469bde320eb52a83b73dc31b0e76f661..9d31940f86f5ae38f5e1d888e1bf3d2e8b98e477 100644
--- a/image_expressions_cpp.py
+++ b/image_expressions_cpp.py
@@ -58,8 +58,8 @@ public:
MyExpr():
Expression('''+str(im_dim)*(im_type=="grad")+''')'''+(''',
- dynamic_scaling_a(1.),
- dynamic_scaling_b(0.),
+ dynamic_scaling_a(1.), // should not be needed
+ dynamic_scaling_b(0.), // should not be needed
UX('''+str(im_dim)+'''),
x(3)''')*(im_is_def)+(''',
X3D(3)''')*(not im_is_def)*(im_dim==2)+'''
@@ -91,7 +91,7 @@ public:
void init_image(
const char* filename,
const char* interpol_mode="'''+('''linear''')*(im_type=="im")+('''linear''')*(im_type=="grad")+'''",
- const double &interpol_out_value='''+('''0.''')*(im_type=="im")+('''1.''')*(im_type=="grad")+(''',
+ const double &interpol_out_value='''+('''0.''')*(im_type=="im")+('''0.''')*(im_type=="grad")+(''',
const double &Z=0.''')*(im_dim==2)+''')
{
vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New();
@@ -206,3 +206,76 @@ public:
}'''
#print ExprIm_cpp
return ExprIm_cpp
+
+def get_ExprCharFuncIm_cpp(
+ im_dim,
+ verbose=0):
+
+ assert (im_dim in (2,3))
+
+ ExprCharFuncIm_cpp = '''\
+#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;
+
+public:
+
+ MyExpr():
+ 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)+'''
+ }
+
+ 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)+('''
+
+ if ((X[0] >= xmin) && (X[0] <= xmax) && (X[1] >= ymin) && (X[1] <= ymax) && (X[2] >= zmin) && (X[2] <= zmax))''')*(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
diff --git a/plot_strains.py b/plot_strains.py
index 00cf3cb2ad8ce3274dca19a34edc7cf1ea8805e3..b0d2cccc07a82971265f89e74cfeb9447a72f332 100644
--- a/plot_strains.py
+++ b/plot_strains.py
@@ -93,6 +93,7 @@ set ylabel "'''+comp_name+''' strain (%)"
if ("-" in comp_name):
plotfile.write('''\
set yrange [-10:10]
+
''')
else:
plotfile.write('''\
@@ -104,8 +105,8 @@ plot 0 linecolor rgb "black" notitle,\\
''')
if (ref_folder is not None) and (ref_basename is not None):
plotfile.write('''\
- "'''+ref_folder+'''/'''+ref_basename+'''-strains.dat" using ($1):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linecolor "black" linewidth 5 notitle,\
- "'''+ref_folder+'''/'''+ref_basename+'''-strains.dat" using ($1):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linecolor "black" pointtype 1 linewidth 1 notitle'''+(len(working_basenames)>0)*(''',\\
+ "'''+ref_folder+'''/'''+ref_basename+'''-strains.dat" using ($1):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linecolor "black" linewidth 5 notitle,\\
+ "'''+ref_folder+'''/'''+ref_basename+'''-strains.dat" using ($1):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linecolor "black" pointtype 1 linewidth 1 notitle'''+(len(working_basenames)>0)*(''',\\
''')+(len(working_basenames)==0)*('''
'''))
@@ -113,8 +114,8 @@ plot 0 linecolor rgb "black" notitle,\\
working_basename = working_basenames[k_basename]
working_scaling = working_scalings[k_basename]
plotfile.write('''\
- "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ('''+str(working_scaling)+'''*($1)+'''+str(k_basename)+'''./10):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linestyle '''+str(k_basename+1)+''' linewidth 5 title "'''+working_basename+'''",\
- "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ('''+str(working_scaling)+'''*($1)+'''+str(k_basename)+'''./10):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linestyle '''+str(k_basename+1)+''' linewidth 1 notitle'''+(k_basename<len(working_basenames)-1)*(''',\\
+ "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ('''+str(working_scaling)+'''*($1)+'''+str(k_basename)+'''./10):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linestyle '''+str(k_basename+1)+''' linewidth 5 title "'''+working_basename+'''",\\
+ "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ('''+str(working_scaling)+'''*($1)+'''+str(k_basename)+'''./10):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linestyle '''+str(k_basename+1)+''' linewidth 1 notitle'''+(k_basename<len(working_basenames)-1)*(''',\\
''')+(k_basename==len(working_basenames)-1)*('''
'''))