Mentions légales du service

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

Adding I1old tangent

parent 3dd63a56
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 115 additions and 100 deletions
fedic.py
+ 115
− 100
View file @ 5612b0a7
......@@ -25,13 +25,14 @@ def fedic(
images_basename,
mesh_folder,
mesh_basename,
images_dimension=3, # 2 or 3
images_zfill=2,
images_k_ref=0,
mesh_degree=1,
images_dimension=3,
k_ref=0,
penalty=0.9,
use_I0_tangent=0,
relax_type="const",
relax_init=0.9,
tangent_type="I1", # I0, I1old or I1
relax_type="const", # const, aitken or manual
relax_init=1.0,
relax_n_iter=10,
tol_res=None,
tol_disp=None,
......@@ -44,140 +45,134 @@ def fedic(
os.mkdir(working_folder)
print "Loading mesh…"
mesh = dolfin.Mesh(mesh_folder+"/"+mesh_basename+".xml")
mesh_filename = mesh_folder+"/"+mesh_basename+".xml"
assert os.path.exists(mesh_filename)
mesh = dolfin.Mesh(mesh_filename)
dX = dolfin.dx(mesh)
mesh_V = dolfin.assemble(dolfin.Constant(1)*dX)
fs = dolfin.VectorFunctionSpace(
vfs = dolfin.VectorFunctionSpace(
mesh=mesh,
family="Lagrange",
degree=mesh_degree)
dX = dolfin.dx(mesh)
print "Computing quadrature degree…"
degree = myFEniCS.compute_quadrature_degree(
image_filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
print "Computing quadrature degree for images…"
ref_image_filename = images_folder+"/"+images_basename+"_"+str(images_k_ref).zfill(images_zfill)+".vti"
images_quadrature = myFEniCS.compute_quadrature_degree(
image_filename=ref_image_filename,
mesh=mesh,
image_dimension=images_dimension)
print "degree = " + str(degree)
#images_quadrature = 1
print "images_quadrature = " + str(images_quadrature)
fe = dolfin.FiniteElement(
family="Quadrature",
cell=mesh.ufl_cell(),
degree=degree,
degree=images_quadrature,
quad_scheme="default")
ve = dolfin.VectorElement(
family="Quadrature",
cell=mesh.ufl_cell(),
degree=degree,
degree=images_quadrature,
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",
filename=ref_image_filename,
element=fe)
DI0 = myFEniCS.ExprGradIm2(
filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
filename=ref_image_filename,
element=ve)
elif (images_dimension == 3):
I0 = myFEniCS.ExprIm3(
filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
filename=ref_image_filename,
element=fe)
DI0 = myFEniCS.ExprGradIm3(
filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
filename=ref_image_filename,
element=ve)
else:
assert (0), "images_dimension must be 2 or 3. Aborting."
I0_norm = dolfin.assemble(I0**2 * dX, form_compiler_parameters={'quadrature_degree':degree})**(1./2)
I0_norm = (dolfin.assemble(I0**2 * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})/mesh_V)**(1./2)
assert (I0_norm > 0.), "I0_norm = " + str(I0_norm)
print "I0_norm = " + str(I0_norm)
print "Defining functions…"
penalty = dolfin.Constant(penalty)
U = dolfin.Function(
fs,
vfs,
name="displacement")
U_old = dolfin.Function(
fs,
name="old displacement")
vfs,
name="previous displacement")
dU = dolfin.Function(
fs,
name="ddisplacement")
ddU = dolfin.TrialFunction(fs)
V = dolfin.TestFunction(fs)
vfs,
name="displacement correction")
ddU = dolfin.TrialFunction(vfs)
V = dolfin.TestFunction(vfs)
print "Printing initial solution…"
file_pvd = dolfin.File(working_folder+"/"+working_basename+"_.pvd")
for vtu in glob.glob(working_folder+"/"+working_basename+"_*.vtu"):
os.remove(vtu)
file_pvd << (U, float(k_ref))
pvd_basename = working_folder+"/"+working_basename+"_"
file_pvd = dolfin.File(pvd_basename+".pvd")
for vtu_filename in glob.glob(pvd_basename+"*.vtu"):
os.remove(vtu_filename)
file_pvd << (U, float(images_k_ref))
print "Defining variational forms…"
if (images_dimension == 2):
I1 = myFEniCS.ExprDefIm2(
U=U,
filename=ref_image_filename,
element=fe)
DI1 = myFEniCS.ExprGradDefIm2(
U=U,
filename=ref_image_filename,
element=ve)
if (use_I0_tangent):
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(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(DI1, V))*dX\
- (1.-penalty) * dolfin.inner(dolfin.grad(U),
dolfin.grad(V))*dX
b0 = penalty * dolfin.inner(I0,
dolfin.dot(DI0, V))*dX
elif (images_dimension == 3):
I1 = myFEniCS.ExprDefIm3(
U=U,
filename=ref_image_filename,
element=fe)
DI1 = myFEniCS.ExprGradDefIm3(
U=U,
filename=ref_image_filename,
element=ve)
if (use_I0_tangent):
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(DI1, ddU),
dolfin.dot(DI1, V))*dX\
+ (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
dolfin.grad( V))*dX
b = penalty * dolfin.inner(I0-I1,
else:
assert (0), "images_dimension must be 2 or 3. Aborting."
if (tangent_type == "I0"):
a = penalty * dolfin.inner(dolfin.dot(DI0, ddU),
dolfin.dot(DI0, V))*dX
elif (tangent_type == "I1") or (tangent_type == "I1old"):
a = penalty * dolfin.inner(dolfin.dot(DI1, ddU),
dolfin.dot(DI1, V))*dX
a += (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
dolfin.grad( V))*dX
b = penalty * dolfin.inner(I0-I1,
dolfin.dot(DI1, V))*dX\
- (1.-penalty) * dolfin.inner(dolfin.grad(U),
dolfin.grad(V))*dX
b0 = dolfin.inner(I0,
b0 = dolfin.inner(I0,
dolfin.dot(DI0, V))*dX
else:
assert (0), "images_dimension must be 2 or 3. Aborting."
B0 = dolfin.assemble(b0, form_compiler_parameters={'quadrature_degree':degree})
B0 = dolfin.assemble(b0, form_compiler_parameters={'quadrature_degree':images_quadrature})
B0_norm = numpy.linalg.norm(B0)
assert (B0_norm > 0.), "B0_norm = " + str(B0_norm)
print "B0_norm = " + str(B0_norm)
if (use_I0_tangent):
A = dolfin.assemble(a, form_compiler_parameters={'quadrature_degree':degree})
else:
A = None
# linear system
A = None
if (tangent_type == "I0"):
A = dolfin.assemble(a, tensor=A, form_compiler_parameters={'quadrature_degree':images_quadrature})
B = None
print "Checking number of frames…"
n_frames = len(glob.glob(images_folder+"/"+images_basename+"_*.vti"))
assert (n_frames > 1), "n_frames = " + str(n_frames)
n_frames = 2
assert (n_frames > 1), "n_frames = " + str(n_frames) + " <= 1. Aborting."
print "n_frames = " + str(n_frames)
assert (abs(k_ref) < n_frames), "k_ref = " + str(k_ref)
k_ref = k_ref%n_frames
assert (abs(images_k_ref) < n_frames), "abs(images_k_ref) = " + str(images_k_ref) + " >= n_frames. Aborting."
images_k_ref = images_k_ref%n_frames
print "images_k_ref = " + str(images_k_ref)
n_iter_tot = 0
global_success = True
......@@ -185,36 +180,53 @@ def fedic(
print "forward_or_backward = " + forward_or_backward
if (forward_or_backward == "forward"):
k_frames = range(k_ref+1, n_frames, +1)
k_frames = range(images_k_ref+1, n_frames, +1)
elif (forward_or_backward == "backward"):
k_frames = range(k_ref-1, -1, -1)
k_frames = range(images_k_ref-1, -1, -1)
print "k_frames = " + str(k_frames)
if (forward_or_backward == "backward"):
U.vector()[:] = 0.
I1.init_image( filename=ref_image_filename)
DI1.init_image(filename=ref_image_filename)
for k_frame in k_frames:
print "k_frame = " + str(k_frame)
if (print_iterations):
if (os.path.exists(working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(2)+".pdf")):
os.remove(working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(2)+".pdf")
file_dat_iter = open(working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(2)+".dat", "w")
file_pvd_iter = dolfin.File(working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(2)+"_.pvd")
for vtu in glob.glob(working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(2)+"_*.vtu"):
os.remove(vtu)
file_pvd_iter << (U, 0.)
frame_basename = working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(images_zfill)
if (os.path.exists(frame_basename+".pdf")):
os.remove(frame_basename+".pdf")
file_dat_frame = open(frame_basename+".dat", "w")
file_dat_frame.write("#k_iter B_norm err_res relax dU_norm U_norm err_disp I1I0_norm err_im err_im_rel\n")
file_pvd_frame = dolfin.File(frame_basename+"_.pvd")
for vtu_filename in glob.glob(frame_basename+"_*.vtu"):
os.remove(vtu_filename)
file_pvd_frame << (U, 0.)
# linear system: matrix
if (tangent_type == "I1old"):
A = dolfin.assemble(a, tensor=A, form_compiler_parameters={'quadrature_degree':images_quadrature})
#print "A = " + str(A.array())
#A_norm = numpy.linalg.norm(A.array())
#print "A_norm = " + str(A_norm)
print "Loading image and image gradient…"
if (images_dimension == 2):
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")
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."
image_filename = images_folder+"/"+images_basename+"_"+str(k_frame).zfill(images_zfill)+".vti"
I1.init_image( filename=image_filename)
DI1.init_image(filename=image_filename)
if (print_iterations):
U_old.vector()[:] = U.vector()[:]
U.vector()[:] = 0.
I1I0_norm = (dolfin.assemble((I1-I0)**2 * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})/mesh_V)**(1./2)
err_im = I1I0_norm/I0_norm
file_dat_frame.write(" ".join([str(val) for val in [-2, None, None, None, None, None, None, I1I0_norm, err_im, None]])+"\n")
U.vector()[:] = U_old.vector()[:]
I1I0_norm = (dolfin.assemble((I1-I0)**2 * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})/mesh_V)**(1./2)
err_im = I1I0_norm/I0_norm
file_dat_frame.write(" ".join([str(val) for val in [-1, None, None, None, None, None, None, I1I0_norm, err_im, None]])+"\n")
print "Running registration…"
k_iter = 0
......@@ -224,8 +236,8 @@ def fedic(
n_iter_tot += 1
# linear system: matrix
if not (use_I0_tangent):
A = dolfin.assemble(a, tensor=A, form_compiler_parameters={'quadrature_degree':degree})
if (tangent_type == "I1"):
A = dolfin.assemble(a, tensor=A, form_compiler_parameters={'quadrature_degree':images_quadrature})
#print "A = " + str(A.array())
#A_norm = numpy.linalg.norm(A.array())
#print "A_norm = " + str(A_norm)
......@@ -237,7 +249,7 @@ def fedic(
elif (k_iter > 1):
B_old[:] = B[:]
#print "B_old = " + str(B_old[0])
B = dolfin.assemble(b, tensor=B, form_compiler_parameters={'quadrature_degree':degree})
B = dolfin.assemble(b, tensor=B, form_compiler_parameters={'quadrature_degree':images_quadrature})
#print "B = " + str(B.array())
# residual error
......@@ -265,29 +277,32 @@ def fedic(
dB[:] = B[:] - B_old[:]
relax *= (-1.) * B_old.inner(dB) / dB.inner(dB)
print "relax = " + str(relax)
elif (relax_type == "automatic"):
elif (relax_type == "manual"):
B_relax = numpy.empty(relax_n_iter)
for relax_k in xrange(relax_n_iter):
#print "relax_k = " + str(relax_k)
relax = float(relax_k+1)/relax_n_iter
U.vector()[:] = U_old.vector()[:] + relax * dU.vector()[:]
B = dolfin.assemble(b, tensor=B, form_compiler_parameters={'quadrature_degree':degree})
B = dolfin.assemble(b, tensor=B, form_compiler_parameters={'quadrature_degree':images_quadrature})
B_relax[relax_k] = numpy.linalg.norm(B)
#print "B_relax = " + str(B_relax[relax_k])
#print "B_relax = " + str(B_relax)
if (print_iterations):
open(working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(2)+"-iter="+str(k_iter).zfill(3)+".dat", "w").write("\n".join([" ".join([str(val) for val in [float(relax_k+1)/relax_n_iter, B_relax[relax_k]]]) for relax_k in xrange(relax_n_iter)]))
os.system("gnuplot -e \"set terminal pdf; set output '"+working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(2)+"-iter="+str(k_iter).zfill(3)+".pdf'; plot '"+working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(2)+"-iter="+str(k_iter).zfill(3)+".dat' u 1:2 w l notitle\"")
iter_basename = frame_basename+"-iter="+str(k_iter).zfill(3)
open(iter_basename+".dat", "w").write("\n".join([" ".join([str(val) for val in [float(relax_k+1)/relax_n_iter, B_relax[relax_k]]]) for relax_k in xrange(relax_n_iter)]))
os.system("gnuplot -e \"set terminal pdf; set output '"+iter_basename+".pdf'; plot '"+iter_basename+".dat' u 1:2 w l notitle\"")
relax_k = numpy.argmin(B_relax)
relax = float(relax_k+1)/relax_n_iter
print "relax = " + str(relax)
else:
assert (0), "relax_type must be \"constant\", \"aitken\" or \"manual\". Aborting."
# solution update
U.vector()[:] = U_old.vector()[:] + relax * dU.vector()[:]
if (print_iterations):
#print "U = " + str(U.vector().array())
file_pvd_iter << (U, float(k_iter+1))
file_pvd_frame << (U, float(k_iter+1))
# displacement error
dU_norm = numpy.linalg.norm(dU.vector().array())
......@@ -301,7 +316,7 @@ def fedic(
# image error
if (k_iter > 0):
I1I0_norm_old = I1I0_norm
I1I0_norm = dolfin.assemble((I1-I0)**2 * dX, form_compiler_parameters={'quadrature_degree':degree})**(1./2)
I1I0_norm = (dolfin.assemble((I1-I0)**2 * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})/mesh_V)**(1./2)
#print "I1I0_norm = " + str(I1I0_norm)
err_im = I1I0_norm/I0_norm
print "err_im = " + str(err_im)
......@@ -312,7 +327,7 @@ def fedic(
print "err_im_rel = " + str(err_im_rel)
if (print_iterations):
file_dat_iter.write(" ".join([str(val) for val in [k_iter, B_norm, err_res, relax, dU_norm, U_norm, err_disp, I1I0_norm, err_im, err_im_rel]])+"\n")
file_dat_frame.write(" ".join([str(val) for val in [k_iter, B_norm, err_res, relax, dU_norm, U_norm, err_disp, I1I0_norm, err_im, err_im_rel]])+"\n")
# exit test
if (tol_disp is not None) and (tol_res is not None):
......@@ -337,9 +352,9 @@ def fedic(
k_iter += 1
if (print_iterations):
os.remove(working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(2)+"_.pvd")
file_dat_iter.close()
os.system("gnuplot -e \"set terminal pdf; set output '"+working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(2)+".pdf'; set logscale y; set yrange [1e-4:1e1]; plot '"+working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(2)+".dat' u 1:3 pt 1 lw 3 title 'err_res', '' u 1:4 pt 1 lw 3 title 'relax', '' u 1:6 pt 1 lw 3 title 'err_disp', '' using 1:8 pt 1 lw 3 title 'err_im', '' using 1:9 pt 1 lw 3 title 'err_im_rel', 1e-3 lt 0 notitle\"")
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-4:1e1]; plot '"+frame_basename+".dat' u 1:3 pt 1 lw 3 title 'err_res', '' u 1:4 pt 1 lw 3 title 'relax', '' u 1:7 pt 1 lw 3 title 'err_disp', '' using 1:9 pt 1 lw 3 title 'err_im', '' using 1:10 pt 1 lw 3 title 'err_im_rel', "+str(tol_res)+" lt -1 notitle\"")
if not (success) and not (continue_after_fail):
break
......@@ -352,6 +367,6 @@ def fedic(
print "n_iter_tot = " + str(n_iter_tot)
os.remove(working_folder+"/"+working_basename+"_.pvd")
os.remove(pvd_basename+".pvd")
return global_success
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