possibility to pick the image quadrarture degree from list

compute_quadrature_degree.py

@@ -22,18 +22,22 @@ def compute_quadrature_degree_from_points_count(
         image_filename,
         mesh_filebasename,
         mesh_ext="vtk",
+        compute_n_quad=False,
         deg_min=1,
         deg_max=20,
         verbose=1):
 
     image = myvtk.readImage(
         filename=image_filename,
-        verbose=verbose)
+        verbose=verbose-1)
     n_points = image.GetNumberOfPoints()
+    image_dimension = myvtk.getImageDimensionality(
+        image=image,
+        verbose=verbose-1)
 
     mesh = myvtk.readUGrid(
         filename=mesh_filebasename+"."+mesh_ext,
-        verbose=verbose)
+        verbose=verbose-1)
     n_cells = mesh.GetNumberOfCells()
 
     (cell_locator,
@@ -43,7 +47,7 @@ def compute_quadrature_degree_from_points_count(
      subId,
      dist) = myvtk.getCellLocator(
         mesh=mesh,
-        verbose=verbose)
+        verbose=verbose-1)
 
     point = numpy.empty(3)
     n_pixels_per_cell = numpy.zeros(n_cells)
@@ -56,26 +60,35 @@ def compute_quadrature_degree_from_points_count(
     n_pixels_per_cell_max = int(max(n_pixels_per_cell))
     n_pixels_per_cell_avg = int(sum(n_pixels_per_cell)/n_cells)
 
-    if (verbose):
-        #print "n_pixels_per_cell = "+str(n_pixels_per_cell)
-        #print "sum(n_pixels_per_cell) = "+str(sum(n_pixels_per_cell))
-        print "n_pixels_per_cell_max = "+str(n_pixels_per_cell_max)
-        print "n_pixels_per_cell_avg = "+str(n_pixels_per_cell_avg)
-
-    mesh = dolfin.Mesh(mesh_filebasename+"."+"xml")
-
-    for degree in xrange(deg_min,deg_max+1):
-        if (verbose): print "degree = "+str(degree)
-        n_quad = len(dolfin.FunctionSpace(
-            mesh,
-            dolfin.FiniteElement(
-                family="Quadrature",
-                cell=mesh.ufl_cell(),
-                degree=degree,
-                quad_scheme="default")).dofmap().dofs())/len(mesh.cells())
-        if (verbose): print "n_quad = "+str(n_quad)
-        #if (n_quad > n_pixels_per_cell_max): break
-        if (n_quad > n_pixels_per_cell_avg): break
+    #if (verbose): print "n_pixels_per_cell = "+str(n_pixels_per_cell)
+    #if (verbose): print "sum(n_pixels_per_cell) = "+str(sum(n_pixels_per_cell))
+    #if (verbose): print "n_pixels_per_cell_max = "+str(n_pixels_per_cell_max)
+    #if (verbose): print "n_pixels_per_cell_avg = "+str(n_pixels_per_cell_avg)
+
+    if (compute_n_quad):
+        mesh = dolfin.Mesh(mesh_filebasename+"."+"xml")
+        #n_quads = []
+        for degree in xrange(deg_min,deg_max+1):
+            if (verbose): print "degree = "+str(degree)
+            n_quad = len(dolfin.FunctionSpace(
+                mesh,
+                dolfin.FiniteElement(
+                    family="Quadrature",
+                    cell=mesh.ufl_cell(),
+                    degree=degree,
+                    quad_scheme="default")).dofmap().dofs())/len(mesh.cells())
+            if (verbose): print "n_quad = "+str(n_quad)
+            #if (n_quad > n_pixels_per_cell_max): break
+            if (n_quad > n_pixels_per_cell_avg): break
+            #n_quads.append(n_quad)
+            #print n_quads
+    else:
+        if (image_dimension == 2):
+            n_quads = [1, 3, 6, 6, 7, 12, 16, 25, 25, 36, 36, 49, 49, 64, 64, 81, 81, 100, 100, 121]
+        elif (image_dimension == 3):
+            n_quads = [1, 4, 5, 14, 15, 24, 64, 125, 125, 216, 216, 343, 343, 512, 512, 729, 729]
+        #degree = numpy.searchsorted(n_quads, n_pixels_per_cell_max)+1
+        degree = numpy.searchsorted(n_quads, n_pixels_per_cell_avg)+1
 
     return degree
 
@@ -92,10 +105,10 @@ def compute_quadrature_degree_from_integral(
 
     image = myvtk.readImage(
         filename=image_filename,
-        verbose=verbose)
+        verbose=verbose-1)
     image_dimension = myvtk.getImageDimensionality(
         image=image,
-        verbose=verbose)
+        verbose=verbose-1)
     if (verbose): print "image_dimension = " + str(image_dimension)
 
     dX = dolfin.dx(mesh)