Various fixes

__init__.py

 #this file was generated by __init__.py.py
 
-from fedic import *
-from plot_strains_vs_radius import *
-from generate_undersampled_images import *
-from plot_strains import *
-from image_expressions_cpp import *
-from compute_strains import *
 from compute_warped_mesh import *
-from generate_images import *
-from plot_displacement_error import *
-from compute_displacements import *
-from compute_unwarped_images import *
 from image_expressions_py import *
+from compute_warped_images import *
+from compute_unwarped_images import *
+from fedic import *
+from generate_images import *
 from plot_binned_strains_vs_radius import *
 from plot_regional_strains import *
+from generate_undersampled_images import *
+from compute_strains import *
+from compute_displacements import *
 from compute_displacement_error import *
-from compute_warped_images import *
 from compute_quadrature_degree import *
+from image_expressions_cpp import *
+from plot_strains import *
+from plot_displacement_error import *
+from plot_strains_vs_radius import *

fedic.py

@@ -182,9 +182,9 @@ def fedic(
         cell=mesh.ufl_cell(),
         degree=images_quadrature,
         quad_scheme="default")
-    te._quad_scheme = "default"                       # should not be needed
-    for k in xrange(mesh_dimension**2):             # should not be needed
-        te.sub_elements()[k]._quad_scheme = "default" # should not be needed
+    te._quad_scheme = "default"              # should not be needed
+    for sub_element in te.sub_elements():    # should not be needed
+        sub_element._quad_scheme = "default" # should not be needed
     if (images_expressions_type == "cpp"):
         Iref = dolfin.Expression(
             cppcode=ddic.get_ExprIm_cpp(
@@ -301,12 +301,12 @@ def fedic(
         mypy.print_str("Defining regularization energy…",tab)
         E     = dolfin.Constant(1.0)
         nu    = dolfin.Constant(regul_poisson)
-        kappa = E/3/(1-2*nu)                          # = E/3 if nu = 0
-        lmbda = E*nu/(1+nu)/(1-(mesh_dimension-1)*nu) # = 0   if nu = 0
-        mu    = E/2/(1+nu)                            # = E/2 if nu = 0
-        C1    = mu/2                                  # = E/4 if nu = 0
-        C2    = mu/2                                  # = E/4 if nu = 0
-        D1    = kappa/2                               # = E/6 if nu = 0
+        lmbda = E*nu/(1+nu)/(1-2*nu)                        # = 0   if nu = 0
+        mu    = E/2/(1+nu)                                  # = E/2 if nu = 0
+        kappa = (mesh_dimension*lmbda+2*mu)/mesh_dimension  # = E/3 (3D) or E/2 (2D) if nu = 0
+        C1    = mu/2                                        # = E/4 if nu = 0
+        C2    = mu/2                                        # = E/4 if nu = 0
+        D1    = kappa/2                                     # = E/6 (3D) or E/4 (2D) if nu = 0
 
         if (regul_model == "linear"): # <- super bad
             e = dolfin.sym(dolfin.grad(U))
@@ -384,7 +384,7 @@ def fedic(
         regul_quadrature = None
 
         form_compiler_parameters_for_regul = {}
-        form_compiler_parameters_for_regul["representation"] = "uflacs"
+        form_compiler_parameters_for_regul["representation"] = "uflacs" # MG20180327: Is that needed?
         if (regul_quadrature is not None):
             form_compiler_parameters_for_regul["quadrature_degree"] = regul_quadrature
 

generate_images.py

@@ -509,12 +509,12 @@ def generateImages(
     else:
         assert (0), "n_dim must be \"1\", \"2\" or \"3\". Aborting."
 
-    spacing = numpy.array(images["L"])/numpy.array(images["n_voxels"])
+    #spacing = numpy.array(images["L"])/numpy.array(images["n_voxels"])
     if   (images["n_dim"] == 1):
         spacing = numpy.array([images["L"][0]/images["n_voxels"][0], 1., 1.])
     elif (images["n_dim"] == 2):
         spacing = numpy.array([images["L"][0]/images["n_voxels"][0], images["L"][1]/images["n_voxels"][1], 1.])
-    elif (images["n_dim"] == 2):
+    elif (images["n_dim"] == 3):
         spacing = numpy.array([images["L"][0]/images["n_voxels"][0], images["L"][1]/images["n_voxels"][1], images["L"][2]/images["n_voxels"][2]])
     vtk_image.SetSpacing(spacing)
 
@@ -544,16 +544,16 @@ def generateImages(
     if not os.path.exists(images["folder"]):
         os.mkdir(images["folder"])
 
-    x0   = numpy.empty(3)
-    x    = numpy.empty(3)
-    X    = numpy.empty(3)
+    x0 = numpy.empty(3)
+    x  = numpy.empty(3)
+    X  = numpy.empty(3)
     if (generate_image_gradient):
         F    = numpy.empty((3,3))
         Finv = numpy.empty((3,3))
     else:
         F    = None
         Finv = None
-    dx   = spacing[0:images["n_dim"]]/images["n_integration"][0:images["n_dim"]]
+    dx = spacing[0:images["n_dim"]]/images["n_integration"][0:images["n_dim"]]
     global_min = float("+Inf")
     global_max = float("-Inf")
     I = numpy.empty(1)
@@ -627,7 +627,7 @@ def generateImages(
             if (I[0] > global_max): global_max = I[0]
         myvtk.writeImage(
             image=vtk_image,
-            filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+".vti",
+            filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+".vtk",
             verbose=verbose-1)
 
     if (images["data_type"] in ("float")):