Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 6a06c102 authored by Martin Genet's avatar Martin Genet
Browse files

Adding image characteristic function for when the mesh extends outside the image

parent 5f045696
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
fedic.py
+ 10
− 5
View file @ 6a06c102
...@@ -330,7 +330,7 @@ def fedic( ...@@ -330,7 +330,7 @@ def fedic(
Div_P) Div_P)
N = dolfin.FacetNormal(mesh) N = dolfin.FacetNormal(mesh)
Jump_P_N = dolfin.jump(P_m, N) Jump_P_N = dolfin.jump(P_m, N)
h = mesh.hmin() h = dolfin.Constant(mesh.hmin())
psi_m_F = dolfin.dot(Jump_P_N, psi_m_F = dolfin.dot(Jump_P_N,
Jump_P_N)/h Jump_P_N)/h
#P_N = P_m * N #P_N = P_m * N
...@@ -456,9 +456,14 @@ def fedic( ...@@ -456,9 +456,14 @@ def fedic(
assert (0), "\"images_expressions_type\" (="+str(images_expressions_type)+") must be \"cpp\" or \"py\". Aborting." assert (0), "\"images_expressions_type\" (="+str(images_expressions_type)+") must be \"cpp\" or \"py\". Aborting."
mypy.print_str("Defining correlation energy…",tab) mypy.print_str("Defining correlation energy…",tab)
psi_c = (Idef - Iref)**2/2 phi_Iref = dolfin.Expression(
Dpsi_c = (Idef - Iref) * dolfin.dot(DIdef, dU_test) cppcode=ddic.get_ExprCharFuncIm_cpp(
DDpsi_c = dolfin.dot(DIdef, dU_trial) * dolfin.dot(DIdef, dU_test) 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): if ("-wHess" in tangent_type):
DDpsi_c += (Idef - Iref) * dolfin.inner(dolfin.dot(DDIdef, dU_trial), dU_test) DDpsi_c += (Idef - Iref) * dolfin.inner(dolfin.dot(DDIdef, dU_trial), dU_test)
...@@ -831,7 +836,7 @@ def fedic( ...@@ -831,7 +836,7 @@ def fedic(
if (print_iterations): if (print_iterations):
os.remove(frame_basename+"_.pvd") os.remove(frame_basename+"_.pvd")
file_dat_frame.close() 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): if not (success) and not (continue_after_fail):
break break
......
image_expressions_cpp.py
+ 76
− 3
View file @ 6a06c102
...@@ -58,8 +58,8 @@ public: ...@@ -58,8 +58,8 @@ public:
MyExpr(): MyExpr():
Expression('''+str(im_dim)*(im_type=="grad")+''')'''+(''', Expression('''+str(im_dim)*(im_type=="grad")+''')'''+(''',
dynamic_scaling_a(1.), dynamic_scaling_a(1.), // should not be needed
dynamic_scaling_b(0.), dynamic_scaling_b(0.), // should not be needed
UX('''+str(im_dim)+'''), UX('''+str(im_dim)+'''),
x(3)''')*(im_is_def)+(''', x(3)''')*(im_is_def)+(''',
X3D(3)''')*(not im_is_def)*(im_dim==2)+''' X3D(3)''')*(not im_is_def)*(im_dim==2)+'''
...@@ -91,7 +91,7 @@ public: ...@@ -91,7 +91,7 @@ public:
void init_image( void init_image(
const char* filename, const char* filename,
const char* interpol_mode="'''+('''linear''')*(im_type=="im")+('''linear''')*(im_type=="grad")+'''", 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)+''') const double &Z=0.''')*(im_dim==2)+''')
{ {
vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New(); vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New();
...@@ -206,3 +206,76 @@ public: ...@@ -206,3 +206,76 @@ public:
}''' }'''
#print ExprIm_cpp #print ExprIm_cpp
return 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
plot_strains.py
+ 5
− 4
View file @ 6a06c102
...@@ -93,6 +93,7 @@ set ylabel "'''+comp_name+''' strain (%)" ...@@ -93,6 +93,7 @@ set ylabel "'''+comp_name+''' strain (%)"
if ("-" in comp_name): if ("-" in comp_name):
plotfile.write('''\ plotfile.write('''\
set yrange [-10:10] set yrange [-10:10]
''') ''')
else: else:
plotfile.write('''\ plotfile.write('''\
...@@ -104,8 +105,8 @@ plot 0 linecolor rgb "black" notitle,\\ ...@@ -104,8 +105,8 @@ plot 0 linecolor rgb "black" notitle,\\
''') ''')
if (ref_folder is not None) and (ref_basename is not None): if (ref_folder is not None) and (ref_basename is not None):
plotfile.write('''\ 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 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 errorbars linecolor "black" pointtype 1 linewidth 1 notitle'''+(len(working_basenames)>0)*(''',\\
''')+(len(working_basenames)==0)*(''' ''')+(len(working_basenames)==0)*('''
''')) '''))
...@@ -113,8 +114,8 @@ plot 0 linecolor rgb "black" notitle,\\ ...@@ -113,8 +114,8 @@ plot 0 linecolor rgb "black" notitle,\\
working_basename = working_basenames[k_basename] working_basename = working_basenames[k_basename]
working_scaling = working_scalings[k_basename] working_scaling = working_scalings[k_basename]
plotfile.write('''\ 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 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 errorbars linestyle '''+str(k_basename+1)+''' linewidth 1 notitle'''+(k_basename<len(working_basenames)-1)*(''',\\
''')+(k_basename==len(working_basenames)-1)*(''' ''')+(k_basename==len(working_basenames)-1)*('''
''')) '''))
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment