Mentions légales du service

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

use vector valued quadrature elements for image gradients

parent c6c13f1b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
compute_quadrature_degree.py
+ 17
− 12
View file @ 09b6808a
......@@ -16,7 +16,7 @@ import myFEniCSPythonLibrary as myFEniCS
def compute_quadrature_degree(
image_filename,
dX,
mesh,
image_dimension=3,
deg_min=1,
deg_max=10,
......@@ -24,13 +24,17 @@ def compute_quadrature_degree(
n_under_tol=1,
verbose=1):
dX = dolfin.dx(mesh)
first_time = True
k_under_tol = 0
for degree in xrange(deg_min,deg_max+1):
if (verbose): print "degree = " + str(degree)
fe = dolfin.FiniteElement(
family="Quadrature",
degree=degree)
cell=mesh.ufl_cell(),
degree=degree,
quad_scheme="default")
if (image_dimension == 2):
I0 = myFEniCS.ExprIm2(
filename=image_filename,
......@@ -41,18 +45,19 @@ def compute_quadrature_degree(
element=fe)
else:
assert (0), "image_dimension must be 2 or 3. Aborting."
I0_norm = dolfin.assemble(I0**2 * dX)**(1./2)
if not (first_time):
I0_norm_old = I0_norm
I0_norm = dolfin.assemble(I0**2 * dX, form_compiler_parameters={'quadrature_degree':degree})**(1./2)
if (verbose): print "I0_norm = " + str(I0_norm)
if (first_time):
first_time = False
continue
I0_norm_err = abs(I0_norm-I0_norm_old)/I0_norm_old
if (verbose): print "I0_norm_err = " + str(I0_norm_err)
if (I0_norm_err < tol):
k_under_tol += 1
else:
I0_norm_err = abs(I0_norm-I0_norm_old)/I0_norm_old
if (verbose): print "I0_norm_err = " + str(I0_norm_err)
if (I0_norm_err < tol):
k_under_tol += 1
else:
k_under_tol = 0
if (k_under_tol >= n_under_tol):
break
I0_norm_old = I0_norm
k_under_tol = 0
if (k_under_tol >= n_under_tol):
break
return degree
fedic.py
+ 43
− 55
View file @ 09b6808a
......@@ -55,41 +55,39 @@ def fedic(
print "Computing quadrature degree…"
degree = myFEniCS.compute_quadrature_degree(
image_filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
dX=dX,
mesh=mesh,
image_dimension=images_dimension)
print "degree = " + str(degree)
fe = dolfin.FiniteElement(
family="Quadrature",
degree=degree)
cell=mesh.ufl_cell(),
degree=degree,
quad_scheme="default")
ve = dolfin.VectorElement(
family="Quadrature",
cell=mesh.ufl_cell(),
degree=degree,
quad_scheme="default")
print "Loading reference image…"
if (images_dimension == 2):
I0 = myFEniCS.ExprIm2(
filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
element=fe)
DXI0 = myFEniCS.ExprGradXIm2(
DI0 = myFEniCS.ExprGradIm2(
filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
element=fe)
DYI0 = myFEniCS.ExprGradYIm2(
filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
element=fe)
element=ve)
elif (images_dimension == 3):
I0 = myFEniCS.ExprIm3(
filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
element=fe)
DXI0 = myFEniCS.ExprGradXIm3(
DI0 = myFEniCS.ExprGradIm3(
filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
element=fe)
DYI0 = myFEniCS.ExprGradYIm3(
filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
element=fe)
DZI0 = myFEniCS.ExprGradZIm3(
filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
element=fe)
element=ve)
else:
assert (0), "images_dimension must be 2 or 3. Aborting."
I0_norm = dolfin.assemble(I0**2 * dX)**(1./2)
I0_norm = dolfin.assemble(I0**2 * dX, form_compiler_parameters={'quadrature_degree':degree})**(1./2)
assert (I0_norm > 0.), "I0_norm = " + str(I0_norm)
print "I0_norm = " + str(I0_norm)
......@@ -98,7 +96,9 @@ def fedic(
U = dolfin.Function(
fs,
name="displacement")
U_old = dolfin.Function(U)
U_old = dolfin.Function(
fs,
name="old displacement")
dU = dolfin.Function(
fs,
name="ddisplacement")
......@@ -116,67 +116,58 @@ def fedic(
I1 = myFEniCS.ExprDefIm2(
U=U,
element=fe)
DXI1 = myFEniCS.ExprGradXDefIm2(
U=U,
element=fe)
DYI1 = myFEniCS.ExprGradYDefIm2(
DI1 = myFEniCS.ExprGradDefIm2(
U=U,
element=fe)
element=ve)
if (use_I0_tangent):
a = penalty * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI0, DYI0)), ddU),
dolfin.dot(dolfin.as_vector((DXI0, DYI0)), V))*dX\
a = penalty * dolfin.inner(dolfin.dot(DI0, ddU),
dolfin.dot(DI0, V))*dX\
+ (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
dolfin.grad( V))*dX
else:
a = penalty * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI1, DYI1)), ddU),
dolfin.dot(dolfin.as_vector((DXI1, DYI1)), V))*dX\
a = penalty * dolfin.inner(dolfin.dot(DI1, ddU),
dolfin.dot(DI1, V))*dX\
+ (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
dolfin.grad( V))*dX
b = penalty * dolfin.inner(I0-I1,
dolfin.dot(dolfin.as_vector((DXI1, DYI1)), V))*dX\
dolfin.dot(DI1, V))*dX\
- (1.-penalty) * dolfin.inner(dolfin.grad(U),
dolfin.grad(V))*dX
b0 = penalty * dolfin.inner(I0,
dolfin.dot(dolfin.as_vector((DXI0, DYI0)), V))*dX
dolfin.dot(DI0, V))*dX
elif (images_dimension == 3):
I1 = myFEniCS.ExprDefIm3(
U=U,
element=fe)
DXI1 = myFEniCS.ExprGradXDefIm3(
DI1 = myFEniCS.ExprGradDefIm3(
U=U,
element=fe)
DYI1 = myFEniCS.ExprGradYDefIm3(
U=U,
element=fe)
DZI1 = myFEniCS.ExprGradZDefIm3(
U=U,
element=fe)
element=ve)
if (use_I0_tangent):
a = penalty * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI0, DYI0, DZI0)), ddU),
dolfin.dot(dolfin.as_vector((DXI0, DYI0, DZI0)), V))*dX\
a = penalty * dolfin.inner(dolfin.dot(DI0, ddU),
dolfin.dot(DI0, V))*dX\
+ (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
dolfin.grad( V))*dX
else:
a = penalty * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI1, DYI1, DZI1)), ddU),
dolfin.dot(dolfin.as_vector((DXI1, DYI1, DZI1)), V))*dX\
a = penalty * dolfin.inner(dolfin.dot(DI1, ddU),
dolfin.dot(DI1, V))*dX\
+ (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
dolfin.grad( V))*dX
b = penalty * dolfin.inner(I0-I1,
dolfin.dot(dolfin.as_vector((DXI1, DYI1, DZI1)), V))*dX\
dolfin.dot(DI1, V))*dX\
- (1.-penalty) * dolfin.inner(dolfin.grad(U),
dolfin.grad(V))*dX
b0 = penalty * dolfin.inner(I0,
dolfin.dot(dolfin.as_vector((DXI0, DYI0, DZI0)), V))*dX
dolfin.dot(DI0, V))*dX
else:
assert (0), "images_dimension must be 2 or 3. Aborting."
B0 = dolfin.assemble(b0)
B0 = dolfin.assemble(b0, form_compiler_parameters={'quadrature_degree':degree})
B0_norm = numpy.linalg.norm(B0)
assert (B0_norm > 0.), "B0_norm = " + str(B0_norm)
print "B0_norm = " + str(B0_norm)
# linear system
if (use_I0_tangent):
A = dolfin.assemble(a)
A = dolfin.assemble(a, form_compiler_parameters={'quadrature_degree':degree})
else:
A = None
B = None
......@@ -218,14 +209,11 @@ def fedic(
print "Loading image and image gradient…"
if (images_dimension == 2):
I1.init_image( filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(2)+".vti")
DXI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(2)+".vti")
DYI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(2)+".vti")
I1.init_image( filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(2)+".vti")
DI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(2)+".vti")
elif (images_dimension == 3):
I1.init_image( filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(2)+".vti")
DXI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(2)+".vti")
DYI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(2)+".vti")
DZI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(2)+".vti")
I1.init_image( filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(2)+".vti")
DI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(2)+".vti")
else:
assert (0), "images_dimension must be 2 or 3. Aborting."
......@@ -238,7 +226,7 @@ def fedic(
# linear system: matrix
if not (use_I0_tangent):
A = dolfin.assemble(a, tensor=A)
A = dolfin.assemble(a, tensor=A, form_compiler_parameters={'quadrature_degree':degree})
#print "A = " + str(A.array())
#A_norm = numpy.linalg.norm(A.array())
#print "A_norm = " + str(A_norm)
......@@ -250,7 +238,7 @@ def fedic(
elif (k_iter > 1):
B_old[:] = B[:]
#print "B_old = " + str(B_old[0])
B = dolfin.assemble(b, tensor=B)
B = dolfin.assemble(b, tensor=B, form_compiler_parameters={'quadrature_degree':degree})
#print "B = " + str(B.array())
B_norm = numpy.linalg.norm(B.array())
#print "B_norm = " + str(B_norm)
......@@ -282,7 +270,7 @@ def fedic(
#print "k_relax = " + str(k_relax)
relax = float(k_relax+1)/relax_n_iter
U.vector()[:] = U_old.vector()[:] + relax * dU.vector()[:]
B = dolfin.assemble(b, tensor=B)
B = dolfin.assemble(b, tensor=B, form_compiler_parameters={'quadrature_degree':degree})
B_relax[k_relax] = numpy.linalg.norm(B)
#print "B_relax = " + str(B_relax[k_relax])
#print "B_relax = " + str(B_relax)
......@@ -309,7 +297,7 @@ def fedic(
# image error
if (k_iter > 0):
I1I0_norm_old = I1I0_norm
I1I0_norm = dolfin.assemble((I1-I0)**2 * dX)**(1./2)
I1I0_norm = dolfin.assemble((I1-I0)**2 * dX, form_compiler_parameters={'quadrature_degree':degree})**(1./2)
#print "I1I0_norm = " + str(I1I0_norm)
err_im = I1I0_norm/I0_norm
print "err_im = " + str(err_im)
......
image_expressions.py
+ 18
− 306
View file @ 09b6808a
......@@ -36,14 +36,14 @@ class ExprIm2(dolfin.Expression):
filename=filename,
verbose=0)
self.s = getScalingFactor(self.image.GetScalarTypeAsString())
#print self.s
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def eval(self, Im, X):
#print "Im"
self.X[0] = X[0]
self.X[1] = X[1]
self.X[0:2] = X[0:2]
#print " X = " + str(X)
self.interpolator.Interpolate(self.X, Im)
#print " Im = " + str(Im)
......@@ -62,6 +62,7 @@ class ExprIm3(dolfin.Expression):
filename=filename,
verbose=0)
self.s = getScalingFactor(self.image.GetScalarTypeAsString())
#print self.s
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
......@@ -75,7 +76,7 @@ class ExprIm3(dolfin.Expression):
#print " Im = " + str(Im)
#self.evalcount += 1
class ExprGradIm2(dolfin.Expression): # for some reason this does not work with quadrature finite elements, hence the following scalar gradients definition
class ExprGradIm2(dolfin.Expression):
def __init__(self, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
......@@ -97,60 +98,11 @@ class ExprGradIm2(dolfin.Expression): # for some reason this does not work with
return (2,)
def eval(self, GradIm, X):
self.X[0] = X[0]
self.X[1] = X[1]
self.X[0:2] = X[0:2]
self.interpolator.Interpolate(self.X, GradIm)
GradIm /= self.s
class ExprGradXIm2(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(self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def eval(self, GradXIm, X):
self.X[0] = X[0]
self.X[1] = X[1]
GradXIm[0] = self.interpolator.Interpolate(self.X[0], self.X[1], self.X[2], 0)
GradXIm /= self.s
class ExprGradYIm2(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(self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def eval(self, GradYIm, X):
self.X[0] = X[0]
self.X[1] = X[1]
GradYIm[0] = self.interpolator.Interpolate(self.X[0], self.X[1], self.X[2], 1)
GradYIm /= self.s
class ExprGradIm3(dolfin.Expression): # for some reason this does not work with quadrature finite elements, hence the following scalar gradients definition
class ExprGradIm3(dolfin.Expression):
def __init__(self, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
......@@ -174,75 +126,12 @@ class ExprGradIm3(dolfin.Expression): # for some reason this does not work with
self.interpolator.Interpolate(X, GradIm)
GradIm /= self.s
class ExprGradXIm3(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(self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def eval(self, GradXIm, X):
GradXIm[0] = self.interpolator.Interpolate(X[0], X[1], X[2], 0)
GradXIm /= self.s
class ExprGradYIm3(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(self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def eval(self, GradYIm, X):
GradYIm[0] = self.interpolator.Interpolate(X[0], X[1], X[2], 1)
GradYIm /= self.s
class ExprGradZIm3(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(self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def eval(self, GradZIm, X):
GradZIm[0] = self.interpolator.Interpolate(X[0], X[1], X[2], 2)
GradZIm /= 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.array([float()]*2)
self.UX = numpy.empty(2)
self.x = numpy.array([float()]*2+[Z])
def init_image(self, filename):
......@@ -261,8 +150,7 @@ class ExprDefIm2(dolfin.Expression):
#print " UX = " + str(self.UX)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[0] = X[0] + self.UX[0]
self.x[1] = X[1] + self.UX[1]
self.x[0:2] = X[0:2] + self.UX[0:2]
#print " x = " + str(self.x)
#print " DefIm = " + str(DefIm)
self.interpolator.Interpolate(self.x, DefIm)
......@@ -275,14 +163,15 @@ class ExprDefIm3(dolfin.Expression):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.array([float()]*3)
self.x = numpy.array([float()]*3)
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(self.image.GetScalarTypeAsString())
#print self.s
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
......@@ -302,12 +191,12 @@ class ExprDefIm3(dolfin.Expression):
DefIm /= self.s
#print " DefIm = " + str(DefIm)
class ExprGradDefIm2(dolfin.Expression): # for some reason this does not work with quadrature finite elements, hence the following scalar gradients definition
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.array([float()]*2)
self.UX = numpy.empty(2)
self.x = numpy.array([float()]*2+[Z])
def init_image(self, filename):
......@@ -332,8 +221,7 @@ class ExprGradDefIm2(dolfin.Expression): # for some reason this does not work wi
#print " UX = " + str(self.UX)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[0] = X[0] + self.UX[0]
self.x[1] = X[1] + self.UX[1]
self.x[0:2] = X[0:2] + self.UX[0:2]
#print " x = " + str(self.x)
#print " GradDefIm = " + str(GradDefIm)
self.interpolator.Interpolate(self.x, GradDefIm)
......@@ -341,91 +229,20 @@ class ExprGradDefIm2(dolfin.Expression): # for some reason this does not work wi
GradDefIm /= self.s
#print " GradDefIm = " + str(GradDefIm)
class ExprGradXDefIm2(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.array([float()]*2)
self.x = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def eval(self, GradXDefIm, X):
#print "GradXDefIm"
#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] = X[0] + self.UX[0]
self.x[1] = X[1] + self.UX[1]
#print " x = " + str(self.x)
#print " GradXDefIm = " + str(GradXDefIm)
GradXDefIm[0] = self.interpolator.Interpolate(self.x[0], self.x[1], self.x[2], 0)
#print " GradXDefIm = " + str(GradXDefIm)
GradXDefIm /= self.s
#print " GradXDefIm = " + str(GradXDefIm)
class ExprGradYDefIm2(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.array([float()]*2)
self.x = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def eval(self, GradYDefIm, X):
#print "GradYDefIm"
#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] = X[0] + self.UX[0]
self.x[1] = X[1] + self.UX[1]
#print " x = " + str(self.x)
#print " GradYDefIm = " + str(GradYDefIm)
GradYDefIm[0] = self.interpolator.Interpolate(self.x[0], self.x[1], self.x[2], 1)
#print " GradYDefIm = " + str(GradYDefIm)
GradYDefIm /= self.s
#print " GradYDefIm = " + str(GradYDefIm)
class ExprGradDefIm3(dolfin.Expression): # for some reason this does not work with quadrature finite elements, hence the following scalar gradients definition
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.array([float()]*3)
self.x = numpy.array([float()]*3)
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(self.image.GetScalarTypeAsString())
#print self.s
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
......@@ -451,111 +268,6 @@ class ExprGradDefIm3(dolfin.Expression): # for some reason this does not work wi
GradDefIm /= self.s
#print " GradDefIm = " + str(GradDefIm)
class ExprGradXDefIm3(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.array([float()]*3)
self.x = numpy.array([float()]*3)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def eval(self, GradXDefIm, X):
#print "GradXDefIm"
#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 " GradXDefIm = " + str(GradXDefIm)
GradXDefIm[0] = self.interpolator.Interpolate(self.x[0], self.x[1], self.x[2], 0)
#print " GradXDefIm = " + str(GradXDefIm)
GradXDefIm /= self.s
#print " GradXDefIm = " + str(GradXDefIm)
class ExprGradYDefIm3(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.array([float()]*3)
self.x = numpy.array([float()]*3)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def eval(self, GradYDefIm, X):
#print "GradYDefIm"
#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 " GradYDefIm = " + str(GradYDefIm)
GradYDefIm[0] = self.interpolator.Interpolate(self.x[0], self.x[1], self.x[2], 1)
#print " GradYDefIm = " + str(GradYDefIm)
GradYDefIm /= self.s
#print " GradYDefIm = " + str(GradYDefIm)
class ExprGradZDefIm3(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.array([float()]*3)
self.x = numpy.array([float()]*3)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def eval(self, GradZDefIm, X):
#print "GradZDefIm"
#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 " GradZDefIm = " + str(GradZDefIm)
GradZDefIm[0] = self.interpolator.Interpolate(self.x[0], self.x[1], self.x[2], 2)
#print " GradZDefIm = " + str(GradZDefIm)
GradZDefIm /= self.s
#print " GradZDefIm = " + str(GradZDefIm)
cppcode_ExprIm='''
#include <vtkSmartPointer.h>
#include <vtkXMLImageDataReader.h>
......
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