Mentions légales du service

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

New print commands

parent ab2383d6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 69 additions and 48 deletions
fedic.py
+ 69
− 48
View file @ c870881c
......@@ -18,6 +18,15 @@ import myVTKPythonLibrary as myVTK
########################################################################
def print_str(tab, string):
print " | "*tab + string
def print_var(tab, name, val):
print " | "*tab + name + " = " + str(val)
def print_sci(tab, name, val):
print " | "*tab + name.ljust(12) + " = " + format(val,".4e")
def fedic(
working_folder,
working_basename,
......@@ -40,22 +49,24 @@ def fedic(
print_iterations=0,
continue_after_fail=0):
print "Loading mesh…"
tab = 0
print_str(tab,"Loading mesh…")
mesh_filename = mesh_folder+"/"+mesh_basename+".xml"
assert os.path.exists(mesh_filename), "No mesh in "+mesh_filename+". Aborting."
mesh = dolfin.Mesh(mesh_filename)
dX = dolfin.dx(mesh)
mesh_volume = dolfin.assemble(dolfin.Constant(1)*dX)
print "Computing quadrature degree for images…"
print_str(tab,"Computing quadrature degree for images…")
ref_image_filename = images_folder+"/"+images_basename+"_"+str(images_k_ref).zfill(images_zfill)+".vti"
if (images_quadrature is None):
images_quadrature = myFEniCS.compute_quadrature_degree(
image_filename=ref_image_filename,
mesh=mesh)
print "images_quadrature = " + str(images_quadrature)
mesh=mesh,
verbose=0)
print_var(tab+1,"images_quadrature",images_quadrature)
print "Loading reference image…"
print_str(tab,"Loading reference image…")
images_dimension = myVTK.computeImageDimensionality(
image_filename=ref_image_filename,
verbose=0)
......@@ -75,6 +86,9 @@ def fedic(
cell=mesh.ufl_cell(),
degree=images_quadrature,
quad_scheme="default")
te._quad_scheme = "default"
for k in xrange(images_dimension**2):
te.sub_elements()[k]._quad_scheme = "default"
if (images_dimension == 2):
Iref = myFEniCS.ExprIm2(
filename=ref_image_filename,
......@@ -91,9 +105,9 @@ def fedic(
element=ve)
Iref_norm = (dolfin.assemble(Iref**2 * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})/mesh_volume)**(1./2)
assert (Iref_norm > 0.), "Iref_norm = "+str(Iref_norm)+" <= 0. Aborting."
print "Iref_norm = " + str(Iref_norm)
print_var(tab+1,"Iref_norm",Iref_norm)
print "Defining functions…"
print_str(tab,"Defining functions…")
vfs = dolfin.VectorFunctionSpace(
mesh=mesh,
family="Lagrange",
......@@ -113,7 +127,7 @@ def fedic(
dU_ = dolfin.TrialFunction(vfs)
dV_ = dolfin.TestFunction(vfs)
print "Defining variational forms…"
print_str(tab,"Defining variational forms…")
penalty = dolfin.Constant(penalty)
if (images_dimension == 2):
Idef = myFEniCS.ExprDefIm2(
......@@ -159,7 +173,7 @@ def fedic(
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)+" <= 0. Aborting."
print "B0_norm = " + str(B0_norm)
print_var(tab+1,"B0_norm",B0_norm)
# linear system
A = None
......@@ -167,17 +181,17 @@ def fedic(
A = dolfin.assemble(a, tensor=A, form_compiler_parameters={'quadrature_degree':images_quadrature})
B = None
print "Checking number of frames…"
print_str(tab,"Checking number of frames…")
n_frames = len(glob.glob(images_folder+"/"+images_basename+"_"+"?"*images_zfill+".vti"))
#n_frames = 2
assert (n_frames > 1), "n_frames = "+str(n_frames)+" <= 1. Aborting."
print "n_frames = " + str(n_frames)
print_var(tab+1,"n_frames",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)
print_var(tab+1,"images_k_ref",images_k_ref)
print "Printing initial solution…"
print_str(tab,"Printing initial solution…")
if not os.path.exists(working_folder):
os.mkdir(working_folder)
pvd_basename = working_folder+"/"+working_basename+"_"
......@@ -186,16 +200,17 @@ def fedic(
os.remove(vtu_filename)
file_pvd << (U, float(images_k_ref))
print_str(tab,"Looping over frames…")
n_iter_tot = 0
global_success = True
for forward_or_backward in ["forward", "backward"]: # not sure if this still works…
print "forward_or_backward = " + forward_or_backward
print_var(tab,"forward_or_backward",forward_or_backward)
if (forward_or_backward == "forward"):
k_frames = range(images_k_ref+1, n_frames, +1)
elif (forward_or_backward == "backward"):
k_frames = range(images_k_ref-1, -1, -1)
print "k_frames = " + str(k_frames)
print_var(tab,"k_frames",k_frames)
if (forward_or_backward == "backward"):
U.vector()[:] = 0.
......@@ -205,8 +220,9 @@ def fedic(
DIdef.init_image(filename=ref_image_filename)
DDIdef.init_image(filename=ref_image_filename)
tab += 1
for k_frame in k_frames:
print "k_frame = " + str(k_frame)
print_var(tab-1,"k_frame",k_frame)
if (print_iterations):
frame_basename = working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(images_zfill)
......@@ -223,11 +239,11 @@ def fedic(
# linear system: matrix
if (tangent_type.startswith("Idef_old")):
A = dolfin.assemble(a, tensor=A, form_compiler_parameters={'quadrature_degree':images_quadrature})
#print "A = " + str(A.array())
#print_var(tab,"A",A.array())
#A_norm = numpy.linalg.norm(A.array())
#print "A_norm = " + str(A_norm)
#print_var(tab,"A_norm",A_norm)
print "Loading image and image gradient…"
print_str(tab,"Loading image and image gradient…")
image_filename = images_folder+"/"+images_basename+"_"+str(k_frame).zfill(images_zfill)+".vti"
Idef.init_image( filename=image_filename)
DIdef.init_image(filename=image_filename)
......@@ -243,19 +259,20 @@ def fedic(
err_im = dI_norm/Iref_norm
file_dat_frame.write(" ".join([str(val) for val in [-1, None, None, None, None, None, None, dI_norm, err_im, None]])+"\n")
print "Running registration…"
print_str(tab,"Running registration…")
tab += 1
k_iter = 0
success = False
while (True):
print "k_iter = " + str(k_iter)
print_var(tab-1,"k_iter",k_iter)
n_iter_tot += 1
# linear system: matrix
if (tangent_type.startswith("Idef")):
A = dolfin.assemble(a, tensor=A, form_compiler_parameters={'quadrature_degree':images_quadrature})
#print "A = " + str(A.array())
#print_var(tab,"A",A.array())
#A_norm = numpy.linalg.norm(A.array())
#print "A_norm = " + str(A_norm)
#print_sci(tab,"A_norm",A_norm)
# linear system: residual
if (relax_type == "aitken"):
......@@ -263,19 +280,18 @@ def fedic(
B_old = B.copy()
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':images_quadrature})
#print "B = " + str(B.array())
#print_var(tab,"B",B.array())
# residual error
B_norm = numpy.linalg.norm(B.array())
#print "B_norm = " + str(B_norm)
#print_sci(tab,"B_norm",B_norm)
err_res = B_norm/B0_norm
print "err_res = " + str(err_res)
print_sci(tab,"err_res",err_res)
# linear system: solve
dolfin.solve(A, dU.vector(), B)
#print "dU = " + str(dU.vector().array())
#print_var(tab,"dU",dU.vector().array())
# relaxation
if (relax_type == "constant"):
......@@ -290,24 +306,25 @@ def fedic(
elif (k_iter > 1):
dB[:] = B[:] - B_old[:]
relax *= (-1.) * B_old.inner(dB) / dB.inner(dB)
print "relax = " + str(relax)
print_var(tab,"relax",relax)
elif (relax_type == "manual"):
B_relax = numpy.empty(relax_n_iter)
B_norm_relax = numpy.empty(relax_n_iter)
tab += 1
for relax_k in xrange(relax_n_iter):
#print "relax_k = " + str(relax_k)
print_var(tab-1,"relax_k",relax_k)
U.vector()[:] += (1./relax_n_iter) * dU.vector()[:]
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])
B_norm_relax[relax_k] = numpy.linalg.norm(B)
print_sci(tab,"B_norm_relax",B_norm_relax[relax_k])
U.vector()[:] -= dU.vector()[:]
#print "B_relax = " + str(B_relax)
#print_var(tab,"B_norm_relax",B_norm_relax)
if (print_iterations):
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)]))
open(iter_basename+".dat", "w").write("\n".join([" ".join([str(val) for val in [float(relax_k+1)/relax_n_iter, B_norm_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_k = numpy.argmin(B_norm_relax)
relax = float(relax_k+1)/relax_n_iter
print "relax = " + str(relax)
print_var(tab,"relax",relax)
else:
assert (0), "relax_type must be \"constant\", \"aitken\" or \"manual\". Aborting."
......@@ -316,30 +333,30 @@ def fedic(
U.vector()[:] += relax * dU.vector()[:]
if (print_iterations):
#print "U = " + str(U.vector().array())
#print_var(tab,"U",U.vector().array())
file_pvd_frame << (U, float(k_iter+1))
# displacement error
dU_norm = numpy.linalg.norm(dU.vector().array())
print "dU_norm = " + str(dU_norm)
print_sci(tab,"dU_norm",dU_norm)
U_norm = numpy.linalg.norm(U.vector().array())
assert (U_norm > 0.), "U_norm = "+str(U_norm)+" <= 0. Aborting."
print "U_norm = " + str(U_norm)
print_sci(tab,"U_norm",U_norm)
err_disp = dU_norm/U_norm
print "err_disp = " + str(err_disp)
print_sci(tab,"err_disp",err_disp)
# image error
if (k_iter > 0):
dI_norm_old = dI_norm
dI_norm = (dolfin.assemble((Idef-Iref)**2 * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})/mesh_volume)**(1./2)
#print "dI_norm = " + str(dI_norm)
#print_sci(tab,"dI_norm",dI_norm)
err_im = dI_norm/Iref_norm
print "err_im = " + str(err_im)
print_sci(tab,"err_im",err_im)
if (k_iter == 0):
err_im_rel = 1.
else:
err_im_rel = abs(dI_norm-dI_norm_old)/dI_norm_old
print "err_im_rel = " + str(err_im_rel)
print_sci(tab,"err_im_rel",err_im_rel)
if (print_iterations):
file_dat_frame.write(" ".join([str(val) for val in [k_iter, B_norm, err_res, relax, dU_norm, U_norm, err_disp, dI_norm, err_im, err_im_rel]])+"\n")
......@@ -355,17 +372,19 @@ def fedic(
if (err_res < tol_res):
success = True
if (success):
print "Nonlinear solver converged…"
print_str(tab,"Nonlinear solver converged…")
break
if (k_iter == n_iter_max-1):
global_success = False
print "Warning! Nonlinear solver failed to converge… (k_frame = "+str(k_frame)+")"
print_str(tab,"Warning! Nonlinear solver failed to converge… (k_frame = "+str(k_frame)+")")
break
# increment counter
k_iter += 1
tab -= 1
if (print_iterations):
os.remove(frame_basename+"_.pvd")
file_dat_frame.close()
......@@ -378,13 +397,15 @@ def fedic(
U_old.vector()[:] += DU_old.vector()[:]
DU_old.vector()[:] = 0.
print "Printing solution…"
print_str(tab,"Printing solution…")
file_pvd << (U, float(k_frame))
tab -= 1
if not (success) and not (continue_after_fail):
break
print "n_iter_tot = " + str(n_iter_tot)
print_var(tab,"n_iter_tot",n_iter_tot)
os.remove(pvd_basename+".pvd")
......
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