Mentions légales du service

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

Adding cpp expressions...finally!

parent ef3d0831
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
__init__.py
+ 3
− 1
View file @ 1a3648e0
......@@ -3,8 +3,10 @@
from fedic import *
from plot_strains_vs_radius import *
from plot_strains import *
from image_expressions_cpp import *
from compute_strains import *
from print_tools import *
from plot_displacement_error import *
from image_expressions_py import *
from compute_displacement_error import *
from compute_quadrature_degree import *
from image_expressions import *
compute_strains.py
+ 4
− 3
View file @ 1a3648e0
......@@ -18,6 +18,7 @@ import myVTKPythonLibrary as myVTK
def compute_strains(
sol_folder,
sol_basename,
sol_zfill=6,
sol_ext="vtu",
disp_array_name="displacement",
ref_folder=None,
......@@ -48,12 +49,12 @@ def compute_strains(
ref_mesh = None
n_sector_ids = 0
n_frames = len(glob.glob(sol_folder+"/"+sol_basename+"_*."+sol_ext))
if (verbose): print "n_frames = " + str(n_frames)
n_zfill = len(glob.glob(sol_folder+"/"+sol_basename+"_*."+sol_ext)[0].rsplit("_")[-1].split(".")[0])
if (verbose): print "n_zfill = " + str(n_zfill)
n_frames = len(glob.glob(sol_folder+"/"+sol_basename+"_"+"[0-9]"*sol_zfill+"."+sol_ext))
if (verbose): print "n_frames = " + str(n_frames)
if (write_strains):
strain_file = open(sol_folder+"/"+sol_basename+"-strains.dat", "w")
strain_file.write("#t Err_avg Err_std Ecc_avg Ecc_std Ell_avg Ell_std Erc_avg Erc_std Erl_avg Erl_std Ecl_avg Ecl_std\n")
......
fedic.py
+ 312
− 146
View file @ 1a3648e0
This diff is collapsed.
image_expressions_cpp.py
+ 104
− 392
View file @ 1a3648e0
......@@ -16,412 +16,124 @@ from myVTKPythonLibrary.mat_vec_tools import *
########################################################################
def getScalingFactor(scalar_type_as_string):
if (scalar_type_as_string == 'unsigned char' ): return float(2**8 -1)
elif (scalar_type_as_string == 'unsigned short'): return float(2**16-1)
elif (scalar_type_as_string == 'unsigned int' ): return float(2**32-1)
elif (scalar_type_as_string == 'unsigned long' ): return float(2**64-1)
elif (scalar_type_as_string == 'float' ): return 1.
elif (scalar_type_as_string == 'double' ): return 1.
else: assert (0), "Wrong image scalar type. Aborting."
def get_ExprIm_cpp(
im_dim,
im_type="im",
im_is_def=False,
verbose=False):
class ExprIm2(dolfin.Expression):
def __init__(self, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.X = numpy.array([float()]*2+[Z])
assert (im_dim in (2,3))
assert (im_type in ("im","grad"))
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=0.,
verbose=0)
ExprIm_cpp = '''\
#include <string.h>
def eval(self, Expr, X):
#print "Expr"
self.X[0:2] = X[0:2]
#print " X = " + str(X)
self.interpolator.Interpolate(self.X, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprIm3(dolfin.Expression):
def __init__(self, filename, **kwargs):
if filename is not None:
self.init_image(filename=filename)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=0.,
verbose=0)
def eval(self, Expr, X):
#print "Expr"
#print " X = " + str(X)
self.interpolator.Interpolate(X, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprGradIm2(dolfin.Expression):
def __init__(self, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.X = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=self.s,
verbose=0)
def value_shape(self):
return (2,)
def eval(self, Expr, X):
self.X[0:2] = X[0:2]
self.interpolator.Interpolate(self.X, Expr)
Expr /= self.s
class ExprGradIm3(dolfin.Expression):
def __init__(self, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=self.s,
verbose=0)
def value_shape(self):
return (3,)
def eval(self, Expr, X):
self.interpolator.Interpolate(X, Expr)
Expr /= self.s
class ExprHessIm2(dolfin.Expression):
def __init__(self, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.X = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=self.s,
verbose=0)
def value_shape(self):
return (2,2)
def eval(self, Expr, X):
self.X[0:2] = X[0:2]
self.interpolator.Interpolate(self.X, Expr)
Expr /= self.s
class ExprHessIm3(dolfin.Expression):
def __init__(self, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=self.s,
verbose=0)
def value_shape(self):
return (3,3)
def eval(self, Expr, X):
self.interpolator.Interpolate(X, Expr)
Expr /= self.s
class ExprDefIm2(dolfin.Expression):
def __init__(self, U, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(2)
self.x = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=0.,
verbose=0)
def eval(self, Expr, X):
#print "Expr"
#print " U = " + str(U.vector().array())
#print " X = " + str(X)
#print " UX = " + str(self.UX)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[0:2] = X[0:2] + self.UX[0:2]
#print " x = " + str(self.x)
#print " Expr = " + str(Expr)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprDefIm3(dolfin.Expression):
def __init__(self, U, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(3)
self.x = numpy.empty(3)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=0.,
verbose=0)
def eval(self, Expr, X):
#print "Expr"
#print " U = " + str(U.vector().array())
#print " X = " + str(X)
#print " UX = " + str(self.UX)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[:] = X + self.UX
#print " x = " + str(self.x)
#print " Expr = " + str(Expr)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprGradDefIm2(dolfin.Expression):
def __init__(self, U, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(2)
self.x = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=self.s,
verbose=0)
def value_shape(self):
return (2,)
def eval(self, Expr, X):
#print "Expr"
#print " U = " + str(U.vector().array())
#print " X = " + str(X)
#print " UX = " + str(self.UX)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[0:2] = X[0:2] + self.UX[0:2]
#print " x = " + str(self.x)
#print " Expr = " + str(Expr)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprGradDefIm3(dolfin.Expression):
def __init__(self, U, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(3)
self.x = numpy.empty(3)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=self.s,
verbose=0)
def value_shape(self):
return (3,)
def eval(self, Expr, X):
#print "Expr"
#print " U = " + str(U.vector().array())
#print " X = " + str(X)
#print " UX = " + str(self.UX)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[:] = X + self.UX
#print " x = " + str(self.x)
#print " Expr = " + str(Expr)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprHessDefIm2(dolfin.Expression):
def __init__(self, U, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(2)
self.x = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=self.s,
verbose=0)
def value_shape(self):
return (2,2)
def eval(self, Expr, X):
self.U.eval(self.UX, X)
self.x[0:2] = X[0:2] + self.UX[0:2]
self.interpolator.Interpolate(self.x, Expr)
Expr /= self.s
class ExprHessDefIm3(dolfin.Expression):
def __init__(self, U, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(3)
self.x = numpy.empty(3)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
out_value=self.s,
verbose=0)
def value_shape(self):
return (3,3)
def eval(self, Expr, X):
self.U.eval(self.UX, X)
self.x[:] = X + self.UX
self.interpolator.Interpolate(self.x, Expr)
Expr /= self.s
cppcode_ExprIm='''
#include <vtkSmartPointer.h>
#include <vtkXMLImageDataReader.h>
#include <vtkImageData.h>
#include <vtkImageData.h>'''+('''
#include <vtkImageGradient.h>''')*(im_type=="grad")+'''
#include <vtkImageInterpolator.h>
namespace dolfin
{
class MyFun : public Expression
class MyExpr : public Expression
{
public:
MyFun(): Expression()
vtkSmartPointer<vtkImageInterpolator> interpolator;'''+('''
const Function* U;
mutable Array<double> UX;
mutable Array<double> x;''')*(im_is_def)+('''
mutable Array<double> X;''')*(not im_is_def)*(im_dim==2)+'''
public:
MyExpr():
Expression('''+('''2''')*(im_type=="grad")+''')'''+(''',
UX('''+str(im_dim)+'''),
x(3)''')*(im_is_def)+(''',
X(3)''')*(not im_is_def)*(im_dim==2)+'''
{
}
void init('''+('''
const Function &UU,''')*(im_is_def)+('''
const double &Z=0.''')*(im_dim==2)+''')
{'''+('''
U = &UU;''')*(im_is_def)+('''
x[2] = Z;''')*(im_is_def)*(im_dim==2)+('''
X[2] = Z;''')*(not im_is_def)*(im_dim==2)+'''
}
void init_image(
const char* filename,
const char* interpol_mode="'''+('''linear''')*(im_type=="im")+('''nearest''')*(im_type=="grad")+'''",
const double &interpol_out_value='''+('''0.''')*(im_type=="im")+('''1.''')*(im_type=="grad")+''')
{
vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New();
reader->SetFileName(filename);
reader->Update();'''+('''
vtkSmartPointer<vtkImageGradient> gradient = vtkSmartPointer<vtkImageGradient>::New();
#if VTK_MAJOR_VERSION >= 6
gradient->SetInputData(reader->GetOutput());
#else
gradient->SetInput(reader->GetOutput());
#endif
gradient->SetDimensionality('''+str(im_dim)+''');
gradient->Update();''')*(im_type=="grad")+'''
interpolator = vtkSmartPointer<vtkImageInterpolator>::New();
if (strcmp(interpol_mode, "nearest") == 0)
{
vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New();
reader->SetFileName("<<filename>>");
reader->Update();
interpolator->Initialize(reader->GetOutput());
interpolator->Update();
};
void eval(Array<double>& values, const Array<double>& x) const
interpolator->SetInterpolationModeToNearest();
}
else if (strcmp(interpol_mode, "linear") == 0)
{
interpolator->Interpolate(x.data(), values.data());
interpolator->SetInterpolationModeToLinear();
}
private:
vtkSmartPointer<vtkImageInterpolator> interpolator = vtkSmartPointer<vtkImageInterpolator>::New();
else if (strcmp(interpol_mode, "cubic") == 0)
{
interpolator->SetInterpolationModeToCubic();
}
else
{
std::cout << "Interpolator interpol_mode (" << interpol_mode << ") must be \\"nearest\\", \\"linear\\" or \\"cubic\\". Aborting." << std::endl;
assert(0);
}
interpolator->SetOutValue(interpol_out_value);
interpolator->Initialize('''+('''reader->GetOutput()''')*(im_type=="im")+('''gradient->GetOutput()''')*(im_type=="grad")+''');
interpolator->Update();
}
void eval(Array<double>& values, const Array<double>& position) const
{'''+('''
std::cout << "position = " << position.str(1) << std::endl;''')*(verbose)+('''
U->eval(UX, position);'''+('''
std::cout << "UX = " << UX.str(1) << std::endl;''')*(verbose)+'''
x[0] = position[0] + UX[0];
x[1] = position[1] + UX[1];'''+('''
std::cout << "x = " << x.str(1) << std::endl;''')*(verbose)+'''
interpolator->Interpolate(x.data(), values.data());''')*(im_is_def)*(im_dim==2)+('''
U->eval(UX, position);'''+('''
std::cout << "UX = " << UX.str(1) << std::endl;''')*(verbose)+'''
x = position + UX;'''+('''
std::cout << "x = " << x.str(1) << std::endl;''')*(verbose)+'''
interpolator->Interpolate(x.data(), values.data());''')*(im_is_def)*(im_dim==3)+('''
X[0] = position[0];
X[1] = position[1];'''+('''
std::cout << "X = " << X.str(1) << std::endl;''')*(verbose)+'''
interpolator->Interpolate(X.data(), values.data());''')*(not im_is_def)*(im_dim==2)+('''
interpolator->Interpolate(position.data(), values.data());''')*(not im_is_def)*(im_dim==3)+('''
std::cout << "values = " << values.str(1) << std::endl;''')*(verbose)+'''
}
};
}'''
return ExprIm_cpp
image_expressions_py.py 0 → 100644
+ 412
− 0
View file @ 1a3648e0
#coding=utf8
########################################################################
### ###
### Created by Martin Genet, 2016 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
########################################################################
import dolfin
import numpy
import myVTKPythonLibrary as myVTK
from myVTKPythonLibrary.mat_vec_tools import *
########################################################################
def getScalingFactor(scalar_type_as_string):
if (scalar_type_as_string == 'unsigned char' ): return float(2**8 -1)
elif (scalar_type_as_string == 'unsigned short'): return float(2**16-1)
elif (scalar_type_as_string == 'unsigned int' ): return float(2**32-1)
elif (scalar_type_as_string == 'unsigned long' ): return float(2**64-1)
elif (scalar_type_as_string == 'float' ): return 1.
elif (scalar_type_as_string == 'double' ): return 1.
else: assert (0), "Wrong image scalar type. Aborting."
class ExprIm2(dolfin.Expression):
def __init__(self, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.X = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="linear",
out_value=0.,
verbose=0)
def eval(self, Expr, X):
#print " X = " + str(X)
self.X[0:2] = X[0:2]
#print " X = " + str(self.X)
self.interpolator.Interpolate(self.X, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprIm3(dolfin.Expression):
def __init__(self, filename, **kwargs):
if filename is not None:
self.init_image(filename=filename)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="linear",
out_value=0.,
verbose=0)
def eval(self, Expr, X):
#print " X = " + str(X)
self.interpolator.Interpolate(X, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprGradIm2(dolfin.Expression):
def __init__(self, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.X = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
verbose=0)
def value_shape(self):
return (2,)
def eval(self, Expr, X):
#print " X = " + str(X)
self.X[0:2] = X[0:2]
#print " X = " + str(self.X)
self.interpolator.Interpolate(self.X, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
class ExprGradIm3(dolfin.Expression):
def __init__(self, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
verbose=0)
def value_shape(self):
return (3,)
def eval(self, Expr, X):
#print " X = " + str(X)
self.interpolator.Interpolate(X, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
class ExprHessIm2(dolfin.Expression):
def __init__(self, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.X = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
verbose=0)
def value_shape(self):
return (2,2)
def eval(self, Expr, X):
self.X[0:2] = X[0:2]
self.interpolator.Interpolate(self.X, Expr)
Expr /= self.s
class ExprHessIm3(dolfin.Expression):
def __init__(self, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
verbose=0)
def value_shape(self):
return (3,3)
def eval(self, Expr, X):
self.interpolator.Interpolate(X, Expr)
Expr /= self.s
class ExprDefIm2(dolfin.Expression):
def __init__(self, U, filename=None, scaling=[1.,0.], Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(2)
self.x = numpy.array([float()]*2+[Z])
self.scaling = scaling
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="linear",
out_value=0.,
verbose=0)
def eval(self, Expr, X):
#print " X = " + str(X)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[0:2] = X[0:2] + self.UX[0:2]
#print " x = " + str(self.x)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
Expr *= self.scaling[0]
Expr += self.scaling[1]
#print " Expr = " + str(Expr)
class ExprDefIm3(dolfin.Expression):
def __init__(self, U, filename=None, scaling=[1.,0.], **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(3)
self.x = numpy.empty(3)
self.scaling = scaling
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="linear",
out_value=0.,
verbose=0)
def eval(self, Expr, X):
#print " X = " + str(X)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[:] = X + self.UX
#print " x = " + str(self.x)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
Expr *= self.scaling[0]
Expr += self.scaling[1]
#print " Expr = " + str(Expr)
class ExprGradDefIm2(dolfin.Expression):
def __init__(self, U, filename=None, scaling=[1.,0.], Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(2)
self.x = numpy.array([float()]*2+[Z])
self.scaling = scaling
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
verbose=0)
def value_shape(self):
return (2,)
def eval(self, Expr, X):
#print " X = " + str(X)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[0:2] = X[0:2] + self.UX[0:2]
#print " x = " + str(self.x)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
Expr *= self.scaling[0]
#print " Expr = " + str(Expr)
class ExprGradDefIm3(dolfin.Expression):
def __init__(self, U, filename=None, scaling=[1.,0.], **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(3)
self.x = numpy.empty(3)
self.scaling = scaling
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
verbose=0)
def value_shape(self):
return (3,)
def eval(self, Expr, X):
#print " X = " + str(X)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[:] = X + self.UX
#print " x = " + str(self.x)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
Expr *= self.scaling[0]
#print " Expr = " + str(Expr)
class ExprHessDefIm2(dolfin.Expression):
def __init__(self, U, filename=None, scaling=[1.,0.], Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(2)
self.x = numpy.array([float()]*2+[Z])
self.scaling = scaling
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
verbose=0)
def value_shape(self):
return (2,2)
def eval(self, Expr, X):
self.U.eval(self.UX, X)
self.x[0:2] = X[0:2] + self.UX[0:2]
self.interpolator.Interpolate(self.x, Expr)
Expr /= self.s
Expr *= self.scaling[0]
class ExprHessDefIm3(dolfin.Expression):
def __init__(self, U, filename=None, scaling=[1.,0.], **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(3)
self.x = numpy.empty(3)
self.scaling = scaling
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
verbose=0)
def value_shape(self):
return (3,3)
def eval(self, Expr, X):
self.U.eval(self.UX, X)
self.x[:] = X + self.UX
self.interpolator.Interpolate(self.x, Expr)
Expr /= self.s
Expr *= self.scaling[0]
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