From 09b6808a4b3e41831d11bb65db1cce246b465f8a Mon Sep 17 00:00:00 2001
From: Martin Genet <martin.genet@polytechnique.edu>
Date: Mon, 21 Mar 2016 14:47:09 +0100
Subject: [PATCH] use vector valued quadrature elements for image gradients
---
compute_quadrature_degree.py | 29 ++--
fedic.py | 98 +++++------
image_expressions.py | 324 ++---------------------------------
3 files changed, 78 insertions(+), 373 deletions(-)
diff --git a/compute_quadrature_degree.py b/compute_quadrature_degree.py
index 539098d..7d5c344 100644
--- a/compute_quadrature_degree.py
+++ b/compute_quadrature_degree.py
@@ -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
diff --git a/fedic.py b/fedic.py
index c8612a7..c5ebdee 100644
--- a/fedic.py
+++ b/fedic.py
@@ -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)
diff --git a/image_expressions.py b/image_expressions.py
index 156165a..4f2f94c 100644
--- a/image_expressions.py
+++ b/image_expressions.py
@@ -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>
--
GitLab