diff --git a/api/GenApi.py b/api/GenApi.py
index a0400d1e06eb2c81dacecc4871b5f16c970d4daa..b3cd86168cc53d44f7b0e8c6d19bf3120b148283 100644
--- a/api/GenApi.py
+++ b/api/GenApi.py
@@ -11,9 +11,10 @@ import string
 import os
 import re
 
+
 class arg:
-    def __init__(self, name, value, python_value, julia_value,
-                 cpp_type, c_type, out):
+    def __init__(self, name, value, python_value, julia_value, cpp_type,
+                 c_type, out):
         self.name = name
         self.value = value
         self.out = out
@@ -35,59 +36,68 @@ class arg:
         self.julia_pre = ""
         self.julia_post = ""
         self.julia_return = name
-        self.texi = name + ((" = " + self.python_value) if self.python_value else "")
+        self.texi = name + (
+            (" = " + self.python_value) if self.python_value else "")
+
 
 # input types
 
+
 def ibool(name, value=None, python_value=None, julia_value=None):
-    a = arg(name, value, python_value, julia_value,
-            "const bool", "const int", False)
+    a = arg(name, value, python_value, julia_value, "const bool", "const int",
+            False)
     a.python_arg = "c_int(bool(" + name + "))"
     a.cwrap_arg = "(int)" + name
     a.julia_ctype = "Cint"
     return a
 
+
 def iint(name, value=None, python_value=None, julia_value=None):
-    a = arg(name, value, python_value, julia_value,
-            "const int", "const int", False)
+    a = arg(name, value, python_value, julia_value, "const int", "const int",
+            False)
     a.python_arg = "c_int(" + name + ")"
     a.julia_ctype = "Cint"
     return a
 
+
 def isize(name, value=None, python_value=None, julia_value=None):
-    a = arg(name, value, python_value, julia_value,
-            "const std::size_t", "const size_t", False)
+    a = arg(name, value, python_value, julia_value, "const std::size_t",
+            "const size_t", False)
     a.python_arg = "c_size_t(" + name + ")"
     a.julia_ctype = "Csize_t"
     return a
 
+
 def idouble(name, value=None, python_value=None, julia_value=None):
-    a = arg(name, value, python_value, julia_value,
-            "const double", "const double", False)
+    a = arg(name, value, python_value, julia_value, "const double",
+            "const double", False)
     a.python_arg = "c_double(" + name + ")"
     a.julia_ctype = "Cdouble"
     return a
 
+
 def istring(name, value=None, python_value=None, julia_value=None):
-    a = arg(name, value, python_value, julia_value,
-            "const std::string &", "const char *", False)
+    a = arg(name, value, python_value, julia_value, "const std::string &",
+            "const char *", False)
     a.python_arg = "c_char_p(" + name + ".encode())"
     a.cwrap_arg = name + ".c_str()"
     a.julia_ctype = "Ptr{Cchar}"
     return a
 
+
 def ivoidstar(name, value=None, python_value=None, julia_value=None):
-    a = arg(name, value, python_value, julia_value,
-            "const void *", "const void *", False)
+    a = arg(name, value, python_value, julia_value, "const void *",
+            "const void *", False)
     a.python_arg = "c_void_p(" + name + ")"
     a.julia_ctype = "Ptr{Cvoid}"
     return a
 
+
 def ivectorint(name, value=None, python_value=None, julia_value=None):
     if julia_value == "[]":
         julia_value = "Cint[]"
-    a = arg(name, value, python_value, julia_value,
-            "const std::vector<int> &", "const int *", False)
+    a = arg(name, value, python_value, julia_value, "const std::vector<int> &",
+            "const int *", False)
     api_name = "api_" + name + "_"
     api_name_n = "api_" + name + "_n_"
     a.c_pre = ("    std::vector<int> " + api_name + "(" + name + ", " + name +
@@ -95,7 +105,8 @@ def ivectorint(name, value=None, python_value=None, julia_value=None):
     a.c_arg = api_name
     a.c = "int * " + name + ", size_t " + name + "_n"
     a.cwrap_pre = ("int *" + api_name + "; size_t " + api_name_n + "; " +
-                   "vector2ptr(" + name + ", &" + api_name + ", &" + api_name_n + ");\n")
+                   "vector2ptr(" + name + ", &" + api_name + ", &" +
+                   api_name_n + ");\n")
     a.cwrap_arg = api_name + ", " + api_name_n
     a.cwrap_post = ns + "Free(" + api_name + ");\n"
     a.python_pre = api_name + ", " + api_name_n + " = _ivectorint(" + name + ")"
@@ -104,6 +115,7 @@ def ivectorint(name, value=None, python_value=None, julia_value=None):
     a.julia_arg = "convert(Vector{Cint}, " + name + "), length(" + name + ")"
     return a
 
+
 def ivectorsize(name, value=None, python_value=None, julia_value=None):
     if julia_value == "[]":
         julia_value = "Csize_t[]"
@@ -111,12 +123,13 @@ def ivectorsize(name, value=None, python_value=None, julia_value=None):
             "const std::vector<std::size_t> &", "const size_t *", False)
     api_name = "api_" + name + "_"
     api_name_n = "api_" + name + "_n_"
-    a.c_pre = ("    std::vector<std::size_t> " + api_name + "(" + name + ", " + name +
-               " + " + name + "_n);\n")
+    a.c_pre = ("    std::vector<std::size_t> " + api_name + "(" + name + ", " +
+               name + " + " + name + "_n);\n")
     a.c_arg = api_name
     a.c = "size_t * " + name + ", size_t " + name + "_n"
     a.cwrap_pre = ("size_t *" + api_name + "; size_t " + api_name_n + "; " +
-                   "vector2ptr(" + name + ", &" + api_name + ", &" + api_name_n + ");\n")
+                   "vector2ptr(" + name + ", &" + api_name + ", &" +
+                   api_name_n + ");\n")
     a.cwrap_arg = api_name + ", " + api_name_n
     a.cwrap_post = ns + "Free(" + api_name + ");\n"
     a.python_pre = api_name + ", " + api_name_n + " = _ivectorsize(" + name + ")"
@@ -125,6 +138,7 @@ def ivectorsize(name, value=None, python_value=None, julia_value=None):
     a.julia_arg = "convert(Vector{Csize_t}, " + name + "), length(" + name + ")"
     return a
 
+
 def ivectordouble(name, value=None, python_value=None, julia_value=None):
     if julia_value == "[]":
         julia_value = "Cdouble[]"
@@ -137,7 +151,8 @@ def ivectordouble(name, value=None, python_value=None, julia_value=None):
     a.c_arg = api_name
     a.c = "double * " + name + ", size_t " + name + "_n"
     a.cwrap_pre = ("double *" + api_name + "; size_t " + api_name_n + "; " +
-                   "vector2ptr(" + name + ", &" + api_name + ", &" + api_name_n + ");\n")
+                   "vector2ptr(" + name + ", &" + api_name + ", &" +
+                   api_name_n + ");\n")
     a.cwrap_arg = api_name + ", " + api_name_n
     a.cwrap_post = ns + "Free(" + api_name + ");\n"
     a.python_pre = api_name + ", " + api_name_n + " = _ivectordouble(" + name + ")"
@@ -146,6 +161,7 @@ def ivectordouble(name, value=None, python_value=None, julia_value=None):
     a.julia_arg = "convert(Vector{Cdouble}, " + name + "), length(" + name + ")"
     return a
 
+
 def ivectorstring(name, value=None, python_value=None, julia_value=None):
     a = arg(name, value, python_value, julia_value,
             "const std::vector<std::string> &", "char **", False)
@@ -156,17 +172,19 @@ def ivectorstring(name, value=None, python_value=None, julia_value=None):
     a.c_arg = api_name
     a.c = "char ** " + name + ", size_t " + name + "_n"
     a.cwrap_pre = ("char **" + api_name + "; size_t " + api_name_n + "; " +
-                   "vectorstring2charptrptr(" + name + ", &" + api_name + ", &" + api_name_n + ");\n")
+                   "vectorstring2charptrptr(" + name + ", &" + api_name +
+                   ", &" + api_name_n + ");\n")
     a.cwrap_arg = api_name + ", " + api_name_n
-    a.cwrap_post = ("for(size_t i = 0; i < " + api_name_n + "; ++i){ " +
-                    ns + "Free(" + api_name + "[i]); } " +
-                    ns + "Free(" + api_name + ");\n")
+    a.cwrap_post = ("for(size_t i = 0; i < " + api_name_n + "; ++i){ " + ns +
+                    "Free(" + api_name + "[i]); } " + ns + "Free(" + api_name +
+                    ");\n")
     a.python_pre = api_name + ", " + api_name_n + " = _ivectorstring(" + name + ")"
     a.python_arg = api_name + ", " + api_name_n
     a.julia_ctype = "Ptr{Ptr{Cchar}}, Csize_t"
     a.julia_arg = name + ", length(" + name + ")"
     return a
 
+
 def ivectorpair(name, value=None, python_value=None, julia_value=None):
     if julia_value == "[]":
         julia_value = "Tuple{Cint,Cint}[]"
@@ -174,11 +192,11 @@ def ivectorpair(name, value=None, python_value=None, julia_value=None):
             "const " + ns + "::vectorpair &", "const int *", False)
     api_name = "api_" + name + "_"
     api_name_n = "api_" + name + "_n_"
-    a.c_pre = ("    " + ns + "::vectorpair " + api_name + "(" + name + "_n/2);\n" +
-               "    for(size_t i = 0; i < " + name + "_n/2; ++i){\n" +
-               "      " + api_name + "[i].first = " + name + "[i * 2 + 0];\n" +
-               "      " + api_name + "[i].second = " + name + "[i * 2 + 1];\n" +
-               "    }\n")
+    a.c_pre = ("    " + ns + "::vectorpair " + api_name + "(" + name +
+               "_n/2);\n" + "    for(size_t i = 0; i < " + name +
+               "_n/2; ++i){\n" + "      " + api_name + "[i].first = " + name +
+               "[i * 2 + 0];\n" + "      " + api_name + "[i].second = " +
+               name + "[i * 2 + 1];\n" + "    }\n")
     a.c_arg = api_name
     a.c = "int * " + name + ", size_t " + name + "_n"
     a.cwrap_pre = ("int *" + api_name + "; size_t " + api_name_n + "; " +
@@ -189,11 +207,12 @@ def ivectorpair(name, value=None, python_value=None, julia_value=None):
     a.python_pre = api_name + ", " + api_name_n + " = _ivectorpair(" + name + ")"
     a.python_arg = api_name + ", " + api_name_n
     a.julia_ctype = "Ptr{Cint}, Csize_t"
-    a.julia_pre = (api_name + " = collect(Cint, Iterators.flatten(" + name + "))\n    " +
-                   api_name_n + " = length(" + api_name + ")")
+    a.julia_pre = (api_name + " = collect(Cint, Iterators.flatten(" + name +
+                   "))\n    " + api_name_n + " = length(" + api_name + ")")
     a.julia_arg = (api_name + ", " + api_name_n)
     return a
 
+
 def ivectorvectorint(name, value=None, python_value=None, julia_value=None):
     if julia_value == "[]":
         julia_value = "Vector{Cint}[]"
@@ -202,108 +221,117 @@ def ivectorvectorint(name, value=None, python_value=None, julia_value=None):
     api_name = "api_" + name + "_"
     api_name_n = "api_" + name + "_n_"
     api_name_nn = "api_" + name + "_nn_"
-    a.c_pre = ("    std::vector<std::vector<int> > " + api_name +
-               "(" + name + "_nn);\n" +
-               "    for(size_t i = 0; i < " + name + "_nn; ++i)\n" +
-               "      " + api_name + "[i] = std::vector<int>(" + name + "[i], " +
-               name + "[i] + " + name + "_n[i]);\n")
+    a.c_pre = ("    std::vector<std::vector<int> > " + api_name + "(" + name +
+               "_nn);\n" + "    for(size_t i = 0; i < " + name +
+               "_nn; ++i)\n" + "      " + api_name +
+               "[i] = std::vector<int>(" + name + "[i], " + name + "[i] + " +
+               name + "_n[i]);\n")
     a.c_arg = api_name
     a.c = ("const int ** " + name + ", const size_t * " + name + "_n, " +
            "size_t " + name + "_nn")
     a.cwrap_pre = ("int **" + api_name + "; size_t *" + api_name_n + ", " +
                    api_name_nn + "; " + "vectorvector2ptrptr(" + name + ", &" +
-                   api_name + ", &" + api_name_n + ", &" + api_name_nn + ");\n")
+                   api_name + ", &" + api_name_n + ", &" + api_name_nn +
+                   ");\n")
     a.cwrap_arg = "(const int **)" + api_name + ", " + api_name_n + ", " + api_name_nn
-    a.cwrap_post = ("for(size_t i = 0; i < " + api_name_nn + "; ++i){ " +
-                    ns + "Free(" + api_name + "[i]); } " +
-                    ns + "Free(" + api_name + "); " + ns + "Free(" + api_name_n + ");\n")
-    a.python_pre = (api_name + ", " + api_name_n + ", " +
-                    api_name_nn + " = _ivectorvectorint(" + name + ")")
+    a.cwrap_post = ("for(size_t i = 0; i < " + api_name_nn + "; ++i){ " + ns +
+                    "Free(" + api_name + "[i]); } " + ns + "Free(" + api_name +
+                    "); " + ns + "Free(" + api_name_n + ");\n")
+    a.python_pre = (api_name + ", " + api_name_n + ", " + api_name_nn +
+                    " = _ivectorvectorint(" + name + ")")
     a.python_arg = api_name + ", " + api_name_n + ", " + api_name_nn
     a.julia_ctype = "Ptr{Ptr{Cint}}, Ptr{Csize_t}, Csize_t"
-    a.julia_pre = (api_name_n + " = [ length(" + name + "[i]) for i in 1:length(" +
-                   name + ") ]")
-    a.julia_arg = ("convert(Vector{Vector{Cint}}," + name + "), " + api_name_n +
-                   ", length(" + name + ")")
+    a.julia_pre = (api_name_n + " = [ length(" + name +
+                   "[i]) for i in 1:length(" + name + ") ]")
+    a.julia_arg = ("convert(Vector{Vector{Cint}}," + name + "), " +
+                   api_name_n + ", length(" + name + ")")
     return a
 
+
 def ivectorvectorsize(name, value=None, python_value=None, julia_value=None):
     if julia_value == "[]":
         julia_value = "Vector{Csize_t}[]"
     a = arg(name, value, python_value, julia_value,
-            "const std::vector<std::vector<std::size_t> > &", "const size_t **", False)
+            "const std::vector<std::vector<std::size_t> > &",
+            "const size_t **", False)
     api_name = "api_" + name + "_"
     api_name_n = "api_" + name + "_n_"
     api_name_nn = "api_" + name + "_nn_"
-    a.c_pre = ("    std::vector<std::vector<std::size_t> > " + api_name +
-               "(" + name + "_nn);\n" +
-               "    for(size_t i = 0; i < " + name + "_nn; ++i)\n" +
-               "      " + api_name + "[i] = std::vector<std::size_t>(" + name + "[i], " +
-               name + "[i] + " + name + "_n[i]);\n")
+    a.c_pre = ("    std::vector<std::vector<std::size_t> > " + api_name + "(" +
+               name + "_nn);\n" + "    for(size_t i = 0; i < " + name +
+               "_nn; ++i)\n" + "      " + api_name +
+               "[i] = std::vector<std::size_t>(" + name + "[i], " + name +
+               "[i] + " + name + "_n[i]);\n")
     a.c_arg = api_name
     a.c = ("const size_t ** " + name + ", const size_t * " + name + "_n, " +
            "size_t " + name + "_nn")
     a.cwrap_pre = ("size_t **" + api_name + "; size_t *" + api_name_n + ", " +
                    api_name_nn + "; " + "vectorvector2ptrptr(" + name + ", &" +
-                   api_name + ", &" + api_name_n + ", &" + api_name_nn + ");\n")
+                   api_name + ", &" + api_name_n + ", &" + api_name_nn +
+                   ");\n")
     a.cwrap_arg = "(const size_t **)" + api_name + ", " + api_name_n + ", " + api_name_nn
-    a.cwrap_post = ("for(size_t i = 0; i < " + api_name_nn + "; ++i){ " +
-                    ns + "Free(" + api_name + "[i]); } " +
-                    ns + "Free(" + api_name + "); " + ns + "Free(" + api_name_n + ");\n")
-    a.python_pre = (api_name + ", " + api_name_n + ", " +
-                    api_name_nn + " = _ivectorvectorsize(" + name + ")")
+    a.cwrap_post = ("for(size_t i = 0; i < " + api_name_nn + "; ++i){ " + ns +
+                    "Free(" + api_name + "[i]); } " + ns + "Free(" + api_name +
+                    "); " + ns + "Free(" + api_name_n + ");\n")
+    a.python_pre = (api_name + ", " + api_name_n + ", " + api_name_nn +
+                    " = _ivectorvectorsize(" + name + ")")
     a.python_arg = api_name + ", " + api_name_n + ", " + api_name_nn
     a.julia_ctype = "Ptr{Ptr{Csize_t}}, Ptr{Csize_t}, Csize_t"
-    a.julia_pre = (api_name_n + " = [ length(" + name + "[i]) for i in 1:length(" +
-                   name + ") ]")
-    a.julia_arg = ("convert(Vector{Vector{Csize_t}}," + name + "), " + api_name_n +
-                   ", length(" + name + ")")
+    a.julia_pre = (api_name_n + " = [ length(" + name +
+                   "[i]) for i in 1:length(" + name + ") ]")
+    a.julia_arg = ("convert(Vector{Vector{Csize_t}}," + name + "), " +
+                   api_name_n + ", length(" + name + ")")
     return a
 
+
 def ivectorvectordouble(name, value=None, python_value=None, julia_value=None):
     if julia_value == "[]":
         julia_value = "Vector{Cdouble}[]"
     a = arg(name, value, python_value, julia_value,
-            "const std::vector<std::vector<double> > &", "const double**", False)
+            "const std::vector<std::vector<double> > &", "const double**",
+            False)
     api_name = "api_" + name + "_"
     api_name_n = "api_" + name + "_n_"
     api_name_nn = "api_" + name + "_nn_"
-    a.c_pre = ("    std::vector<std::vector<double> > " + api_name +
-               "(" + name + "_nn);\n" +
-               "    for(size_t i = 0; i < " + name + "_nn; ++i)\n" +
-               "      " + api_name + "[i] = std::vector<double>(" + name + "[i], " +
-               name + "[i] + " + name + "_n[i]);\n")
+    a.c_pre = ("    std::vector<std::vector<double> > " + api_name + "(" +
+               name + "_nn);\n" + "    for(size_t i = 0; i < " + name +
+               "_nn; ++i)\n" + "      " + api_name +
+               "[i] = std::vector<double>(" + name + "[i], " + name +
+               "[i] + " + name + "_n[i]);\n")
     a.c_arg = api_name
     a.c = ("const double ** " + name + ", const size_t * " + name + "_n, " +
            "size_t " + name + "_nn")
     a.cwrap_pre = ("double **" + api_name + "; size_t *" + api_name_n + ", " +
-                   api_name_nn + "; " +
-                   "vectorvector2ptrptr(" + name + ", &" + api_name + ", &" +
-                   api_name_n + ", &" + api_name_nn + ");\n")
+                   api_name_nn + "; " + "vectorvector2ptrptr(" + name + ", &" +
+                   api_name + ", &" + api_name_n + ", &" + api_name_nn +
+                   ");\n")
     a.cwrap_arg = "(const double **)" + api_name + ", " + api_name_n + ", " + api_name_nn
-    a.cwrap_post = ("for(size_t i = 0; i < " + api_name_nn + "; ++i){ " +
-                    ns + "Free(" + api_name + "[i]); } " +
-                    ns + "Free(" + api_name + "); " + ns + "Free(" + api_name_n + ");\n")
+    a.cwrap_post = ("for(size_t i = 0; i < " + api_name_nn + "; ++i){ " + ns +
+                    "Free(" + api_name + "[i]); } " + ns + "Free(" + api_name +
+                    "); " + ns + "Free(" + api_name_n + ");\n")
     a.python_pre = (api_name + ", " + api_name_n + ", " + api_name_nn +
                     " = _ivectorvectordouble(" + name + ")")
     a.python_arg = api_name + ", " + api_name_n + ", " + api_name_nn
     a.julia_ctype = "Ptr{Ptr{Cdouble}}, Ptr{Csize_t}, Csize_t"
-    a.julia_pre = (api_name_n + " = [ length(" + name + "[i]) for i in 1:length(" +
-                   name + ") ]")
-    a.julia_arg = ("convert(Vector{Vector{Cdouble}}," + name + "), " + api_name_n +
-                   ", length(" + name + ")")
+    a.julia_pre = (api_name_n + " = [ length(" + name +
+                   "[i]) for i in 1:length(" + name + ") ]")
+    a.julia_arg = ("convert(Vector{Vector{Cdouble}}," + name + "), " +
+                   api_name_n + ", length(" + name + ")")
     return a
 
+
 # output types
 
+
 class oint(arg):
     rcpp_type = "int"
     rc_type = "int"
     rtexi_type = "integer value"
     rjulia_type = "Cint"
+
     def __init__(self, name, value=None, python_value=None, julia_value=None):
-        arg.__init__(self, name, value, python_value, julia_value,
-                     "int &", "int *", True)
+        arg.__init__(self, name, value, python_value, julia_value, "int &",
+                     "int *", True)
         api_name = "api_" + name + "_"
         self.c_arg = "*" + name
         self.cwrap_arg = "&" + name
@@ -315,11 +343,13 @@ class oint(arg):
         self.julia_arg = api_name
         self.julia_return = api_name + "[]"
 
+
 class osize(arg):
     rcpp_type = "std::size_t"
     rc_type = "size_t"
     rtexi_type = "size value"
     rjulia_type = "Csize_t"
+
     def __init__(self, name, value=None, python_value=None, julia_value=None):
         arg.__init__(self, name, value, python_value, julia_value,
                      "std::size_t &", "size_t *", True)
@@ -334,14 +364,16 @@ class osize(arg):
         self.julia_arg = api_name
         self.julia_return = api_name + "[]"
 
+
 class odouble(arg):
     rcpp_type = "double"
     rc_type = "double"
     rtexi_type = "floating point value"
     rjulia_type = "Cdouble"
+
     def __init__(self, name, value=None, python_value=None, julia_value=None):
-        arg.__init__(self, name, value, python_value, julia_value,
-                     "double &", "double *", True)
+        arg.__init__(self, name, value, python_value, julia_value, "double &",
+                     "double *", True)
         api_name = "api_" + name + "_"
         self.c_arg = "*" + name
         self.cwrap_arg = "&" + name
@@ -353,17 +385,18 @@ class odouble(arg):
         self.julia_arg = api_name
         self.julia_return = api_name + "[]"
 
+
 def ostring(name, value=None, python_value=None, julia_value=None):
-    a = arg(name, value, python_value, julia_value,
-            "std::string &", "char **", True)
+    a = arg(name, value, python_value, julia_value, "std::string &", "char **",
+            True)
     api_name = "api_" + name + "_"
     a.c_pre = "    std::string " + api_name + ";\n"
     a.c_arg = api_name
     a.c_post = "    *" + name + " = strdup(" + api_name + ".c_str());\n"
     a.cwrap_pre = "char *" + api_name + ";\n"
     a.cwrap_arg = "&" + api_name
-    a.cwrap_post = (name + " = std::string(" + api_name + "); " +
-                    ns + "Free(" + api_name + ");\n")
+    a.cwrap_post = (name + " = std::string(" + api_name + "); " + ns +
+                    "Free(" + api_name + ");\n")
     a.python_pre = api_name + " = c_char_p()"
     a.python_arg = "byref(" + api_name + ")"
     a.python_return = "_ostring(" + api_name + ")"
@@ -373,9 +406,10 @@ def ostring(name, value=None, python_value=None, julia_value=None):
     a.julia_post = name + " = unsafe_string(" + api_name + "[])"
     return a
 
+
 def ovectorint(name, value=None, python_value=None, julia_value=None):
-    a = arg(name, value, python_value, julia_value,
-            "std::vector<int> &", "int **", True)
+    a = arg(name, value, python_value, julia_value, "std::vector<int> &",
+            "int **", True)
     api_name = "api_" + name + "_"
     api_name_n = api_name + "n_"
     a.c_pre = "    std::vector<int> " + api_name + ";\n"
@@ -390,13 +424,14 @@ def ovectorint(name, value=None, python_value=None, julia_value=None):
     a.python_arg = "byref(" + api_name + "), byref(" + api_name_n + ")"
     a.python_return = "_ovectorint(" + api_name + ", " + api_name_n + ".value)"
     a.julia_ctype = "Ptr{Ptr{Cint}}, Ptr{Csize_t}"
-    a.julia_pre = (api_name + " = Ref{Ptr{Cint}}()\n    " +
-                   api_name_n + " = Ref{Csize_t}()")
+    a.julia_pre = (api_name + " = Ref{Ptr{Cint}}()\n    " + api_name_n +
+                   " = Ref{Csize_t}()")
     a.julia_arg = api_name + ", " + api_name_n
     a.julia_post = (name + " = unsafe_wrap(Array, " + api_name + "[], " +
                     api_name_n + "[], own=true)")
     return a
 
+
 def ovectorsize(name, value=None, python_value=None, julia_value=None):
     a = arg(name, value, python_value, julia_value,
             "std::vector<std::size_t> &", "size_t **", True)
@@ -414,16 +449,17 @@ def ovectorsize(name, value=None, python_value=None, julia_value=None):
     a.python_arg = "byref(" + api_name + "), byref(" + api_name_n + ")"
     a.python_return = "_ovectorsize(" + api_name + ", " + api_name_n + ".value)"
     a.julia_ctype = "Ptr{Ptr{Csize_t}}, Ptr{Csize_t}"
-    a.julia_pre = (api_name + " = Ref{Ptr{Csize_t}}()\n    " +
-                   api_name_n + " = Ref{Csize_t}()")
+    a.julia_pre = (api_name + " = Ref{Ptr{Csize_t}}()\n    " + api_name_n +
+                   " = Ref{Csize_t}()")
     a.julia_arg = api_name + ", " + api_name_n
     a.julia_post = (name + " = unsafe_wrap(Array, " + api_name + "[], " +
                     api_name_n + "[], own=true)")
     return a
 
+
 def ovectordouble(name, value=None, python_value=None, julia_value=None):
-    a = arg(name, value, python_value, julia_value,
-            "std::vector<double> &", "double *", True)
+    a = arg(name, value, python_value, julia_value, "std::vector<double> &",
+            "double *", True)
     api_name = "api_" + name + "_"
     api_name_n = api_name + "n_"
     a.c_pre = "    std::vector<double> " + api_name + ";\n"
@@ -438,13 +474,14 @@ def ovectordouble(name, value=None, python_value=None, julia_value=None):
     a.python_arg = "byref(" + api_name + "), byref(" + api_name_n + ")"
     a.python_return = "_ovectordouble(" + api_name + ", " + api_name_n + ".value)"
     a.julia_ctype = "Ptr{Ptr{Cdouble}}, Ptr{Csize_t}"
-    a.julia_pre = (api_name + " = Ref{Ptr{Cdouble}}()\n    " +
-                   api_name_n + " = Ref{Csize_t}()")
+    a.julia_pre = (api_name + " = Ref{Ptr{Cdouble}}()\n    " + api_name_n +
+                   " = Ref{Csize_t}()")
     a.julia_arg = api_name + ", " + api_name_n
     a.julia_post = (name + " = unsafe_wrap(Array, " + api_name + "[], " +
                     api_name_n + "[], own=true)")
     return a
 
+
 def ovectorstring(name, value=None, python_value=None, julia_value=None):
     a = arg(name, value, python_value, julia_value,
             "std::vector<std::string> &", "char **", True)
@@ -452,32 +489,32 @@ def ovectorstring(name, value=None, python_value=None, julia_value=None):
     api_name_n = api_name + "n_"
     a.c_pre = "    std::vector<std::string> " + api_name + ";\n"
     a.c_arg = api_name
-    a.c_post = ("    vectorstring2charptrptr(" + api_name + ", " + name + ", " +
-                name + "_n);\n")
+    a.c_post = ("    vectorstring2charptrptr(" + api_name + ", " + name +
+                ", " + name + "_n);\n")
     a.c = "char *** " + name + ", size_t * " + name + "_n"
     a.cwrap_pre = "char **" + api_name + "; size_t " + api_name_n + ";\n"
     a.cwrap_arg = "&" + api_name + ", " + "&" + api_name_n
     a.cwrap_post = (name + ".resize(" + api_name_n + "); " +
-                    "for(size_t i = 0; i < " + api_name_n + "; ++i){ " +
-                    name + "[i] = std::string(" + api_name + "[i]); " +
-                    ns + "Free(" + api_name + "[i]); } " +
-                    ns + "Free(" + api_name + ");\n")
+                    "for(size_t i = 0; i < " + api_name_n + "; ++i){ " + name +
+                    "[i] = std::string(" + api_name + "[i]); " + ns + "Free(" +
+                    api_name + "[i]); } " + ns + "Free(" + api_name + ");\n")
     a.python_pre = api_name + ", " + api_name_n + " = POINTER(POINTER(c_char))(), c_size_t()"
     a.python_arg = "byref(" + api_name + "), byref(" + api_name_n + ")"
     a.python_return = "_ovectorstring(" + api_name + ", " + api_name_n + ".value)"
     a.julia_ctype = "Ptr{Ptr{Ptr{Cchar}}}, Ptr{Csize_t}"
-    a.julia_pre = (api_name + " = Ref{Ptr{Ptr{Cchar}}}()\n    " +
-                   api_name_n + " = Ref{Csize_t}()")
+    a.julia_pre = (api_name + " = Ref{Ptr{Ptr{Cchar}}}()\n    " + api_name_n +
+                   " = Ref{Csize_t}()")
     a.julia_arg = api_name + ", " + api_name_n
-    a.julia_post = ("tmp_" + api_name + " = unsafe_wrap(Array, " + api_name + "[], " +
-                    api_name_n + "[], own=true)\n    " +
-                    name + " = [unsafe_string(tmp_" + api_name + "[i]) for i in 1:length(tmp_" +
-                    api_name + ") ]")
+    a.julia_post = ("tmp_" + api_name + " = unsafe_wrap(Array, " + api_name +
+                    "[], " + api_name_n + "[], own=true)\n    " + name +
+                    " = [unsafe_string(tmp_" + api_name +
+                    "[i]) for i in 1:length(tmp_" + api_name + ") ]")
     return a
 
+
 def ovectorpair(name, value=None, python_value=None, julia_value=None):
-    a = arg(name, value, python_value, julia_value,
-            ns + "::vectorpair &", "int **", True)
+    a = arg(name, value, python_value, julia_value, ns + "::vectorpair &",
+            "int **", True)
     api_name = "api_" + name + "_"
     api_name_n = api_name + "n_"
     a.c_pre = "    " + ns + "::vectorpair " + api_name + ";\n"
@@ -488,22 +525,23 @@ def ovectorpair(name, value=None, python_value=None, julia_value=None):
     a.cwrap_arg = "&" + api_name + ", " + "&" + api_name_n
     a.cwrap_post = (name + ".resize(" + api_name_n + " / 2); " +
                     "for(size_t i = 0; i < " + api_name_n + " / 2; ++i){ " +
-                    name + "[i].first = " + api_name + "[i * 2 + 0]; " +
-                    name + "[i].second = " + api_name + "[i * 2 + 1]; } " +
-                    ns + "Free(" + api_name + ");\n")
+                    name + "[i].first = " + api_name + "[i * 2 + 0]; " + name +
+                    "[i].second = " + api_name + "[i * 2 + 1]; } " + ns +
+                    "Free(" + api_name + ");\n")
     a.python_pre = api_name + ", " + api_name_n + " = POINTER(c_int)(), c_size_t()"
     a.python_arg = "byref(" + api_name + "), byref(" + api_name_n + ")"
     a.python_return = "_ovectorpair(" + api_name + ", " + api_name_n + ".value)"
     a.julia_ctype = "Ptr{Ptr{Cint}}, Ptr{Csize_t}"
-    a.julia_pre = (api_name + " = Ref{Ptr{Cint}}()\n    " +
-                   api_name_n + " = Ref{Csize_t}()")
+    a.julia_pre = (api_name + " = Ref{Ptr{Cint}}()\n    " + api_name_n +
+                   " = Ref{Csize_t}()")
     a.julia_arg = api_name + ", " + api_name_n
-    a.julia_post = ("tmp_" + api_name + " = unsafe_wrap(Array, " + api_name + "[], " +
-                    api_name_n + "[], own=true)\n    " +
-                    name + " = [ (tmp_" + api_name + "[i], tmp_" + api_name + "[i+1]) " +
-                    "for i in 1:2:length(tmp_" + api_name + ") ]")
+    a.julia_post = ("tmp_" + api_name + " = unsafe_wrap(Array, " + api_name +
+                    "[], " + api_name_n + "[], own=true)\n    " + name +
+                    " = [ (tmp_" + api_name + "[i], tmp_" + api_name +
+                    "[i+1]) " + "for i in 1:2:length(tmp_" + api_name + ") ]")
     return a
 
+
 def ovectorvectorint(name, value=None, python_value=None, julia_value=None):
     a = arg(name, value, python_value, julia_value,
             "std::vector<std::vector<int> > &", "int **", True)
@@ -514,34 +552,37 @@ def ovectorvectorint(name, value=None, python_value=None, julia_value=None):
     a.c_arg = api_name
     a.c_post = ("    vectorvector2ptrptr(" + api_name + ", " + name + ", " +
                 name + "_n, " + name + "_nn);\n")
-    a.c  = "int *** " + name + ", size_t ** " + name + "_n, size_t *" + name + "_nn"
+    a.c = "int *** " + name + ", size_t ** " + name + "_n, size_t *" + name + "_nn"
     a.cwrap_pre = "int **" + api_name + "; size_t *" + api_name_n + ", " + api_name_nn + ";\n"
     a.cwrap_arg = "&" + api_name + ", " + "&" + api_name_n + ", " + "&" + api_name_nn
     a.cwrap_post = (name + ".resize(" + api_name_nn + "); " +
                     "for(size_t i = 0; i < " + api_name_nn + "; ++i){ " +
-                    name + "[i].assign(" + api_name + "[i], " + api_name + "[i] + " +
-                    api_name_n + "[i]); " + ns + "Free(" + api_name + "[i]); } " +
-                    ns + "Free(" + api_name + "); " + ns + "Free(" + api_name_n + ");\n")
-    a.python_pre = (api_name + ", " + api_name_n + ", " + api_name_nn +
-                    " = POINTER(POINTER(c_int))(), POINTER(c_size_t)(), c_size_t()")
-    a.python_arg = ("byref(" + api_name + "), byref(" + api_name_n + "), byref(" +
-                    api_name_nn + ")")
-    a.python_return = ("_ovectorvectorint(" + api_name + ", " + api_name_n + ", " +
-                       api_name_nn + ")")
+                    name + "[i].assign(" + api_name + "[i], " + api_name +
+                    "[i] + " + api_name_n + "[i]); " + ns + "Free(" +
+                    api_name + "[i]); } " + ns + "Free(" + api_name + "); " +
+                    ns + "Free(" + api_name_n + ");\n")
+    a.python_pre = (
+        api_name + ", " + api_name_n + ", " + api_name_nn +
+        " = POINTER(POINTER(c_int))(), POINTER(c_size_t)(), c_size_t()")
+    a.python_arg = ("byref(" + api_name + "), byref(" + api_name_n +
+                    "), byref(" + api_name_nn + ")")
+    a.python_return = ("_ovectorvectorint(" + api_name + ", " + api_name_n +
+                       ", " + api_name_nn + ")")
     a.julia_ctype = "Ptr{Ptr{Ptr{Cint}}}, Ptr{Ptr{Csize_t}}, Ptr{Csize_t}"
-    a.julia_pre = (api_name + " = Ref{Ptr{Ptr{Cint}}}()\n    " +
-                   api_name_n + " = Ref{Ptr{Csize_t}}()\n    " +
-                   api_name_nn + " = Ref{Csize_t}()")
+    a.julia_pre = (api_name + " = Ref{Ptr{Ptr{Cint}}}()\n    " + api_name_n +
+                   " = Ref{Ptr{Csize_t}}()\n    " + api_name_nn +
+                   " = Ref{Csize_t}()")
     a.julia_arg = api_name + ", " + api_name_n + ", " + api_name_nn
-    a.julia_post = ("tmp_" + api_name + " = unsafe_wrap(Array, " + api_name + "[], " +
-                    api_name_nn + "[], own=true)\n    " +
-                    "tmp_" + api_name_n + " = unsafe_wrap(Array, " + api_name_n + "[], " +
-                    api_name_nn + "[], own=true)\n    " +
-                    name + " = [ unsafe_wrap(Array, tmp_" + api_name + "[i], " +
+    a.julia_post = ("tmp_" + api_name + " = unsafe_wrap(Array, " + api_name +
+                    "[], " + api_name_nn + "[], own=true)\n    " + "tmp_" +
+                    api_name_n + " = unsafe_wrap(Array, " + api_name_n +
+                    "[], " + api_name_nn + "[], own=true)\n    " + name +
+                    " = [ unsafe_wrap(Array, tmp_" + api_name + "[i], " +
                     "tmp_" + api_name_n + "[i], own=true) for i in 1:" +
                     api_name_nn + "[] ]")
     return a
 
+
 def ovectorvectorsize(name, value=None, python_value=None, julia_value=None):
     a = arg(name, value, python_value, julia_value,
             "std::vector<std::vector<std::size_t> > &", "size_t **", True)
@@ -552,34 +593,37 @@ def ovectorvectorsize(name, value=None, python_value=None, julia_value=None):
     a.c_arg = api_name
     a.c_post = ("    vectorvector2ptrptr(" + api_name + ", " + name + ", " +
                 name + "_n, " + name + "_nn);\n")
-    a.c  = "size_t *** " + name + ", size_t ** " + name + "_n, size_t *" + name + "_nn"
+    a.c = "size_t *** " + name + ", size_t ** " + name + "_n, size_t *" + name + "_nn"
     a.cwrap_pre = "size_t **" + api_name + "; size_t *" + api_name_n + ", " + api_name_nn + ";\n"
     a.cwrap_arg = "&" + api_name + ", " + "&" + api_name_n + ", " + "&" + api_name_nn
     a.cwrap_post = (name + ".resize(" + api_name_nn + "); " +
                     "for(size_t i = 0; i < " + api_name_nn + "; ++i){ " +
-                    name + "[i].assign(" + api_name + "[i], " + api_name + "[i] + " +
-                    api_name_n + "[i]); " + ns + "Free(" + api_name + "[i]); } " +
-                    ns + "Free(" + api_name + "); " + ns + "Free(" + api_name_n + ");\n")
-    a.python_pre = (api_name + ", " + api_name_n + ", " + api_name_nn +
-                    " = POINTER(POINTER(c_size_t))(), POINTER(c_size_t)(), c_size_t()")
-    a.python_arg = ("byref(" + api_name + "), byref(" + api_name_n + "), byref(" +
-                    api_name_nn + ")")
-    a.python_return = ("_ovectorvectorsize(" + api_name + ", " + api_name_n + ", " +
-                       api_name_nn + ")")
+                    name + "[i].assign(" + api_name + "[i], " + api_name +
+                    "[i] + " + api_name_n + "[i]); " + ns + "Free(" +
+                    api_name + "[i]); } " + ns + "Free(" + api_name + "); " +
+                    ns + "Free(" + api_name_n + ");\n")
+    a.python_pre = (
+        api_name + ", " + api_name_n + ", " + api_name_nn +
+        " = POINTER(POINTER(c_size_t))(), POINTER(c_size_t)(), c_size_t()")
+    a.python_arg = ("byref(" + api_name + "), byref(" + api_name_n +
+                    "), byref(" + api_name_nn + ")")
+    a.python_return = ("_ovectorvectorsize(" + api_name + ", " + api_name_n +
+                       ", " + api_name_nn + ")")
     a.julia_ctype = "Ptr{Ptr{Ptr{Csize_t}}}, Ptr{Ptr{Csize_t}}, Ptr{Csize_t}"
     a.julia_pre = (api_name + " = Ref{Ptr{Ptr{Csize_t}}}()\n    " +
-                   api_name_n + " = Ref{Ptr{Csize_t}}()\n    " +
-                   api_name_nn + " = Ref{Csize_t}()")
+                   api_name_n + " = Ref{Ptr{Csize_t}}()\n    " + api_name_nn +
+                   " = Ref{Csize_t}()")
     a.julia_arg = api_name + ", " + api_name_n + ", " + api_name_nn
-    a.julia_post = ("tmp_" + api_name + " = unsafe_wrap(Array, " + api_name + "[], " +
-                    api_name_nn + "[], own=true)\n    " +
-                    "tmp_" + api_name_n + " = unsafe_wrap(Array, " + api_name_n + "[], " +
-                    api_name_nn + "[], own=true)\n    " +
-                    name + " = [ unsafe_wrap(Array, tmp_" + api_name + "[i], " +
+    a.julia_post = ("tmp_" + api_name + " = unsafe_wrap(Array, " + api_name +
+                    "[], " + api_name_nn + "[], own=true)\n    " + "tmp_" +
+                    api_name_n + " = unsafe_wrap(Array, " + api_name_n +
+                    "[], " + api_name_nn + "[], own=true)\n    " + name +
+                    " = [ unsafe_wrap(Array, tmp_" + api_name + "[i], " +
                     "tmp_" + api_name_n + "[i], own=true) for i in 1:" +
                     api_name_nn + "[] ]")
     return a
 
+
 def ovectorvectordouble(name, value=None, python_value=None, julia_value=None):
     a = arg(name, value, python_value, julia_value,
             "std::vector<std::vector<double> > &", "double **", True)
@@ -588,37 +632,40 @@ def ovectorvectordouble(name, value=None, python_value=None, julia_value=None):
     api_name_nn = api_name + "nn_"
     a.c_pre = "    std::vector<std::vector<double> > " + api_name + ";\n"
     a.c_arg = api_name
-    a.c_post = ("    vectorvector2ptrptr(" + api_name + ", " + name + ", " + name +
-                "_n, " + name + "_nn);\n")
-    a.c  = "double *** " + name + ", size_t ** " + name + "_n, size_t *" + name + "_nn"
+    a.c_post = ("    vectorvector2ptrptr(" + api_name + ", " + name + ", " +
+                name + "_n, " + name + "_nn);\n")
+    a.c = "double *** " + name + ", size_t ** " + name + "_n, size_t *" + name + "_nn"
     a.cwrap_pre = ("double **" + api_name + "; size_t *" + api_name_n + ", " +
                    api_name_nn + ";\n")
     a.cwrap_arg = "&" + api_name + ", " + "&" + api_name_n + ", " + "&" + api_name_nn
     a.cwrap_post = (name + ".resize(" + api_name_nn + "); " +
                     "for(size_t i = 0; i < " + api_name_nn + "; ++i){ " +
-                    name + "[i].assign(" + api_name + "[i], " + api_name + "[i] + " +
-                    api_name_n + "[i]); " + ns + "Free(" + api_name + "[i]); } " +
-                    ns + "Free(" + api_name + "); " + ns + "Free(" + api_name_n + ");\n")
-    a.python_pre = (api_name + ", " + api_name_n + ", " + api_name_nn +
-                    " = POINTER(POINTER(c_double))(), POINTER(c_size_t)(), c_size_t()")
-    a.python_arg = ("byref(" + api_name + "), byref(" + api_name_n + "), byref(" +
-                    api_name_nn + ")")
-    a.python_return = ("_ovectorvectordouble(" + api_name + ", " + api_name_n + ", " +
-                       api_name_nn + ")")
+                    name + "[i].assign(" + api_name + "[i], " + api_name +
+                    "[i] + " + api_name_n + "[i]); " + ns + "Free(" +
+                    api_name + "[i]); } " + ns + "Free(" + api_name + "); " +
+                    ns + "Free(" + api_name_n + ");\n")
+    a.python_pre = (
+        api_name + ", " + api_name_n + ", " + api_name_nn +
+        " = POINTER(POINTER(c_double))(), POINTER(c_size_t)(), c_size_t()")
+    a.python_arg = ("byref(" + api_name + "), byref(" + api_name_n +
+                    "), byref(" + api_name_nn + ")")
+    a.python_return = ("_ovectorvectordouble(" + api_name + ", " + api_name_n +
+                       ", " + api_name_nn + ")")
     a.julia_ctype = "Ptr{Ptr{Ptr{Cdouble}}}, Ptr{Ptr{Csize_t}}, Ptr{Csize_t}"
     a.julia_pre = (api_name + " = Ref{Ptr{Ptr{Cdouble}}}()\n    " +
-                   api_name_n + " = Ref{Ptr{Csize_t}}()\n    " +
-                   api_name_nn + " = Ref{Csize_t}()")
+                   api_name_n + " = Ref{Ptr{Csize_t}}()\n    " + api_name_nn +
+                   " = Ref{Csize_t}()")
     a.julia_arg = api_name + ", " + api_name_n + ", " + api_name_nn
-    a.julia_post = ("tmp_" + api_name + " = unsafe_wrap(Array, " + api_name + "[], " +
-                    api_name_nn + "[], own=true)\n    " +
-                    "tmp_" + api_name_n + " = unsafe_wrap(Array, " + api_name_n + "[], " +
-                    api_name_nn + "[], own=true)\n    " +
-                    name + " = [ unsafe_wrap(Array, tmp_" + api_name + "[i], " +
+    a.julia_post = ("tmp_" + api_name + " = unsafe_wrap(Array, " + api_name +
+                    "[], " + api_name_nn + "[], own=true)\n    " + "tmp_" +
+                    api_name_n + " = unsafe_wrap(Array, " + api_name_n +
+                    "[], " + api_name_nn + "[], own=true)\n    " + name +
+                    " = [ unsafe_wrap(Array, tmp_" + api_name + "[i], " +
                     "tmp_" + api_name_n + "[i], own=true) for i in 1:" +
                     api_name_nn + "[] ]")
     return a
 
+
 def ovectorvectorpair(name, value=None, python_value=None, julia_value=None):
     a = arg(name, value, python_value, julia_value,
             "std::vector<" + ns + "::vectorpair> &", "int **", True)
@@ -637,33 +684,35 @@ def ovectorvectorpair(name, value=None, python_value=None, julia_value=None):
                     name + "[i].resize(" + api_name_n + "[i] / 2); " +
                     "for(size_t j = 0; j < " + api_name_n + "[i] / 2; ++j){ " +
                     name + "[i][j].first = " + api_name + "[i][j * 2 + 0]; " +
-                    name + "[i][j].second = " + api_name + "[i][j * 2 + 1]; } " +
-                    ns + "Free(" + api_name + "[i]); } " +
-                    ns + "Free(" + api_name + "); " + ns + "Free(" + api_name_n + ");\n")
-    a.python_pre = (api_name + ", " + api_name_n + ", " + api_name_nn +
-                    " = POINTER(POINTER(c_int))(), POINTER(c_size_t)(), c_size_t()")
-    a.python_arg = ("byref(" + api_name + "), byref(" + api_name_n + "), byref(" +
-                    api_name_nn + ")")
-    a.python_return = ("_ovectorvectorpair(" + api_name + ", " + api_name_n + ", " +
-                       api_name_nn + ")")
+                    name + "[i][j].second = " + api_name +
+                    "[i][j * 2 + 1]; } " + ns + "Free(" + api_name +
+                    "[i]); } " + ns + "Free(" + api_name + "); " + ns +
+                    "Free(" + api_name_n + ");\n")
+    a.python_pre = (
+        api_name + ", " + api_name_n + ", " + api_name_nn +
+        " = POINTER(POINTER(c_int))(), POINTER(c_size_t)(), c_size_t()")
+    a.python_arg = ("byref(" + api_name + "), byref(" + api_name_n +
+                    "), byref(" + api_name_nn + ")")
+    a.python_return = ("_ovectorvectorpair(" + api_name + ", " + api_name_n +
+                       ", " + api_name_nn + ")")
     a.julia_ctype = "Ptr{Ptr{Ptr{Cint}}}, Ptr{Ptr{Csize_t}}, Ptr{Csize_t}"
-    a.julia_pre = (api_name + " = Ref{Ptr{Ptr{Cint}}}()\n    " +
-                   api_name_n + " = Ref{Ptr{Csize_t}}()\n    " +
-                   api_name_nn + " = Ref{Csize_t}()")
+    a.julia_pre = (api_name + " = Ref{Ptr{Ptr{Cint}}}()\n    " + api_name_n +
+                   " = Ref{Ptr{Csize_t}}()\n    " + api_name_nn +
+                   " = Ref{Csize_t}()")
     a.julia_arg = api_name + ", " + api_name_n + ", " + api_name_nn
-    a.julia_post = ("tmp_" + api_name + " = unsafe_wrap(Array, " + api_name + "[], " +
-                    api_name_nn + "[], own=true)\n    " +
-                    "tmp_" + api_name_n + " = unsafe_wrap(Array, " + api_name_n + "[], " +
-                    api_name_nn + "[], own=true)\n    " +
-                    name + " = Vector{Tuple{Cint,Cint}}[]\n    " +
-                    "resize!(" + name + ", " + api_name_nn + "[])\n    " +
-                    "for i in 1:" + api_name_nn + "[]\n    " +
-                    "    tmp = unsafe_wrap(Array, tmp_" + api_name + "[i], tmp_" +
-                    api_name_n + "[i], own=true)\n    " +
-                    "    " + name + "[i] = [(tmp[i], tmp[i+1]) for i in 1:2:length(tmp)]\n    " +
-                    "end")
+    a.julia_post = (
+        "tmp_" + api_name + " = unsafe_wrap(Array, " + api_name + "[], " +
+        api_name_nn + "[], own=true)\n    " + "tmp_" + api_name_n +
+        " = unsafe_wrap(Array, " + api_name_n + "[], " + api_name_nn +
+        "[], own=true)\n    " + name + " = Vector{Tuple{Cint,Cint}}[]\n    " +
+        "resize!(" + name + ", " + api_name_nn + "[])\n    " + "for i in 1:" +
+        api_name_nn + "[]\n    " + "    tmp = unsafe_wrap(Array, tmp_" +
+        api_name + "[i], tmp_" + api_name_n + "[i], own=true)\n    " + "    " +
+        name + "[i] = [(tmp[i], tmp[i+1]) for i in 1:2:length(tmp)]\n    " +
+        "end")
     return a
 
+
 def argcargv():
     a = arg("", None, None, None, "", "", False)
     a.cpp = "int argc = 0, char ** argv = 0"
@@ -682,8 +731,8 @@ def argcargv():
     a.texi = "(argc = 0)}, @code{argv = []"
     return a
 
-class Module:
 
+class Module:
     def __init__(self, name, doc):
         self.name = name
         self.doc = doc
@@ -702,7 +751,7 @@ class Module:
         return module
 
 
-cpp_header="""// {0}
+cpp_header = """// {0}
 //
 // See the LICENSE.txt file for license information. Please report all
 // issues on {1}
@@ -754,10 +803,10 @@ namespace {6} {{
 
 """
 
-cpp_footer="""#endif
+cpp_footer = """#endif
 """
 
-c_header="""/*
+c_header = """/*
  * {0}
  *
  * See the LICENSE.txt file for license information. Please report all
@@ -796,11 +845,11 @@ c_header="""/*
 {2}_API void *{6}Malloc(size_t n);
 """
 
-c_footer="""
+c_footer = """
 #endif
 """
 
-c_cpp_header="""// {0}
+c_cpp_header = """// {0}
 //
 // See the LICENSE.txt file for license information. Please report all
 // issues on {1}
@@ -824,7 +873,7 @@ extern \"C\" {{
 }}
 """
 
-c_cpp_utils="""
+c_cpp_utils = """
 void vectorvectorpair2intptrptr(const std::vector<{0}::vectorpair > &v, int ***p, size_t **size, size_t *sizeSize)
 {{
   *p = (int**){0}Malloc(sizeof(int*) * v.size());
@@ -835,7 +884,7 @@ void vectorvectorpair2intptrptr(const std::vector<{0}::vectorpair > &v, int ***p
 }}
 """
 
-cwrap_header="""// {0}
+cwrap_header = """// {0}
 //
 // See the LICENSE.txt file for license information. Please report all
 // issues on {1}
@@ -885,7 +934,7 @@ namespace {6} {{
 
 """
 
-cwrap_utils="""
+cwrap_utils = """
 template<typename t>
 {1}void vector2ptr(const std::vector<t> &v, t **p, size_t *size)
 {{
@@ -927,7 +976,7 @@ template<typename t>
 }}
 """
 
-cwrap_footer="""#endif
+cwrap_footer = """#endif
 """
 
 python_header = """# {0}
@@ -1128,13 +1177,20 @@ julia_header = """# {0}
 # Julia types. See `tutorial/julia' and `demos/api' for examples.
 """
 
-def capi(s): return s[:1].upper() + s[1:]
 
-class API:
+def capi(s):
+    return s[:1].upper() + s[1:]
+
 
-    def __init__(self, version_major, version_minor, namespace="gmsh", code="Gmsh",
-                 copyright="Gmsh - Copyright (C) 1997-2020 C. Geuzaine, J.-F. Remacle",
-                 issues="https://gitlab.onelab.info/gmsh/gmsh/issues."):
+class API:
+    def __init__(
+        self,
+        version_major,
+        version_minor,
+        namespace="gmsh",
+        code="Gmsh",
+        copyright="Gmsh - Copyright (C) 1997-2020 C. Geuzaine, J.-F. Remacle",
+        issues="https://gitlab.onelab.info/gmsh/gmsh/issues."):
         self.version_major = version_major
         self.version_minor = version_minor
         global ns
@@ -1157,21 +1213,29 @@ class API:
             indent += "  "
             for rtype, name, args, doc, special in module.fs:
                 rt = rtype.rcpp_type if rtype else "void"
-                f.write(indent + "// " + cpp_mpath + name + "\n" + indent + "//\n")
-                f.write(indent + "// " + ("\n" + indent + "// ").join(
-                    textwrap.wrap(doc, 80-len(indent))) + "\n")
-
-                fnameapi = indent + ns.upper() + "_API " + rt + " " + name + "(";
+                f.write(indent + "// " + cpp_mpath + name + "\n" + indent +
+                        "//\n")
+                f.write(indent + "// " +
+                        ("\n" + indent +
+                         "// ").join(textwrap.wrap(doc, 80 - len(indent))) +
+                        "\n")
+
+                fnameapi = indent + ns.upper(
+                ) + "_API " + rt + " " + name + "("
                 f.write(fnameapi)
                 if args:
-                    f.write((",\n" + ' ' * len(fnameapi)).join(a.cpp for a in args))
+                    f.write((",\n" + ' ' * len(fnameapi)).join(a.cpp
+                                                               for a in args))
                 f.write(");\n\n")
             for m in module.submodules:
                 write_module(m, indent, cpp_mpath)
             f.write(indent[:-2] + "} // namespace " + module.name + "\n\n")
+
         with open(ns + ".h", "w") as f:
-            f.write(cpp_header.format(self.copyright, self.issues, ns.upper(), self.code,
-                                      self.version_major, self.version_minor, ns))
+            f.write(
+                cpp_header.format(self.copyright, self.issues, ns.upper(),
+                                  self.code, self.version_major,
+                                  self.version_minor, ns))
             for m in self.modules:
                 write_module(m, "", "")
             f.write(cpp_footer)
@@ -1189,57 +1253,63 @@ class API:
             for rtype, name, args, doc, special in module.fs:
                 # *c.h
                 fname = c_namespace + name[0].upper() + name[1:]
-                f.write("\n/* " + "\n * ".join(textwrap.wrap(doc, 75)) + " */\n")
+                f.write("\n/* " + "\n * ".join(textwrap.wrap(doc, 75)) +
+                        " */\n")
                 fnameapi = ns.upper() + "_API " + (rtype.rc_type if rtype else
                                                    "void") + " " + fname + "("
-                f.write(fnameapi
-                        + (",\n" + ' ' * len(fnameapi)).join(
-                            list((a.c for a in args + (oint("ierr"), ))))
-                        + ");\n")
+                f.write(fnameapi + (",\n" + ' ' * len(fnameapi)).join(
+                    list((a.c for a in args + (oint("ierr"), )))) + ");\n")
                 if "rawc" not in special:
                     # *c.cpp
-                    fc.write(ns.upper() + "_API " + (rtype.rc_type if rtype else "void"))
-                    fc.write(" " + fname + "(" +
-                             ", ".join(list((a.c for a in args + (oint("ierr"), )))) +
-                             ")\n{\n")
+                    fc.write(ns.upper() + "_API " +
+                             (rtype.rc_type if rtype else "void"))
+                    fc.write(
+                        " " + fname + "(" +
+                        ", ".join(list((a.c for a in args +
+                                        (oint("ierr"), )))) + ")\n{\n")
                     if rtype:
                         fc.write("  " + rtype.rc_type + " result_api_ = 0;\n")
-                    fc.write("  if(ierr) *ierr = 0;\n");
-                    fc.write("  try {\n");
+                    fc.write("  if(ierr) *ierr = 0;\n")
+                    fc.write("  try {\n")
                     fc.write("".join((a.c_pre for a in args)))
                     fc.write("    ")
                     if rtype:
                         fc.write("result_api_ = ")
-                    fc.write(cpp_namespace + name + "(" + ", ".join(
-                        list((a.c_arg for a in args))) +
-                             ");\n")
+                    fc.write(cpp_namespace + name + "(" +
+                             ", ".join(list((a.c_arg for a in args))) + ");\n")
                     fc.write("".join((a.c_post for a in args)))
                     fc.write("  }\n  catch(int api_ierr_){\n    " +
-                             "if(ierr) *ierr = api_ierr_;\n  }\n");
+                             "if(ierr) *ierr = api_ierr_;\n  }\n")
                     if rtype:
-                        fc.write("  return result_api_;\n");
+                        fc.write("  return result_api_;\n")
                     fc.write("}\n\n")
                 # *.h_cwrap
-                fcwrap.write(indent + "// " + ("\n" + indent + "// ").join
-                             (textwrap.wrap(doc, 80-len(indent))) + "\n")
+                fcwrap.write(indent + "// " +
+                             ("\n" + indent +
+                              "// ").join(textwrap.wrap(doc, 80 -
+                                                        len(indent))) + "\n")
                 rt = rtype.rcpp_type if rtype else "void"
-                fnameapi = indent + "inline " + rt + " " + name + "(";
+                fnameapi = indent + "inline " + rt + " " + name + "("
                 fcwrap.write(fnameapi)
                 if args:
-                    fcwrap.write((",\n" + ' ' * len(fnameapi)).join(a.cpp for a in args))
-                fcwrap.write(")\n" + indent + "{\n" + indent + "  int ierr = 0;\n")
+                    fcwrap.write(
+                        (",\n" + ' ' * len(fnameapi)).join(a.cpp
+                                                           for a in args))
+                fcwrap.write(")\n" + indent + "{\n" + indent +
+                             "  int ierr = 0;\n")
                 for a in args:
                     if a.cwrap_pre:
                         fcwrap.write(indent + "  " + a.cwrap_pre)
                 fcwrap.write(indent + "  ")
                 if rtype:
                     fcwrap.write(rt + " result_api_ = ")
-                fcwrap.write(fname + "(" + ", ".join((a.cwrap_arg for a in args)))
+                fcwrap.write(fname + "(" + ", ".join((a.cwrap_arg
+                                                      for a in args)))
                 if args:
                     fcwrap.write(", &ierr);\n")
                 else:
                     fcwrap.write("&ierr);\n")
-                fcwrap.write(indent + "  " + "if(ierr) throw ierr;\n");
+                fcwrap.write(indent + "  " + "if(ierr) throw ierr;\n")
                 for a in args:
                     if a.cwrap_post:
                         fcwrap.write(indent + "  " + a.cwrap_post)
@@ -1248,22 +1318,28 @@ class API:
                 fcwrap.write(indent + "}\n\n")
             for m in module.submodules:
                 write_module(m, c_namespace, cpp_namespace, indent)
-            fcwrap.write(indent[:-2] + "} // namespace " + module.name + "\n\n")
+            fcwrap.write(indent[:-2] + "} // namespace " + module.name +
+                         "\n\n")
+
         with open(ns + "c.h", "w") as f:
             with open(ns + "c.cpp", "w") as fc:
                 with open(ns + ".h_cwrap", "w") as fcwrap:
-                    f.write(c_header.format(self.copyright, self.issues, ns.upper(),
-                                            self.code, self.version_major,
-                                            self.version_minor, ns))
-                    fc.write(c_cpp_header.format(self.copyright, self.issues, ns,
-                                                 ns.upper()))
+                    f.write(
+                        c_header.format(self.copyright, self.issues,
+                                        ns.upper(), self.code,
+                                        self.version_major, self.version_minor,
+                                        ns))
+                    fc.write(
+                        c_cpp_header.format(self.copyright, self.issues, ns,
+                                            ns.upper()))
                     fc.write(cwrap_utils.format(ns, ""))
                     fc.write(c_cpp_utils.format(ns))
                     fc.write("\n")
-                    fcwrap.write(cwrap_header.format(self.copyright, self.issues,
-                                                     ns.upper(), self.code,
-                                                     self.version_major,
-                                                     self.version_minor, ns))
+                    fcwrap.write(
+                        cwrap_header.format(self.copyright, self.issues,
+                                            ns.upper(), self.code,
+                                            self.version_major,
+                                            self.version_minor, ns))
                     fcwrap.write("namespace " + ns + " {\n")
                     s = cwrap_utils.format(ns, "inline ").split('\n')
                     for line in s:
@@ -1277,6 +1353,7 @@ class API:
     def write_python(self):
         def parg(a):
             return a.name + (("=" + a.python_value) if a.python_value else "")
+
         def write_function(f, fun, c_mpath, py_mpath, indent):
             (rtype, name, args, doc, special) = fun
             if "onlycc++" in special: return
@@ -1285,34 +1362,34 @@ class API:
             f.write("\n")
             if c_mpath != ns:
                 f.write(indent + "@staticmethod\n")
-            f.write(indent + "def " + name + "("
-                   + ", ".join((parg(a) for a in iargs))
-                   + "):\n")
+            f.write(indent + "def " + name + "(" +
+                    ", ".join((parg(a) for a in iargs)) + "):\n")
             indent += "    "
             f.write(indent + '"""\n')
-            f.write(indent + py_mpath + name + "(" + ", ".join(parg(a) for a in iargs)
-                    + ")\n\n")
-            f.write(indent + ("\n" + indent).join(textwrap.wrap(doc, 75)) + "\n")
+            f.write(indent + py_mpath + name + "(" +
+                    ", ".join(parg(a) for a in iargs) + ")\n\n")
+            f.write(indent + ("\n" + indent).join(textwrap.wrap(doc, 75)) +
+                    "\n")
             if rtype or oargs:
                 f.write("\n" + indent + "Return " + ", ".join(
-                    ([("an " if rtype.rtexi_type == "integer value" else "a ") +
-                      rtype.rtexi_type] if rtype else[])
-                   + [("`" + a.name + "'") for a in oargs])
-               + ".\n")
+                    ([("an " if rtype.rtexi_type == "integer value" else "a "
+                       ) + rtype.rtexi_type] if rtype else []) +
+                    [("`" + a.name + "'") for a in oargs]) + ".\n")
             f.write(indent + '"""\n')
             for a in args:
                 if a.python_pre: f.write(indent + a.python_pre + "\n")
             f.write(indent + "ierr = c_int()\n")
-            f.write(indent + "api__result__ = " if ((rtype is oint) or
-                                                    (rtype is odouble)) else (indent))
+            f.write(indent + "api__result__ = " if (
+                (rtype is oint) or (rtype is odouble)) else (indent))
             c_name = c_mpath + name[0].upper() + name[1:]
             f.write("lib." + c_name + "(\n    " + indent +
                     (",\n" + indent + "    ").join(
-                        tuple((a.python_arg for a in args)) + ("byref(ierr)", )) +
-                    ")\n")
+                        tuple((a.python_arg
+                               for a in args)) + ("byref(ierr)", )) + ")\n")
             f.write(indent + "if ierr.value != 0:\n")
             f.write(indent + "    raise ValueError(\n")
-            f.write(indent + '        "' + c_name + ' returned non-zero error code: ",\n')
+            f.write(indent + '        "' + c_name +
+                    ' returned non-zero error code: ",\n')
             f.write(indent + "        ierr.value)\n")
             r = (["api__result__"]) if rtype else []
             r += list((o.python_return for o in oargs))
@@ -1322,6 +1399,7 @@ class API:
                 else:
                     f.write(indent + "return (\n" + indent + "    " +
                             (",\n" + indent + "    ").join(r) + ")\n")
+
         def write_module(f, m, c_mpath, py_mpath, indent):
             if c_mpath:
                 c_mpath += m.name[0].upper() + m.name[1:]
@@ -1335,51 +1413,58 @@ class API:
                 f.write("\n\n" + indent + "class " + module.name + ":\n")
                 indentm = indent + "    "
                 f.write(indentm + '"""\n')
-                f.write(indentm + ("\n" + indentm).join
-                        (textwrap.wrap(capi(module.doc), 75)) + "\n")
+                f.write(indentm +
+                        ("\n" +
+                         indentm).join(textwrap.wrap(capi(module.doc), 75)) +
+                        "\n")
                 f.write(indentm + '"""\n')
                 write_module(f, module, c_mpath, py_mpath, indentm)
+
         with open(ns + ".py", "w") as f:
-            f.write(python_header.format(self.copyright, self.issues, self.code,
-                                         self.version_major, self.version_minor,
-                                         ns.upper(), ns))
+            f.write(
+                python_header.format(self.copyright, self.issues, self.code,
+                                     self.version_major, self.version_minor,
+                                     ns.upper(), ns))
             for module in self.modules:
                 write_module(f, module, "", "", "")
 
     def write_julia(self):
         def parg(a):
             return a.name + ((" = " + a.julia_value) if a.julia_value else "")
+
         def write_function(f, fun, c_mpath, jl_mpath):
             (rtype, name, args, doc, special) = fun
             if "onlycc++" in special: return
             iargs = list(a for a in args if not a.out)
             oargs = list(a for a in args if a.out)
             f.write('\n"""\n    ')
-            f.write(jl_mpath + name + "(" + ", ".join(parg(a) for a in iargs) + ")\n\n")
+            f.write(jl_mpath + name + "(" + ", ".join(parg(a) for a in iargs) +
+                    ")\n\n")
             f.write("\n".join(textwrap.wrap(doc, 80)).replace("'", "`") + "\n")
             if rtype or oargs:
-                f.write("\nReturn " + ", ".join(
-                    ([("an " if rtype.rtexi_type == "integer value" else "a ") +
-                      rtype.rtexi_type] if rtype else[])
-                   + [("`" + a.name + "`") for a in oargs])
-               + ".\n")
+                f.write("\nReturn " + ", ".join((
+                    [("an " if rtype.rtexi_type == "integer value" else "a ") +
+                     rtype.rtexi_type] if rtype else []) +
+                                                [("`" + a.name + "`")
+                                                 for a in oargs]) + ".\n")
             f.write('"""\n')
-            f.write("function " + name + "("
-                   + ", ".join((parg(a) for a in iargs))
-                   + ")\n")
+            f.write("function " + name + "(" +
+                    ", ".join((parg(a) for a in iargs)) + ")\n")
             for a in args:
                 if a.julia_pre: f.write("    " + a.julia_pre + "\n")
             f.write("    ierr = Ref{Cint}()\n    ")
-            f.write("api__result__ = " if ((rtype is oint) or (rtype is odouble)) else "")
+            f.write("api__result__ = " if (
+                (rtype is oint) or (rtype is odouble)) else "")
             c_name = c_mpath + name[0].upper() + name[1:]
             f.write("ccall((:" + c_name + ", " +
                     ("" if c_mpath == ns else ns + ".") + "lib), " +
                     ("Cvoid" if rtype is None else rtype.rjulia_type) + ",\n" +
-                    " " * 10 + "(" + ", ".join(
-                        (tuple(a.julia_ctype for a in args) + ("Ptr{Cint}",))) +
-                    ("," if not len(args) else "") + "),\n" +
-                    " " * 10 + ", ".join(
-                        tuple(a.julia_arg for a in args) + ("ierr",)) + ")\n")
+                    " " * 10 + "(" + ", ".join((tuple(a.julia_ctype
+                                                      for a in args) +
+                                                ("Ptr{Cint}", ))) +
+                    ("," if not len(args) else "") + "),\n" + " " * 10 +
+                    ", ".join(tuple(a.julia_arg
+                                    for a in args) + ("ierr", )) + ")\n")
             f.write('    ierr[] != 0 && error("' + c_name +
                     ' returned non-zero error code: $(ierr[])")\n')
             for a in args:
@@ -1392,6 +1477,7 @@ class API:
             else:
                 f.write(", ".join(r))
             f.write("\nend\n")
+
         def write_module(f, m, c_mpath, jl_mpath, level):
             f.write('\n"""\n    ')
             f.write("module " + jl_mpath + m.name + "\n\n")
@@ -1399,18 +1485,20 @@ class API:
             f.write('"""\n')
             f.write("module " + m.name + "\n\n")
             if level == 1:
-                f.write('const {0}_API_VERSION = "{1}.{2}"\n'.
-                        format(ns.upper(), self.version_major, self.version_minor))
-                f.write('const {0}_API_VERSION_MAJOR = {1}\n'.
-                        format(ns.upper(), self.version_major))
-                f.write('const {0}_API_VERSION_MINOR = {1}\n'.
-                        format(ns.upper(), self.version_minor))
+                f.write('const {0}_API_VERSION = "{1}.{2}"\n'.format(
+                    ns.upper(), self.version_major, self.version_minor))
+                f.write('const {0}_API_VERSION_MAJOR = {1}\n'.format(
+                    ns.upper(), self.version_major))
+                f.write('const {0}_API_VERSION_MINOR = {1}\n'.format(
+                    ns.upper(), self.version_minor))
                 f.write('const libdir = dirname(@__FILE__)\n')
-                f.write('const libname = Sys.iswindows() ? "' + ns +
-                        '-{0}.{1}'.format(self.version_major, self.version_minor) +
-                        '" : "lib' + ns + '"\n')
+                f.write(
+                    'const libname = Sys.iswindows() ? "' + ns +
+                    '-{0}.{1}'.format(self.version_major, self.version_minor) +
+                    '" : "lib' + ns + '"\n')
                 f.write('import Libdl\n')
-                f.write('const lib = Libdl.find_library([libname], [libdir])\n')
+                f.write(
+                    'const lib = Libdl.find_library([libname], [libdir])\n')
             else:
                 f.write("import " + ("." * level) + ns + "\n")
             if c_mpath:
@@ -1424,9 +1512,11 @@ class API:
             for module in m.submodules:
                 write_module(f, module, c_mpath, jl_mpath, level + 1)
             f.write("\nend # end of module " + m.name + "\n")
+
         with open(ns + ".jl", "w") as f:
-            f.write(julia_header.format(self.copyright, self.issues, self.code,
-                                        self.version_major, self.version_minor))
+            f.write(
+                julia_header.format(self.copyright, self.issues, self.code,
+                                    self.version_major, self.version_minor))
             for module in self.modules:
                 write_module(f, module, "", "", 1)
 
@@ -1436,8 +1526,10 @@ class API:
                 return int(s)
             except:
                 return s
+
         def alphanum_key(s):
-            return [ tryint(c) for c in re.split('([0-9]+)', s) ]
+            return [tryint(c) for c in re.split('([0-9]+)', s)]
+
         def get_file_data(path, ext):
             data = []
             for r, d, f in os.walk(path):
@@ -1450,49 +1542,57 @@ class API:
                         data.append([filename, contents])
             data.sort(key=lambda x: alphanum_key(x[0]))
             return data
+
         def find_function(name, data):
             match = []
             for file in data:
                 l = 0
                 for line in file[1]:
-                    l = l+1
+                    l = l + 1
                     if re.search(name, line):
                         match.append((file[0], l))
-                        break # keep only one match per file
+                        break  # keep only one match per file
             return match
-        def write_module(module, path, node, node_next, node_prev, cpp_data, py_data):
-            f.write("@node " + node + ", " + node_next + ", " + node_prev + ", Gmsh API\n");
-            f.write("@section Namespace @code{" + path + "}: " + module.doc + "\n\n");
-            f.write("@ftable @code\n");
+
+        def write_module(module, path, node, node_next, node_prev, cpp_data,
+                         py_data):
+            f.write("@node " + node + ", " + node_next + ", " + node_prev +
+                    ", Gmsh API\n")
+            f.write("@section Namespace @code{" + path + "}: " + module.doc +
+                    "\n\n")
+            f.write("@ftable @code\n")
             for rtype, name, args, doc, special in module.fs:
-                f.write("@item " + path + '/' + name + "\n");
+                f.write("@item " + path + '/' + name + "\n")
                 tdoc = doc.replace("`", "@code{").replace("'", "}")
                 f.write("\n".join(textwrap.wrap(tdoc, 80)) + "\n\n")
-                f.write("@table @asis\n");
+                f.write("@table @asis\n")
                 iargs = list(a for a in args if not a.out)
                 oargs = list(a for a in args if a.out)
-                f.write("@item " + "Input:\n" +
-                        (", ".join(("@code{" + iarg.texi + "}") for iarg in iargs)
-                         if len(iargs) else "-") + "\n")
-                f.write("@item " + "Output:\n" +
-                        (", ".join(("@code{" + oarg.name + "}") for oarg in oargs)
-                         if len(oargs) else "-") + "\n")
+                f.write("@item " + "Input:\n" + (", ".join(
+                    ("@code{" + iarg.texi + "}")
+                    for iarg in iargs) if len(iargs) else "-") + "\n")
+                f.write("@item " + "Output:\n" + (", ".join(
+                    ("@code{" + oarg.name + "}")
+                    for oarg in oargs) if len(oargs) else "-") + "\n")
                 f.write("@item " + "Return:\n" +
                         (rtype.rtexi_type if rtype else "-") + "\n")
                 cpp_name = path.replace('/', '::') + '::' + name + '\('
                 cpp = find_function(cpp_name, cpp_data)
                 py_name = path.replace('/', '.') + '.' + name + '\('
                 py = find_function(py_name, py_data)
+
                 def write_matches(lang, matches):
                     git = 'https://gitlab.onelab.info/gmsh/gmsh/tree/master/'
                     f.write(lang + ' (')
-                    for i in range(min(5, len(matches))): # write max 5 matches
+                    for i in range(min(5,
+                                       len(matches))):  # write max 5 matches
                         if i > 0: f.write(', ')
-                        f.write('@url{' + git + matches[i][0][3:] +
-                                '#L' + str(matches[i][1]) + ',' +
+                        f.write('@url{' + git + matches[i][0][3:] + '#L' +
+                                str(matches[i][1]) + ',' +
                                 os.path.basename(matches[i][0]) + '}')
                     if len(matches) > 5: f.write(', ...')
                     f.write(')')
+
                 if len(cpp) or len(py):
                     f.write("@item " + "Examples:\n")
                     if len(cpp):
@@ -1501,17 +1601,23 @@ class API:
                     if len(py):
                         write_matches("Python", py)
                     f.write("\n")
-                f.write("@end table\n\n");
-            f.write("@end ftable\n\n");
+                f.write("@end table\n\n")
+            f.write("@end ftable\n\n")
+
         with open("api.texi", "w") as f:
-            f.write("@c This file was generated by api/gen.py: do not edit manually!\n\n")
+            f.write(
+                "@c This file was generated by api/gen.py: do not edit manually!\n\n"
+            )
+
             def flatten_module(flat, module, path):
                 p = path + ("/" if len(path) else "") + module.name
                 flat.append((module, p))
                 for m in module.submodules:
                     flatten_module(flat, m, p)
+
             def node_name(n):
                 return "Namespace " + n[1]
+
             flat = []
             for m in self.modules:
                 flatten_module(flat, m, "")
@@ -1526,6 +1632,6 @@ class API:
             py_data.extend(get_file_data('../demos/api', '.py'))
             for i in range(N):
                 write_module(flat[i][0], flat[i][1], node_name(flat[i]),
-                             "" if i == N - 1 else node_name(flat[i+1]),
-                             "" if i == 0 else node_name(flat[i-1]),
+                             "" if i == N - 1 else node_name(flat[i + 1]),
+                             "" if i == 0 else node_name(flat[i - 1]),
                              cpp_data, py_data)
diff --git a/demos/api/adapt_mesh.py b/demos/api/adapt_mesh.py
index 809f6c598232bb321ec939587b348ebb6b5eaaeb..cf0feb9371b3fa02466b8ecd98655b9724857ed9 100644
--- a/demos/api/adapt_mesh.py
+++ b/demos/api/adapt_mesh.py
@@ -4,39 +4,38 @@ import gmsh
 
 
 def triangle_max_edge(x):
-    a = np.sum((x[:,0,:]-x[:,1,:])**2,1)**0.5
-    b = np.sum((x[:,0,:]-x[:,2,:])**2,1)**0.5
-    c = np.sum((x[:,1,:]-x[:,2,:])**2,1)**0.5
-    return np.maximum(a,np.maximum(b,c))
+    a = np.sum((x[:, 0, :] - x[:, 1, :])**2, 1)**0.5
+    b = np.sum((x[:, 0, :] - x[:, 2, :])**2, 1)**0.5
+    c = np.sum((x[:, 1, :] - x[:, 2, :])**2, 1)**0.5
+    return np.maximum(a, np.maximum(b, c))
 
 
 class Mesh:
-
     def __init__(self):
         self.vtags, vxyz, _ = gmsh.model.mesh.getNodes()
-        self.vxyz = vxyz.reshape((-1,3))
-        vmap = dict({j:i for i,j in enumerate(self.vtags)})
+        self.vxyz = vxyz.reshape((-1, 3))
+        vmap = dict({j: i for i, j in enumerate(self.vtags)})
         self.triangles_tags, evtags = gmsh.model.mesh.getElementsByType(2)
         evid = np.array([vmap[j] for j in evtags])
-        self.triangles = evid.reshape((self.triangles_tags.shape[-1],-1))
+        self.triangles = evid.reshape((self.triangles_tags.shape[-1], -1))
 
 
 def my_function(xyz):
-    a = 6*(np.hypot(xyz[...,0]-.5,xyz[...,1]-.5)-.2)
-    f = np.real(np.arctanh(a+0j))
+    a = 6 * (np.hypot(xyz[..., 0] - .5, xyz[..., 1] - .5) - .2)
+    f = np.real(np.arctanh(a + 0j))
     return f
 
 
 def compute_interpolation_error(nodes, triangles, f):
-    uvw,weights = gmsh.model.mesh.getIntegrationPoints(2,"Gauss2")
-    jac,det,pt = gmsh.model.mesh.getJacobians(2,uvw)
-    numcomp,sf = gmsh.model.mesh.getBasisFunctions(2,uvw,"Lagrange")
-    sf = sf.reshape((weights.shape[0],-1))
-    qx = pt.reshape((triangles.shape[0],-1,3))
-    det = np.abs(det.reshape((triangles.shape[0],-1)))
+    uvw, weights = gmsh.model.mesh.getIntegrationPoints(2, "Gauss2")
+    jac, det, pt = gmsh.model.mesh.getJacobians(2, uvw)
+    numcomp, sf = gmsh.model.mesh.getBasisFunctions(2, uvw, "Lagrange")
+    sf = sf.reshape((weights.shape[0], -1))
+    qx = pt.reshape((triangles.shape[0], -1, 3))
+    det = np.abs(det.reshape((triangles.shape[0], -1)))
     f_vert = f(nodes)
-    f_fem = np.dot(f_vert[triangles],sf)
-    err_tri = np.sum((f_fem-f(qx))**2*det*weights,1)
+    f_fem = np.dot(f_vert[triangles], sf)
+    err_tri = np.sum((f_fem - f(qx))**2 * det * weights, 1)
     return f_vert, np.sqrt(err_tri)
 
 
@@ -44,15 +43,20 @@ def compute_size_field(nodes, triangles, err, N):
     x = nodes[triangles]
     a = 2.
     d = 2.
-    fact = (a**((2.+a)/(1.+a)) + a**(1./(1.+a))) * np.sum(err**(2./(1.+a)))
-    ri = err**(2./(2.*(1+a))) * a**(1./(d*(1.+a))) * ((1.+a)*N/fact)**(1./d)
-    return triangle_max_edge(x)/ri
+    fact = (a**((2. + a) /
+                (1. + a)) + a**(1. / (1. + a))) * np.sum(err**(2. / (1. + a)))
+    ri = err**(2. /
+               (2. *
+                (1 + a))) * a**(1. /
+                                (d *
+                                 (1. + a))) * ((1. + a) * N / fact)**(1. / d)
+    return triangle_max_edge(x) / ri
 
 
 print("Usage: adapt_mesh [intial lc] [target #elements] [dump files]")
 
-lc = 0.02;
-N = 10000;
+lc = 0.02
+N = 10000
 dumpfiles = False
 if len(sys.argv) > 1: lc = float(sys.argv[1])
 if len(sys.argv) > 2: N = int(sys.argv[2])
@@ -67,28 +71,29 @@ square = gmsh.model.occ.addRectangle(0, 0, 0, 1, 1)
 gmsh.model.occ.synchronize()
 
 # create intial uniform mesh
-pnts = gmsh.model.getBoundary([(2,square)], True, True, True);
+pnts = gmsh.model.getBoundary([(2, square)], True, True, True)
 gmsh.model.mesh.setSize(pnts, lc)
 gmsh.model.mesh.generate(2)
 if dumpfiles: gmsh.write("mesh.msh")
 mesh = Mesh()
 
 # compute and visualize the interpolation error
-f_nod, err_ele = compute_interpolation_error(mesh.vxyz, mesh.triangles, my_function)
+f_nod, err_ele = compute_interpolation_error(mesh.vxyz, mesh.triangles,
+                                             my_function)
 f_view = gmsh.view.add("nodal function")
-gmsh.view.addModelData(f_view, 0, "square", "NodeData",
-                       mesh.vtags, f_nod[:,None])
+gmsh.view.addModelData(f_view, 0, "square", "NodeData", mesh.vtags,
+                       f_nod[:, None])
 if dumpfiles: gmsh.view.write(f_view, "f.pos")
 err_view = gmsh.view.add("element-wise error")
 gmsh.view.addModelData(err_view, 0, "square", "ElementData",
-                       mesh.triangles_tags, err_ele[:,None])
+                       mesh.triangles_tags, err_ele[:, None])
 if dumpfiles: gmsh.view.write(err_view, "err.pos")
 
 # compute and visualize the remeshing size field
-sf_ele = compute_size_field(mesh.vxyz,mesh.triangles, err_ele, N)
+sf_ele = compute_size_field(mesh.vxyz, mesh.triangles, err_ele, N)
 sf_view = gmsh.view.add("mesh size field")
 gmsh.view.addModelData(sf_view, 0, "square", "ElementData",
-                       mesh.triangles_tags, sf_ele[:,None])
+                       mesh.triangles_tags, sf_ele[:, None])
 if dumpfiles: gmsh.view.write(sf_view, "sf.pos")
 
 # create a new gmsh.model (to remesh the original gmsh.model in-place, the size field
@@ -98,22 +103,23 @@ gmsh.model.occ.addRectangle(0, 0, 0, 1, 1)
 gmsh.model.occ.synchronize()
 
 # mesh the new gmsh.model using the size field
-bg_field = gmsh.model.mesh.field.add("PostView");
-gmsh.model.mesh.field.setNumber(bg_field, "ViewTag", sf_view);
-gmsh.model.mesh.field.setAsBackgroundMesh(bg_field);
+bg_field = gmsh.model.mesh.field.add("PostView")
+gmsh.model.mesh.field.setNumber(bg_field, "ViewTag", sf_view)
+gmsh.model.mesh.field.setAsBackgroundMesh(bg_field)
 gmsh.model.mesh.generate(2)
 if dumpfiles: gmsh.write("mesh2.msh")
 mesh2 = Mesh()
 
 # compute and visualize the interpolation error on the adapted mesh
-f2_nod, err2_ele = compute_interpolation_error(mesh2.vxyz,mesh2.triangles, my_function)
+f2_nod, err2_ele = compute_interpolation_error(mesh2.vxyz, mesh2.triangles,
+                                               my_function)
 f2_view = gmsh.view.add("nodal function on adapted mesh")
-gmsh.view.addModelData(f2_view, 0, "square2", "NodeData",
-                       mesh2.vtags, f2_nod[:,None])
+gmsh.view.addModelData(f2_view, 0, "square2", "NodeData", mesh2.vtags,
+                       f2_nod[:, None])
 if dumpfiles: gmsh.view.write(f2_view, "f2.pos")
 err2_view = gmsh.view.add("element-wise error on adapated mesh")
 gmsh.view.addModelData(err2_view, 0, "square2", "ElementData",
-                       mesh2.triangles_tags, err2_ele[:,None])
+                       mesh2.triangles_tags, err2_ele[:, None])
 if dumpfiles: gmsh.view.write(err2_view, "err2.pos")
 
 # show everything in the gui
diff --git a/demos/api/boolean.py b/demos/api/boolean.py
index f47f59a0f20451ffb7f2173a4ffb00315a9e096d..847c2fa7213c4e09dc69c484a91cdccbd1990038 100644
--- a/demos/api/boolean.py
+++ b/demos/api/boolean.py
@@ -3,40 +3,38 @@
 import gmsh
 import sys
 
-model = gmsh.model
-factory = model.occ
-
 gmsh.initialize(sys.argv)
 
 gmsh.option.setNumber("General.Terminal", 1)
 
-model.add("boolean")
+gmsh.model.add("boolean")
 
 # from http://en.wikipedia.org/wiki/Constructive_solid_geometry
 
-gmsh.option.setNumber("Mesh.Algorithm", 6);
-gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.4);
-gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.4);
+gmsh.option.setNumber("Mesh.Algorithm", 6)
+gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.4)
+gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.4)
 
-R = 1.4; Rs = R*.7; Rt = R*1.25
+R = 1.4
+Rs = R * .7
+Rt = R * 1.25
 
-factory.addBox(-R,-R,-R, 2*R,2*R,2*R, 1)
-factory.addSphere(0,0,0,Rt, 2)
-factory.intersect([(3, 1)], [(3, 2)], 3)
-factory.addCylinder(-2*R,0,0, 4*R,0,0, Rs, 4)
-factory.addCylinder(0,-2*R,0, 0,4*R,0, Rs, 5)
-factory.addCylinder(0,0,-2*R, 0,0,4*R, Rs, 6)
-factory.fuse([(3, 4), (3, 5)], [(3, 6)], 7)
-factory.cut([(3, 3)], [(3, 7)], 8)
+gmsh.model.occ.addBox(-R, -R, -R, 2 * R, 2 * R, 2 * R, 1)
+gmsh.model.occ.addSphere(0, 0, 0, Rt, 2)
+gmsh.model.occ.intersect([(3, 1)], [(3, 2)], 3)
+gmsh.model.occ.addCylinder(-2 * R, 0, 0, 4 * R, 0, 0, Rs, 4)
+gmsh.model.occ.addCylinder(0, -2 * R, 0, 0, 4 * R, 0, Rs, 5)
+gmsh.model.occ.addCylinder(0, 0, -2 * R, 0, 0, 4 * R, Rs, 6)
+gmsh.model.occ.fuse([(3, 4), (3, 5)], [(3, 6)], 7)
+gmsh.model.occ.cut([(3, 3)], [(3, 7)], 8)
 
-factory.synchronize();
+gmsh.model.occ.synchronize()
 
-model.mesh.generate(3)
-#model.mesh.refine()
-#model.mesh.setOrder(2)
-#model.mesh.partition(4)
+gmsh.model.mesh.generate(3)
+#gmsh.model.mesh.refine()
+#gmsh.model.mesh.setOrder(2)
+#gmsh.model.mesh.partition(4)
 
 gmsh.write("boolean.msh")
 
 gmsh.finalize()
-
diff --git a/demos/api/crack.py b/demos/api/crack.py
index 8481c7b9a2b60e392c3a9fd15bf0609835a1835a..fa2a51efd9fc06bd1d7f7411c87f25e12f104358 100644
--- a/demos/api/crack.py
+++ b/demos/api/crack.py
@@ -18,7 +18,7 @@ pt3 = gmsh.model.occ.addPoint(0.6, 0.1, 0)
 pt4 = gmsh.model.occ.addPoint(0.1, 0.3, 0)
 line2 = gmsh.model.occ.addLine(pt3, pt4)
 
-o, m = gmsh.model.occ.fragment([(2,surf1)], [(1,line1), (1,line2)])
+o, m = gmsh.model.occ.fragment([(2, surf1)], [(1, line1), (1, line2)])
 gmsh.model.occ.synchronize()
 
 # m contains, for each input entity (surf1, line1 and line2), the child entities
diff --git a/demos/api/custom_gui.py b/demos/api/custom_gui.py
index 83f7ba3f4bb20b1ce976c047ccb7d85399b59cae..6cc40c991bac3c0f5c32af6a00f755f3c85be5ae 100644
--- a/demos/api/custom_gui.py
+++ b/demos/api/custom_gui.py
@@ -31,10 +31,11 @@ gmsh.onelab.set(parameters)
 # flag that will be set to interrupt a calculation
 stop_computation = False
 
+
 # a computationally expensive routine, that will be run in its own thread
 def compute(arg):
     iterations = gmsh.onelab.getNumber("My App/Iterations")
-    progress= gmsh.onelab.getNumber("My App/Show progress?")
+    progress = gmsh.onelab.getNumber("My App/Show progress?")
     n = int(iterations[0]) if len(iterations) > 0 else 1
     show = True if (len(progress) > 0 and progress[0] == 1) else False
     p = 0
@@ -44,7 +45,7 @@ def compute(arg):
         # stop computation if requested by clicking on "Stop it!"
         if stop_computation:
             break
-        k = math.sin(k) + math.cos(j/45.)
+        k = math.sin(k) + math.cos(j / 45.)
         # show progress in real time?
         if (show == 1) and (n > 1) and (not j % (n / 100)):
             p = p + 1
@@ -64,6 +65,7 @@ def compute(arg):
     gmsh.fltk.awake("update")
     return
 
+
 # create the graphical user interface
 gmsh.fltk.initialize()
 
@@ -89,7 +91,7 @@ while 1:
         n = int(gmsh.onelab.getNumber("My App/Number of threads")[0])
         for i in range(n):
             t = threading.Thread(target=compute,
-                                 args=("My App/Thread {0}".format(i + 1),))
+                                 args=("My App/Thread {0}".format(i + 1), ))
             t.start()
     elif action[0] == "should stop":
         stop_computation = True
diff --git a/demos/api/discrete.py b/demos/api/discrete.py
index 0044b99885b2b1006f9f69c860597c36c56954ce..2c89b1ff5092a3c9be6619af799f8030133b44c5 100644
--- a/demos/api/discrete.py
+++ b/demos/api/discrete.py
@@ -4,25 +4,45 @@ import sys
 gmsh.initialize(sys.argv)
 gmsh.option.setNumber("General.Terminal", 1)
 
-gmsh.model.add("test");
+gmsh.model.add("test")
 
 # add discrete surface with tag 1
 gmsh.model.addDiscreteEntity(2, 1)
 
 # add 4 mesh nodes
-gmsh.model.mesh.addNodes(2, 1,
-                         [1, 2, 3, 4], # node tags: 1, 2, 3, and 4
-                         [0., 0., 0., # coordinates of node 1
-                          1., 0., 0., # coordinates of node 2
-                          1., 1., 0., # ...
-                          0., 1., 0.])
+gmsh.model.mesh.addNodes(
+    2,
+    1,
+    [1, 2, 3, 4],  # node tags: 1, 2, 3, and 4
+    [
+        0.,
+        0.,
+        0.,  # coordinates of node 1
+        1.,
+        0.,
+        0.,  # coordinates of node 2
+        1.,
+        1.,
+        0.,  # ...
+        0.,
+        1.,
+        0.
+    ])
 
 # add 2 triangles
-gmsh.model.mesh.addElements(2, 1,
-                            [2], # single type : 3-node triangle
-                            [[1, 2]], # triangle tags: 1 and 2
-                            [[1, 2, 3, # triangle 1: nodes 1, 2, 3
-                              1, 3, 4]]) # triangle 2: nodes 1, 3, 4
+gmsh.model.mesh.addElements(
+    2,
+    1,
+    [2],  # single type : 3-node triangle
+    [[1, 2]],  # triangle tags: 1 and 2
+    [[
+        1,
+        2,
+        3,  # triangle 1: nodes 1, 2, 3
+        1,
+        3,
+        4
+    ]])  # triangle 2: nodes 1, 3, 4
 
 # export the mesh ; use explore.py to read and examine the mesh
 gmsh.write("test.msh")
diff --git a/demos/api/explore.py b/demos/api/explore.py
index 8352e331b462e563c01ca5762903da8f93854545..badd7ae73e9f2c013a6db0de6a3fc4aee1195b4a 100644
--- a/demos/api/explore.py
+++ b/demos/api/explore.py
@@ -28,10 +28,11 @@ for e in entities:
     print(" - boundary entities " + str(boundary))
     partitions = gmsh.model.getPartitions(e[0], e[1])
     if len(partitions):
-           print(" - Partition tag(s): " + str(partitions) +
-                 " - parent entity " + str(gmsh.model.getParent(e[0], e[1])))
+        print(" - Partition tag(s): " + str(partitions) + " - parent entity " +
+              str(gmsh.model.getParent(e[0], e[1])))
     for t in elemTypes:
-        name, dim, order, numv, parv, _ = gmsh.model.mesh.getElementProperties(t)
+        name, dim, order, numv, parv, _ = gmsh.model.mesh.getElementProperties(
+            t)
         print(" - Element type: " + name + ", order " + str(order) + " (" +
               str(numv) + " nodes in param coord: " + str(parv) + ")")
 
diff --git a/demos/api/flatten.py b/demos/api/flatten.py
index 71e50b3091ec6bd2c99bf2ffd5427d0612eb8ea8..b4f5789a946bfd5fb281e68ce8a3991a5b51e4c1 100644
--- a/demos/api/flatten.py
+++ b/demos/api/flatten.py
@@ -24,7 +24,8 @@ entities = gmsh.model.getEntities()
 # get the nodes and elements
 for e in entities:
     nodeTags[e], nodeCoords[e], _ = gmsh.model.mesh.getNodes(e[0], e[1])
-    elementTypes[e], elementTags[e], elementNodeTags[e] = gmsh.model.mesh.getElements(e[0], e[1])
+    elementTypes[e], elementTags[e], elementNodeTags[
+        e] = gmsh.model.mesh.getElements(e[0], e[1])
 
 # delete the mesh
 gmsh.model.mesh.clear()
@@ -34,7 +35,8 @@ for e in entities:
     for i in range(2, len(nodeCoords[e]), 3):
         nodeCoords[e][i] = 0
     gmsh.model.mesh.addNodes(e[0], e[1], nodeTags[e], nodeCoords[e])
-    gmsh.model.mesh.addElements(e[0], e[1], elementTypes[e], elementTags[e], elementNodeTags[e])
+    gmsh.model.mesh.addElements(e[0], e[1], elementTypes[e], elementTags[e],
+                                elementNodeTags[e])
 
 gmsh.fltk.run()
 #gmsh.write('flat.msh')
diff --git a/demos/api/glue_and_remesh_stl.py b/demos/api/glue_and_remesh_stl.py
index 4579146797b654dd7501c7258a5086016f7dab1d..856c6b04fbdb7d91cd3a8a03c4c1e9b54dedc428 100644
--- a/demos/api/glue_and_remesh_stl.py
+++ b/demos/api/glue_and_remesh_stl.py
@@ -16,7 +16,7 @@ gmsh.model.mesh.removeDuplicateNodes()
 
 # classify surface mesh according to given angle, and create discrete model
 # entities (surfaces, curves and points) accordingly
-gmsh.model.mesh.classifySurfaces(math.pi/2)
+gmsh.model.mesh.classifySurfaces(math.pi / 2)
 
 # Notes:
 #
@@ -42,9 +42,9 @@ gmsh.model.geo.addVolume([l])
 gmsh.model.geo.synchronize()
 
 # mesh
-gmsh.option.setNumber("Mesh.Algorithm", 6);
-gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.4);
-gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.4);
+gmsh.option.setNumber("Mesh.Algorithm", 6)
+gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.4)
+gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.4)
 gmsh.model.mesh.generate(3)
 
 gmsh.fltk.run()
diff --git a/demos/api/gui.py b/demos/api/gui.py
index 9151cdac1ebd32787867b8c412d88bcbabef5326..bc376d3c562b3004df512b72f3a7450b8783fe2b 100644
--- a/demos/api/gui.py
+++ b/demos/api/gui.py
@@ -1,9 +1,6 @@
 import gmsh
 import sys
 
-model = gmsh.model
-factory = model.occ
-
 gmsh.initialize(sys.argv)
 
 # creates the FLTK user interface; this could also be called after the geometry
@@ -11,26 +8,28 @@ gmsh.initialize(sys.argv)
 gmsh.fltk.initialize()
 
 # Copied from boolean.py...
-model.add("boolean")
-gmsh.option.setNumber("Mesh.Algorithm", 6);
-gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.4);
-gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.4);
-R = 1.4; Rs = R*.7; Rt = R*1.25
-factory.addBox(-R,-R,-R, 2*R,2*R,2*R, 1)
-factory.addSphere(0,0,0,Rt, 2)
-factory.intersect([(3, 1)], [(3, 2)], 3)
-factory.addCylinder(-2*R,0,0, 4*R,0,0, Rs, 4)
-factory.addCylinder(0,-2*R,0, 0,4*R,0, Rs, 5)
-factory.addCylinder(0,0,-2*R, 0,0,4*R, Rs, 6)
-factory.fuse([(3, 4), (3, 5)], [(3, 6)], 7)
-factory.cut([(3, 3)], [(3, 7)], 8)
-factory.synchronize();
+gmsh.model.add("boolean")
+gmsh.option.setNumber("Mesh.Algorithm", 6)
+gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.4)
+gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.4)
+R = 1.4
+Rs = R * .7
+Rt = R * 1.25
+gmsh.model.occ.addBox(-R, -R, -R, 2 * R, 2 * R, 2 * R, 1)
+gmsh.model.occ.addSphere(0, 0, 0, Rt, 2)
+gmsh.model.occ.intersect([(3, 1)], [(3, 2)], 3)
+gmsh.model.occ.addCylinder(-2 * R, 0, 0, 4 * R, 0, 0, Rs, 4)
+gmsh.model.occ.addCylinder(0, -2 * R, 0, 0, 4 * R, 0, Rs, 5)
+gmsh.model.occ.addCylinder(0, 0, -2 * R, 0, 0, 4 * R, Rs, 6)
+gmsh.model.occ.fuse([(3, 4), (3, 5)], [(3, 6)], 7)
+gmsh.model.occ.cut([(3, 3)], [(3, 7)], 8)
+gmsh.model.occ.synchronize()
 # ...end of copy
 
 # hide volume
-model.setVisibility(model.getEntities(3),0)
+gmsh.model.setVisibility(gmsh.model.getEntities(3), 0)
 # color all surfaces gold
-model.setColor(model.getEntities(2),249,166,2)
+gmsh.model.setColor(gmsh.model.getEntities(2), 249, 166, 2)
 
 # this would be equivalent to gmsh.fltk.run():
 #
diff --git a/demos/api/hex.py b/demos/api/hex.py
index 4ea3e3a5cad69c4f4049c3e1edd9cc85b77b4645..c5b809028ee747d39d9601701b25f80ff37093b5 100644
--- a/demos/api/hex.py
+++ b/demos/api/hex.py
@@ -1,11 +1,11 @@
 import gmsh
 gmsh.initialize()
-gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) # save in old MSH format
+gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)  # save in old MSH format
 N = 4
-Rec2d = True # tris or quads
-Rec3d = True # tets, prisms or hexas
-p = gmsh.model.geo.addPoint(0,0,0)
-l = gmsh.model.geo.extrude([(0,p)], 1, 0, 0, [N], [1])
+Rec2d = True  # tris or quads
+Rec3d = True  # tets, prisms or hexas
+p = gmsh.model.geo.addPoint(0, 0, 0)
+l = gmsh.model.geo.extrude([(0, p)], 1, 0, 0, [N], [1])
 s = gmsh.model.geo.extrude([l[1]], 0, 1, 0, [N], [1], recombine=Rec2d)
 v = gmsh.model.geo.extrude([s[1]], 0, 0, 1, [N], [1], recombine=Rec3d)
 gmsh.model.geo.synchronize()
diff --git a/demos/api/mesh_from_discrete_curve.py b/demos/api/mesh_from_discrete_curve.py
index b3f3daed7a46abcaa0428ccdb22a59e46267ebfd..5419903978a5c377b9e9ad15bb87b121b1c817de 100644
--- a/demos/api/mesh_from_discrete_curve.py
+++ b/demos/api/mesh_from_discrete_curve.py
@@ -5,18 +5,21 @@ import math
 gmsh.initialize(sys.argv)
 gmsh.option.setNumber("General.Terminal", 1)
 
-gmsh.model.add("2d surface mesh with purely discrete boundary");
+gmsh.model.add("2d surface mesh with purely discrete boundary")
 
 # create a discrete curve with N nodes and N line elements
 gmsh.model.addDiscreteEntity(1, 100)
 N = 50
-dt = 2*math.pi/N
+dt = 2 * math.pi / N
 pts = [[math.cos(i * dt), math.sin(i * dt), 0] for i in range(N)]
 flat_pts = [item for sublist in pts for item in sublist]
-gmsh.model.mesh.addNodes(1, 100, range(1, N+1), flat_pts)
-n = [item for sublist in [[i, i+1] for i in range(1, N+1)] for item in sublist]
+gmsh.model.mesh.addNodes(1, 100, range(1, N + 1), flat_pts)
+n = [
+    item for sublist in [[i, i + 1] for i in range(1, N + 1)]
+    for item in sublist
+]
 n[-1] = 1
-gmsh.model.mesh.addElements(1, 100, [1], [range(1, N+1)], [n])
+gmsh.model.mesh.addElements(1, 100, [1], [range(1, N + 1)], [n])
 
 # create a plane surface from the discrete curve
 gmsh.model.geo.addCurveLoop([100], 101)
diff --git a/demos/api/neighbors.py b/demos/api/neighbors.py
index e448f0e1a4128f24f582eafdbffddbac589d3e58..20db8f2013b34dda9d78d87a215ce95fc7ffe66c 100644
--- a/demos/api/neighbors.py
+++ b/demos/api/neighbors.py
@@ -7,8 +7,8 @@ import sys
 gmsh.initialize(sys.argv)
 gmsh.option.setNumber("General.Terminal", 1)
 
-gmsh.model.add("my test model");
-gmsh.model.occ.addBox(0,0,0, 1,1,1);
+gmsh.model.add("my test model")
+gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1)
 gmsh.model.occ.synchronize()
 gmsh.model.mesh.generate(3)
 
@@ -20,9 +20,9 @@ print "--- computing face x tet incidence"
 faces = []
 fxt = {}
 for i in range(0, len(fnodes), 3):
-    f = tuple(sorted(fnodes[i:i+3]))
+    f = tuple(sorted(fnodes[i:i + 3]))
     faces.append(f)
-    t = tets[i/12]
+    t = tets[i / 12]
     if not f in fxt:
         fxt[f] = [t]
     else:
@@ -32,7 +32,7 @@ print "--- computing neighbors by face"
 txt = {}
 for i in range(0, len(faces)):
     f = faces[i]
-    t = tets[i/4]
+    t = tets[i / 4]
     if not t in txt:
         txt[t] = set()
     for tt in fxt[f]:
diff --git a/demos/api/normals.py b/demos/api/normals.py
index 49b272d3433213cad0efe1fb8a2c3e3bee074122..70ab49ab4efdc290dd0ae2a34975557f9e63e825 100644
--- a/demos/api/normals.py
+++ b/demos/api/normals.py
@@ -26,24 +26,24 @@ for e in ent:
     normals = gmsh.model.getNormal(surf, param)
     # get surface curvature
     curv = gmsh.model.getCurvature(2, surf, param)
-    for i in range(0,len(coord),3):
+    for i in range(0, len(coord), 3):
         nn.append(coord[i])
-        nn.append(coord[i+1])
-        nn.append(coord[i+2])
+        nn.append(coord[i + 1])
+        nn.append(coord[i + 2])
         nn.append(normals[i])
-        nn.append(normals[i+1])
-        nn.append(normals[i+2])
+        nn.append(normals[i + 1])
+        nn.append(normals[i + 2])
         cc.append(coord[i])
-        cc.append(coord[i+1])
-        cc.append(coord[i+2])
-        cc.append(curv[i/3])
+        cc.append(coord[i + 1])
+        cc.append(coord[i + 2])
+        cc.append(curv[i / 3])
 
 t = gmsh.view.add("normals")
-gmsh.view.addListData(t, "VP", len(nn)/6, nn)
+gmsh.view.addListData(t, "VP", len(nn) / 6, nn)
 gmsh.view.write(t, "normals.pos")
 
 t = gmsh.view.add("curvatures")
-gmsh.view.addListData(t, "SP", len(cc)/4, cc)
+gmsh.view.addListData(t, "SP", len(cc) / 4, cc)
 gmsh.view.write(t, "curvatures.pos")
 
 gmsh.finalize()
diff --git a/demos/api/onelab_data.py b/demos/api/onelab_data.py
index 1c06f99ccd4b5ed747270f26755f4a12128ac9fd..1ac4e2782fd6b84101636305af2243be87defb00 100644
--- a/demos/api/onelab_data.py
+++ b/demos/api/onelab_data.py
@@ -11,7 +11,7 @@ gmsh.option.setNumber("General.Terminal", 1)
 gmsh.open(sys.argv[1])
 
 # attempts to run a client selected when opening the file (e.g. a .pro file)
-gmsh.onelab.run();
+gmsh.onelab.run()
 
 json = gmsh.onelab.get()
 print json
diff --git a/demos/api/onelab_test.py b/demos/api/onelab_test.py
index 784b647506cf4480b03cf38f5d5300a7a0becf5c..4c597a3c4798272bb2dfb6c59d048837b53c9920 100644
--- a/demos/api/onelab_test.py
+++ b/demos/api/onelab_test.py
@@ -35,7 +35,7 @@ gmsh.onelab.set("""
 
 # get the full parameter, store it as a python dict, and change an attribute
 p = json.loads(gmsh.onelab.get("check 1"))
-p["attributes"] = {"Highlight":"Blue"}
+p["attributes"] = {"Highlight": "Blue"}
 gmsh.onelab.set(json.dumps(p))
 
 # shorter way to just change the value, without json overhead
diff --git a/demos/api/open.py b/demos/api/open.py
index 1fa194b5d702b28af56f2509bdac68cbe4b39856..62ae0cc8e3af1d4030bcac5ff6f69e7d182e72bb 100644
--- a/demos/api/open.py
+++ b/demos/api/open.py
@@ -11,4 +11,3 @@ gmsh.open(sys.argv[1])
 gmsh.model.mesh.generate(3)
 gmsh.write("test.msh")
 gmsh.finalize()
-
diff --git a/demos/api/partition.py b/demos/api/partition.py
index 633f50eb844dd7c29b12c8870aed70ae9b45e602..d87540826f5e0f7bbe995df30626a5be34ea12f9 100644
--- a/demos/api/partition.py
+++ b/demos/api/partition.py
@@ -40,7 +40,8 @@ entities = gmsh.model.getEntities()
 for e in entities:
     partitions = gmsh.model.getPartitions(e[0], e[1])
     if len(partitions):
-        print("Entity " + str(e) + " of type " + gmsh.model.getType(e[0], e[1]))
+        print("Entity " + str(e) + " of type " +
+              gmsh.model.getType(e[0], e[1]))
         print(" - Partition(s): " + str(partitions))
         print(" - Parent: " + str(gmsh.model.getParent(e[0], e[1])))
         print(" - Boundary: " + str(gmsh.model.getBoundary([e])))
diff --git a/demos/api/periodic.py b/demos/api/periodic.py
index b9e6c7985df6cdacdc6d93f70cec41d48649272e..ef08e6c00e25693a73acfd50051ab34f9a75280f 100644
--- a/demos/api/periodic.py
+++ b/demos/api/periodic.py
@@ -8,15 +8,17 @@ R = 2
 gmsh.model.occ.addBox(0, 0, 0, R, R, R)
 gmsh.model.occ.synchronize()
 
-ent = gmsh.model.getEntities(0);
-gmsh.model.mesh.setSize(ent, 1);
-gmsh.model.mesh.setSize([(0,1)], 0.01);
-gmsh.model.mesh.setPeriodic(2, [2], [1], [1,0,0,R, 0,1,0,0, 0,0,1,0, 0,0,0,1])
+ent = gmsh.model.getEntities(0)
+gmsh.model.mesh.setSize(ent, 1)
+gmsh.model.mesh.setSize([(0, 1)], 0.01)
+gmsh.model.mesh.setPeriodic(2, [2], [1],
+                            [1, 0, 0, R, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
 
 gmsh.model.mesh.generate(2)
 #gmsh.model.mesh.setOrder(2)
 
-masterTag, nodeTags, nodeMasterTags, tfo = gmsh.model.mesh.getPeriodicNodes(2, 2)
+masterTag, nodeTags, nodeMasterTags, tfo = gmsh.model.mesh.getPeriodicNodes(
+    2, 2)
 print masterTag, nodeTags, nodeMasterTags, tfo
 
 gmsh.write("periodic.msh")
diff --git a/demos/api/plugin.py b/demos/api/plugin.py
index 373792209eb1dfdc38ac2ffa560c387eb3b877f9..99b7a1680557dda8030606d88d0482ad5872d4f5 100644
--- a/demos/api/plugin.py
+++ b/demos/api/plugin.py
@@ -5,23 +5,17 @@ gmsh.initialize(sys.argv)
 gmsh.option.setNumber("General.Terminal", 1)
 
 # Copied from discrete.py...
-gmsh.model.add("test");
+gmsh.model.add("test")
 gmsh.model.addDiscreteEntity(2, 1)
 gmsh.model.mesh.addNodes(2, 1, [1, 2, 3, 4],
-                         [0., 0., 0.,
-                          1., 0., 0.,
-                          1., 1., 0.,
-                          0., 1., 0.])
-gmsh.model.mesh.addElements(2, 1, [2], [[1, 2]],
-                            [[1, 2, 3,
-                              1, 3, 4]])
+                         [0., 0., 0., 1., 0., 0., 1., 1., 0., 0., 1., 0.])
+gmsh.model.mesh.addElements(2, 1, [2], [[1, 2]], [[1, 2, 3, 1, 3, 4]])
 # ... end of copy
 
 # create a view with some data
 t = gmsh.view.add("some data")
-gmsh.view.addModelData(t, 0, "test", "NodeData",
-                       [1, 2, 3, 4],
-                       [[1.],[10.],[20.],[1.]])
+gmsh.view.addModelData(t, 0, "test", "NodeData", [1, 2, 3, 4],
+                       [[1.], [10.], [20.], [1.]])
 
 # test getting data back
 dataType, tags, data, time, numComp = gmsh.view.getModelData(t, 0)
diff --git a/demos/api/poisson.py b/demos/api/poisson.py
index 5f68113d2d6600197a9486b7372325540ff4b0b8..5802d640f03e5b03318353b0bb8ba95cafa7d2c7 100644
--- a/demos/api/poisson.py
+++ b/demos/api/poisson.py
@@ -14,51 +14,59 @@ import sys
 DEBUG = 0
 RECOMBINE = 0
 
+
 def debug(*args):
     if not DEBUG: return
     if sys.version_info.major == 3:
-        exec("print( *args )")
+        exec ("print( *args )")
     else:
-        for item in args: exec("print item,")
-        exec("print")
+        for item in args:
+            exec ("print item,")
+        exec ("print")
+
 
 def error(argv):
     debug(argv)
     exit(1)
 
+
 def create_geometry():
-    model.add("poisson")
+    gmsh.model.add("poisson")
     surf = []
-    surf.append((2, factory.addRectangle(0, 0, 0, 1, 1)))
-    surf.append((2, factory.addDisk(0.7, 0.5, 0, 0.1, 0.1)))
-    surf.append((2, factory.addRectangle(0.2, 0.4, 0, 0.2, 0.2)))
-    surf, _ = factory.fragment(surf, [])
-    factory.synchronize()
-
-    bnd = model.getBoundary(surf, combined=True, oriented=False, recursive=False)
+    surf.append((2, gmsh.model.occ.addRectangle(0, 0, 0, 1, 1)))
+    surf.append((2, gmsh.model.occ.addDisk(0.7, 0.5, 0, 0.1, 0.1)))
+    surf.append((2, gmsh.model.occ.addRectangle(0.2, 0.4, 0, 0.2, 0.2)))
+    surf, _ = gmsh.model.occ.fragment(surf, [])
+    gmsh.model.occ.synchronize()
+
+    bnd = gmsh.model.getBoundary(surf,
+                                 combined=True,
+                                 oriented=False,
+                                 recursive=False)
     bnd = np.array(bnd)
 
-    model.addPhysicalGroup(2, [surf[0][1]], 2)
-    model.setPhysicalName(2, 2, 'SourceCircle')
-    model.addPhysicalGroup(2, [surf[1][1]], 3)
-    model.setPhysicalName(2, 3, 'SourceSquare')
-    model.addPhysicalGroup(2, [surf[2][1]], 4)
-    model.setPhysicalName(2, 4, 'Domain')
-    model.addPhysicalGroup(1, bnd[:,1], 11)
-    model.setPhysicalName(1, 11, 'Boundary')
+    gmsh.model.addPhysicalGroup(2, [surf[0][1]], 2)
+    gmsh.model.setPhysicalName(2, 2, 'SourceCircle')
+    gmsh.model.addPhysicalGroup(2, [surf[1][1]], 3)
+    gmsh.model.setPhysicalName(2, 3, 'SourceSquare')
+    gmsh.model.addPhysicalGroup(2, [surf[2][1]], 4)
+    gmsh.model.setPhysicalName(2, 4, 'Domain')
+    gmsh.model.addPhysicalGroup(1, bnd[:, 1], 11)
+    gmsh.model.setPhysicalName(1, 11, 'Boundary')
 
-    model.mesh.setSize(model.getEntities(0), 0.1);
+    gmsh.model.mesh.setSize(gmsh.model.getEntities(0), 0.1)
     return
 
+
 def fem_solve():
-    mshNodes = np.array(model.mesh.getNodes()[0])
+    mshNodes = np.array(gmsh.model.mesh.getNodes()[0])
     numMeshNodes = len(mshNodes)
     maxNodeTag = int(np.amax(mshNodes))
     debug('numMeshNodes =', numMeshNodes, ' maxNodeTag =', maxNodeTag)
 
     # typNodes[tag] = {0:does not exist, 1:internal node, 2:boundary node}
     # Existing node tags are defined here. Boundary node tag are identified later.
-    typNodes = np.zeros(maxNodeTag+1, dtype=np.int32)
+    typNodes = np.zeros(maxNodeTag + 1, dtype=np.int32)
     for tagNode in mshNodes:
         typNodes[tagNode] = 1
 
@@ -66,66 +74,70 @@ def fem_solve():
     # assembly
     matrowflat = np.array([], dtype=np.int32)
     matcolflat = np.array([], dtype=np.int32)
-    matflat    = np.array([], dtype=np.int32)
+    matflat = np.array([], dtype=np.int32)
     rhsrowflat = np.array([], dtype=np.int32)
-    rhsflat    = np.array([], dtype=np.int32)
+    rhsflat = np.array([], dtype=np.int32)
 
     # The read-in mesh is iterated over, looping successively (nested loops) over:
     # Physical groups/Geometrical entities/Element types/Elements
-    vGroups = model.getPhysicalGroups()
+    vGroups = gmsh.model.getPhysicalGroups()
     for iGroup in vGroups:
-        dimGroup = iGroup[0] # 1D, 2D or 3D
+        dimGroup = iGroup[0]  # 1D, 2D or 3D
         tagGroup = iGroup[1]
-        namGroup = model.getPhysicalName(dimGroup, tagGroup)
-        vEntities = model.getEntitiesForPhysicalGroup(dimGroup, tagGroup)
+        namGroup = gmsh.model.getPhysicalName(dimGroup, tagGroup)
+        vEntities = gmsh.model.getEntitiesForPhysicalGroup(dimGroup, tagGroup)
         for tagEntity in vEntities:
             # FIXME dimEntity should be optional when tagEntity given
             dimEntity = dimGroup
-            vElementTypes = model.mesh.getElementTypes(dimEntity,tagEntity)
+            vElementTypes = gmsh.model.mesh.getElementTypes(
+                dimEntity, tagEntity)
             for elementType in vElementTypes:
-                vTags, vNodes = model.mesh.getElementsByType(elementType, tagEntity)
+                vTags, vNodes = gmsh.model.mesh.getElementsByType(
+                    elementType, tagEntity)
                 numElements = len(vTags)
                 numGroupNodes = len(vNodes)
-                enode = np.array(vNodes, dtype=np.int32).reshape((numElements,-1))
+                enode = np.array(vNodes, dtype=np.int32).reshape(
+                    (numElements, -1))
                 numElementNodes = enode.shape[1]
                 debug('\nIn group', namGroup, ' with tag ', tagGroup,
                       ', numElements = e =', numElements)
-                debug('numGroupNodes =', numGroupNodes,', numElementNodes = n =',
-                      numElementNodes)
+                debug('numGroupNodes =', numGroupNodes,
+                      ', numElementNodes = n =', numElementNodes)
                 debug('%enodes (e,n) =', enode.shape)
 
                 # Assembly of stiffness matrix for all 2 dimensional elements
                 # (triangles or quadrangles)
-                if dimEntity==2 :
-                    prop = model.mesh.getElementProperties(elementType)
-                    uvw, weights = model.mesh.getIntegrationPoints(
+                if dimEntity == 2:
+                    prop = gmsh.model.mesh.getElementProperties(elementType)
+                    uvw, weights = gmsh.model.mesh.getIntegrationPoints(
                         elementType, "Gauss" + str(2 * prop[2]))
-                    numcomp, sf = model.mesh.getBasisFunctions(
+                    numcomp, sf, _ = gmsh.model.mesh.getBasisFunctions(
                         elementType, uvw, 'Lagrange')
                     weights = np.array(weights)
                     numGaussPoints = len(weights)
-                    debug('numGaussPoints = g =', numGaussPoints, ', %weights (g) =',
-                          weights.shape)
+                    debug('numGaussPoints = g =', numGaussPoints,
+                          ', %weights (g) =', weights.shape)
                     sf = np.array(sf).reshape((numGaussPoints, -1))
                     debug('%sf (g,n) =', sf.shape)
                     if sf.shape[1] != numElementNodes:
                         error('Something went wrong')
-                    numcomp, dsfdu = model.mesh.getBasisFunctions(
+                    numcomp, dsfdu, _ = gmsh.model.mesh.getBasisFunctions(
                         elementType, uvw, 'GradLagrange')
                     # remove useless dsfdw
                     dsfdu = np.array(dsfdu).reshape(
-                        (numGaussPoints,numElementNodes,3))[:,:,:-1]
+                        (numGaussPoints, numElementNodes, 3))[:, :, :-1]
                     debug('%dsfdu (g,n,u) =', dsfdu.shape)
 
-                    qjac, qdet, qpoint = model.mesh.getJacobians(
+                    qjac, qdet, qpoint = gmsh.model.mesh.getJacobians(
                         elementType, uvw, tagEntity)
                     debug('Gauss integr:', len(qjac), len(qdet), len(qpoint),
-                           '= (9, 1, 3) x', numGaussPoints, 'x', numElements)
-                    qdet = np.array(qdet).reshape((numElements, numGaussPoints))
+                          '= (9, 1, 3) x', numGaussPoints, 'x', numElements)
+                    qdet = np.array(qdet).reshape(
+                        (numElements, numGaussPoints))
                     debug('%qdet (e,g) =', qdet.shape)
                     # remove components of dxdu useless in dimEntity dimensions (here 2D)
                     dxdu = np.array(qjac).reshape(
-                        (numElements, numGaussPoints, 3, 3))[:,:,:-1,:-1]
+                        (numElements, numGaussPoints, 3, 3))[:, :, :-1, :-1]
                     # jacobian stored by row, so dxdu[i][j] = dxdu_ij = dxi/duj
                     debug('%dxdu (e,g,x,u)=', dxdu.shape)
 
@@ -133,17 +145,21 @@ def fem_solve():
                     dudx = np.linalg.inv(dxdu)
                     debug('%dudx (e,g,u,x) =', dudx.shape)
                     # sum over u = dot product
-                    dsfdx  = np.einsum("egxu,gnu->egnx",dudx,dsfdu);
+                    dsfdx = np.einsum("egxu,gnu->egnx", dudx, dsfdu)
                     debug('%dsfdx (e,g,n,x) =', dsfdx.shape)
                     # Gauss integration
-                    localmat = np.einsum("egik,egjk,eg,g->eij", dsfdx, dsfdx, qdet,
-                                         weights)
+                    localmat = np.einsum("egik,egjk,eg,g->eij", dsfdx, dsfdx,
+                                         qdet, weights)
                     debug('%localmat (e,n,n) =', localmat.shape)
 
                     # The next two lines are rather obscure. See explanations at
                     # the bottom of the file.
-                    matcol = np.repeat(enode[:,:,None],numElementNodes,axis=2)
-                    matrow = np.repeat(enode[:,None,:],numElementNodes,axis=1)
+                    matcol = np.repeat(enode[:, :, None],
+                                       numElementNodes,
+                                       axis=2)
+                    matrow = np.repeat(enode[:, None, :],
+                                       numElementNodes,
+                                       axis=1)
 
                     matcolflat = np.append(matcolflat, matcol.flatten())
                     matrowflat = np.append(matrowflat, matrow.flatten())
@@ -154,7 +170,8 @@ def fem_solve():
                             load = -1
                         else:
                             load = 1
-                        localrhs = load * np.einsum("gn,eg,g->en", sf, qdet, weights)
+                        localrhs = load * np.einsum("gn,eg,g->en", sf, qdet,
+                                                    weights)
                         rhsrowflat = np.append(rhsrowflat, enode.flatten())
                         rhsflat = np.append(rhsflat, localrhs.flatten())
 
@@ -166,23 +183,23 @@ def fem_solve():
     # Associate to all mesh nodes a line number in the system matrix, reserving
     # top lines for internal nodes and bottom lines for fixed nodes (boundary
     # nodes).
-    node2unknown = np.zeros(maxNodeTag+1, dtype=np.int32)
+    node2unknown = np.zeros(maxNodeTag + 1, dtype=np.int32)
     index = 0
-    for tagNode,typ in enumerate(typNodes):
-        if  typ == 1: # not fixed
+    for tagNode, typ in enumerate(typNodes):
+        if typ == 1:  # not fixed
             index += 1
             node2unknown[tagNode] = index
     numUnknowns = index
     debug('numUnknowns =', numUnknowns)
-    for tagNode,typ in enumerate(typNodes):
-        if  typ == 2: # fixed
+    for tagNode, typ in enumerate(typNodes):
+        if typ == 2:  # fixed
             index += 1
             node2unknown[tagNode] = index
 
     if index != numMeshNodes:
         error('Something went wrong')
 
-    unknown2node = np.zeros(numMeshNodes+1, dtype=np.int32)
+    unknown2node = np.zeros(numMeshNodes + 1, dtype=np.int32)
     for node, unkn in enumerate(node2unknown):
         unknown2node[unkn] = node
 
@@ -191,39 +208,38 @@ def fem_solve():
     # docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html
     # 'node2unknown-1' are because python numbers rows and columns from 0
     globalmat = scipy.sparse.coo_matrix(
-        (matflat, (node2unknown[matcolflat]-1,node2unknown[matrowflat]-1) ),
+        (matflat,
+         (node2unknown[matcolflat] - 1, node2unknown[matrowflat] - 1)),
         shape=(numMeshNodes, numMeshNodes)).tocsr()
 
     globalrhs = np.zeros(numMeshNodes)
-    for index,node in enumerate(rhsrowflat):
-        globalrhs[node2unknown[node]-1] += rhsflat[index]
+    for index, node in enumerate(rhsrowflat):
+        globalrhs[node2unknown[node] - 1] += rhsflat[index]
 
     debug('%globalmat =', globalmat.shape, ' %globalrhs =', globalrhs.shape)
 
     # Solve linear system Ax=b
-    sol = scipy.sparse.linalg.spsolve(globalmat[:numUnknowns,:numUnknowns],
+    sol = scipy.sparse.linalg.spsolve(globalmat[:numUnknowns, :numUnknowns],
                                       globalrhs[:numUnknowns])
-    sol = np.append(sol, np.zeros(numMeshNodes-numUnknowns))
+    sol = np.append(sol, np.zeros(numMeshNodes - numUnknowns))
     debug('%sol =', sol.shape)
 
     # Export solution
     sview = gmsh.view.add("solution")
     gmsh.view.addModelData(sview, 0, "", "NodeData", unknown2node[1:],
-                           sol[:,None])
+                           sol[:, None])
     return
 
 
-model = gmsh.model
-factory = model.occ
 gmsh.initialize(sys.argv)
 create_geometry()
 
 if RECOMBINE:
-    model.mesh.setRecombine(2,2)
-    model.mesh.setRecombine(2,3)
-    model.mesh.setRecombine(2,4)
+    gmsh.model.mesh.setRecombine(2, 2)
+    gmsh.model.mesh.setRecombine(2, 3)
+    gmsh.model.mesh.setRecombine(2, 4)
 
-model.mesh.generate(2)
+gmsh.model.mesh.generate(2)
 
 if DEBUG:
     gmsh.write('poisson.msh')
diff --git a/demos/api/remesh_stl.py b/demos/api/remesh_stl.py
index 758638212d3b89503588ff2e45d895694b87277f..57e68e35e869f38719c314a8aaa79572a639742e 100644
--- a/demos/api/remesh_stl.py
+++ b/demos/api/remesh_stl.py
@@ -5,9 +5,9 @@ import os
 gmsh.initialize()
 
 gmsh.option.setNumber("General.Terminal", 1)
-gmsh.option.setNumber("Mesh.Algorithm", 6);
-gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.75);
-gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.75);
+gmsh.option.setNumber("Mesh.Algorithm", 6)
+gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.75)
+gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.75)
 
 # load STL file
 path = os.path.dirname(os.path.abspath(__file__))
@@ -15,7 +15,7 @@ gmsh.merge(os.path.join(path, "object.stl"))
 
 # split input surface mesh based on an angle threshold of 40 degrees between
 # triangles, and generate new discrete surfaces suitable for reparametrization
-gmsh.model.mesh.classifySurfaces(40*math.pi/180., True, True)
+gmsh.model.mesh.classifySurfaces(40 * math.pi / 180., True, True)
 
 # create a geometry (through reparametrization) for all discrete curves and
 # discrete surfaces
diff --git a/demos/api/spherical_surf.py b/demos/api/spherical_surf.py
index 5ba950f9a87faa2fcc831a1f5ec5d42fcfb97fc1..ce28deb1d8a1b4969e5e00a1b936ad829e15cf8b 100644
--- a/demos/api/spherical_surf.py
+++ b/demos/api/spherical_surf.py
@@ -6,14 +6,14 @@ gmsh.initialize(sys.argv)
 gmsh.model.add("sphere_cut")
 R = 1
 R1 = 0.95
-sph = gmsh.model.occ.addSphere(0,0,0, R, -1, 0, math.pi/2, math.pi/2)
-b1 = gmsh.model.occ.addBox(R1,0,0, R,R,R)
-b2 = gmsh.model.occ.addBox(0,R1,0, R,R,R)
-b3 = gmsh.model.occ.addBox(0,0,R1, R,R,R)
-gmsh.model.occ.cut([(3,sph)], [(3,b1), (3,b2), (3,b3)])
+sph = gmsh.model.occ.addSphere(0, 0, 0, R, -1, 0, math.pi / 2, math.pi / 2)
+b1 = gmsh.model.occ.addBox(R1, 0, 0, R, R, R)
+b2 = gmsh.model.occ.addBox(0, R1, 0, R, R, R)
+b3 = gmsh.model.occ.addBox(0, 0, R1, R, R, R)
+gmsh.model.occ.cut([(3, sph)], [(3, b1), (3, b2), (3, b3)])
 gmsh.model.occ.synchronize()
 
-gmsh.model.removeEntities([(3,sph)])
-gmsh.model.removeEntities([(2,2), (2,4), (2,6)], True)
+gmsh.model.removeEntities([(3, sph)])
+gmsh.model.removeEntities([(2, 2), (2, 4), (2, 6)], True)
 gmsh.fltk.run()
 gmsh.finalize()
diff --git a/demos/api/spline.py b/demos/api/spline.py
index 522ef3f5cfc40ac24bad2965b9be67e377691a9d..a1b78520d05504310c47b12a836e69c2827ef061 100644
--- a/demos/api/spline.py
+++ b/demos/api/spline.py
@@ -1,41 +1,39 @@
 import gmsh
 import math
 
-model = gmsh.model
-factory = model.occ
-
 gmsh.initialize()
 gmsh.option.setNumber("General.Terminal", 1)
 
-model.add("spline")
+gmsh.model.add("spline")
 
 for i in range(1, 11):
-    factory.addPoint(i, math.sin(i/9.*2.*math.pi), 0, 0.1, i)
+    gmsh.model.occ.addPoint(i, math.sin(i / 9. * 2. * math.pi), 0, 0.1, i)
 
-factory.addSpline(range(1, 11), 1)
-factory.addBSpline(range(1, 11), 2)
-factory.addBezier(range(1, 11), 3)
+gmsh.model.occ.addSpline(range(1, 11), 1)
+gmsh.model.occ.addBSpline(range(1, 11), 2)
+gmsh.model.occ.addBezier(range(1, 11), 3)
 
-factory.addPoint(0.2,-1.6,0,0.1, 101) 
-factory.addPoint(1.2,-1.6,0,0.1, 102) 
-factory.addPoint(1.2,-1.1,0,0.1, 103) 
-factory.addPoint(0.3,-1.1,0,0.1, 104) 
-factory.addPoint(0.7,-1,0,0.1, 105) 
+gmsh.model.occ.addPoint(0.2, -1.6, 0, 0.1, 101)
+gmsh.model.occ.addPoint(1.2, -1.6, 0, 0.1, 102)
+gmsh.model.occ.addPoint(1.2, -1.1, 0, 0.1, 103)
+gmsh.model.occ.addPoint(0.3, -1.1, 0, 0.1, 104)
+gmsh.model.occ.addPoint(0.7, -1, 0, 0.1, 105)
 
 # periodic bspline through the control points
-factory.addSpline([103,102,101,104,105,103], 100)
+gmsh.model.occ.addSpline([103, 102, 101, 104, 105, 103], 100)
 
 # periodic bspline from given control points and default parameters - will
 # create a new vertex
-factory.addBSpline([103,102,101,104,105,103], 101)
+gmsh.model.occ.addBSpline([103, 102, 101, 104, 105, 103], 101)
 
 # general bspline with explicit degree, knots and multiplicities
-factory.addPoint(0,-2,0,0.1,201)
-factory.addPoint(1,-2,0,0.1,202)
-factory.addPoint(1,-3,0,0.1,203)
-factory.addPoint(0,-3,0,0.1,204)
-factory.addBSpline([201,202,203,204], 200, 2, [], [0,0.5,1], [3,1,3])
-
-factory.synchronize()
+gmsh.model.occ.addPoint(0, -2, 0, 0.1, 201)
+gmsh.model.occ.addPoint(1, -2, 0, 0.1, 202)
+gmsh.model.occ.addPoint(1, -3, 0, 0.1, 203)
+gmsh.model.occ.addPoint(0, -3, 0, 0.1, 204)
+gmsh.model.occ.addBSpline([201, 202, 203, 204], 200, 2, [], [0, 0.5, 1],
+                          [3, 1, 3])
+
+gmsh.model.occ.synchronize()
 gmsh.fltk.run()
 gmsh.finalize()
diff --git a/demos/api/step_assembly.py b/demos/api/step_assembly.py
index 4a27c99294943e26cc151bf61fc257cbb1f616a5..e73a600a8418dc3e07a8f67a6236d2c5eccdd902 100644
--- a/demos/api/step_assembly.py
+++ b/demos/api/step_assembly.py
@@ -18,7 +18,7 @@ for e in ent:
               str(gmsh.model.occ.getMass(e[0], e[1])) + ')')
         path = n.split('/')
         if e[0] == 3 and len(path) > 3:
-            if(path[2] not in physicals):
+            if (path[2] not in physicals):
                 physicals[path[2]] = []
             physicals[path[2]].append(e[1])
 
diff --git a/demos/api/step_boundary_colors.py b/demos/api/step_boundary_colors.py
index 2f73e87de57975a6065b9fc8b4f20be9f09fd501..22012eb288afe53a7d108137ec019be85305a87c 100644
--- a/demos/api/step_boundary_colors.py
+++ b/demos/api/step_boundary_colors.py
@@ -10,4 +10,4 @@ gmsh.merge(os.path.join(path, 'step_boundary_colors.stp'))
 
 for tag in gmsh.model.getEntities():
     col = gmsh.model.getColor(*tag)
-    if col != (0,0,255,0): print('entity', tag, 'color', col)
+    if col != (0, 0, 255, 0): print('entity', tag, 'color', col)
diff --git a/demos/api/terrain.py b/demos/api/terrain.py
index 75a21ac4be9aef1b8ca43ad70bb93dfe9713392b..3e2746301e37dbcfcdfa07f8210e1d37f78cf001 100644
--- a/demos/api/terrain.py
+++ b/demos/api/terrain.py
@@ -5,40 +5,47 @@ import sys
 gmsh.initialize(sys.argv)
 gmsh.option.setNumber('General.Terminal', 1)
 
-gmsh.model.add("terrain");
+gmsh.model.add("terrain")
 
 # create the terrain surface from N x N input data points (here simulated using
 # a simple function):
 N = 100
-coords = [] # x, y, z coordinates of all the points
-nodes = [] # tags of corresponding nodes
-tris = [] # connectivities (node tags) of triangle elements
-lin = [[],[],[],[]] # connectivities of boundary line elements
+coords = []  # x, y, z coordinates of all the points
+nodes = []  # tags of corresponding nodes
+tris = []  # connectivities (node tags) of triangle elements
+lin = [[], [], [], []]  # connectivities of boundary line elements
+
+
 def tag(i, j):
-    return (N+1) * i + j + 1
+    return (N + 1) * i + j + 1
+
+
 for i in range(N + 1):
     for j in range(N + 1):
         nodes.append(tag(i, j))
-        coords.extend([float(i)/N, float(j)/N, 0.05 * math.sin(10*float(i+j) / N)])
+        coords.extend([
+            float(i) / N,
+            float(j) / N, 0.05 * math.sin(10 * float(i + j) / N)
+        ])
         if i > 0 and j > 0:
-            tris.extend([tag(i-1, j-1), tag(i, j-1), tag(i-1, j)])
-            tris.extend([tag(i, j-1), tag(i, j), tag(i-1, j)])
+            tris.extend([tag(i - 1, j - 1), tag(i, j - 1), tag(i - 1, j)])
+            tris.extend([tag(i, j - 1), tag(i, j), tag(i - 1, j)])
         if (i == 0 or i == N) and j > 0:
-            lin[3 if i == 0 else 1].extend([tag(i, j-1), tag(i, j)])
+            lin[3 if i == 0 else 1].extend([tag(i, j - 1), tag(i, j)])
         if (j == 0 or j == N) and i > 0:
-            lin[0 if j == 0 else 2].extend([tag(i-1, j), tag(i, j)])
-pnt = [tag(0, 0), tag(N, 0), tag(N, N), tag(0, N)] # corner points element
+            lin[0 if j == 0 else 2].extend([tag(i - 1, j), tag(i, j)])
+pnt = [tag(0, 0), tag(N, 0), tag(N, N), tag(0, N)]  # corner points element
 
 # create 4 corner points
-gmsh.model.geo.addPoint(0, 0, coords[3*tag(0,0)-1], 1)
-gmsh.model.geo.addPoint(1, 0, coords[3*tag(N,0)-1], 2)
-gmsh.model.geo.addPoint(1, 1, coords[3*tag(N,N)-1], 3)
-gmsh.model.geo.addPoint(0, 1, coords[3*tag(0,N)-1], 4)
+gmsh.model.geo.addPoint(0, 0, coords[3 * tag(0, 0) - 1], 1)
+gmsh.model.geo.addPoint(1, 0, coords[3 * tag(N, 0) - 1], 2)
+gmsh.model.geo.addPoint(1, 1, coords[3 * tag(N, N) - 1], 3)
+gmsh.model.geo.addPoint(0, 1, coords[3 * tag(0, N) - 1], 4)
 gmsh.model.geo.synchronize()
 
 # create 4 discrete bounding curves, with their boundary points
 for i in range(4):
-    gmsh.model.addDiscreteEntity(1, i+1, [i+1, i+2 if i < 3 else 1])
+    gmsh.model.addDiscreteEntity(1, i + 1, [i + 1, i + 2 if i < 3 else 1])
 
 # create one discrete surface, with its bounding curves
 gmsh.model.addDiscreteEntity(2, 1, [1, 2, -3, -4])
@@ -100,31 +107,31 @@ c16 = gmsh.model.geo.addLine(3, p7)
 c17 = gmsh.model.geo.addLine(4, p8)
 
 # bottom and top
-ll1 = gmsh.model.geo.addCurveLoop([c1,c2,c3,c4])
+ll1 = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4])
 s1 = gmsh.model.geo.addPlaneSurface([ll1])
-ll2 = gmsh.model.geo.addCurveLoop([c5,c6,c7,c8])
+ll2 = gmsh.model.geo.addCurveLoop([c5, c6, c7, c8])
 s2 = gmsh.model.geo.addPlaneSurface([ll2])
 
 # lower
-ll3 = gmsh.model.geo.addCurveLoop([c1,c11,-1,-c10])
+ll3 = gmsh.model.geo.addCurveLoop([c1, c11, -1, -c10])
 s3 = gmsh.model.geo.addPlaneSurface([ll3])
-ll4 = gmsh.model.geo.addCurveLoop([c2,c12,-2,-c11])
+ll4 = gmsh.model.geo.addCurveLoop([c2, c12, -2, -c11])
 s4 = gmsh.model.geo.addPlaneSurface([ll4])
-ll5 = gmsh.model.geo.addCurveLoop([c3,c13,3,-c12])
+ll5 = gmsh.model.geo.addCurveLoop([c3, c13, 3, -c12])
 s5 = gmsh.model.geo.addPlaneSurface([ll5])
-ll6 = gmsh.model.geo.addCurveLoop([c4,c10,4,-c13])
+ll6 = gmsh.model.geo.addCurveLoop([c4, c10, 4, -c13])
 s6 = gmsh.model.geo.addPlaneSurface([ll6])
 sl1 = gmsh.model.geo.addSurfaceLoop([s1, s3, s4, s5, s6, 1])
 v1 = gmsh.model.geo.addVolume([sl1])
 
 # upper
-ll7 = gmsh.model.geo.addCurveLoop([c5,-c15,-1,c14])
+ll7 = gmsh.model.geo.addCurveLoop([c5, -c15, -1, c14])
 s7 = gmsh.model.geo.addPlaneSurface([ll7])
-ll8 = gmsh.model.geo.addCurveLoop([c6,-c16,-2,c15])
+ll8 = gmsh.model.geo.addCurveLoop([c6, -c16, -2, c15])
 s8 = gmsh.model.geo.addPlaneSurface([ll8])
-ll9 = gmsh.model.geo.addCurveLoop([c7,-c17,3,c16])
+ll9 = gmsh.model.geo.addCurveLoop([c7, -c17, 3, c16])
 s9 = gmsh.model.geo.addPlaneSurface([ll9])
-ll10 = gmsh.model.geo.addCurveLoop([c8,-c14,4,c17])
+ll10 = gmsh.model.geo.addCurveLoop([c8, -c14, 4, c17])
 s10 = gmsh.model.geo.addPlaneSurface([ll10])
 sl2 = gmsh.model.geo.addSurfaceLoop([s2, s7, s8, s9, s10, 1])
 v2 = gmsh.model.geo.addVolume([sl2])
diff --git a/demos/api/terrain_stl.py b/demos/api/terrain_stl.py
index f67e55b6259b02dba68fc913ec8759843d0a637c..2d996cc9ec306958e6b3bfab37b47aa77ca64c4b 100644
--- a/demos/api/terrain_stl.py
+++ b/demos/api/terrain_stl.py
@@ -14,7 +14,7 @@ gmsh.merge(os.path.join(path, 'terrain_stl_data.stl'))
 # classify the surface mesh according to given angle, and create discrete model
 # entities (surfaces, curves and points) accordingly; curveAngle forces bounding
 # curves to be split on sharp corners
-gmsh.model.mesh.classifySurfaces(math.pi, curveAngle=math.pi/3)
+gmsh.model.mesh.classifySurfaces(math.pi, curveAngle=math.pi / 3)
 
 # create a geometry for the discrete curves and surfaces
 gmsh.model.mesh.createGeometry()
@@ -23,7 +23,7 @@ gmsh.model.mesh.createGeometry()
 s = gmsh.model.getEntities(2)
 c = gmsh.model.getBoundary(s)
 
-if(len(c) != 4):
+if (len(c) != 4):
     gmsh.logger.write('Should have 4 boundary curves!', level='error')
 
 z = -1000
@@ -51,16 +51,16 @@ c11 = gmsh.model.geo.addLine(p2, p[1])
 c12 = gmsh.model.geo.addLine(p3, p[2])
 c13 = gmsh.model.geo.addLine(p4, p[3])
 
-ll1 = gmsh.model.geo.addCurveLoop([c1,c2,c3,c4])
+ll1 = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4])
 s1 = gmsh.model.geo.addPlaneSurface([ll1])
 
-ll3 = gmsh.model.geo.addCurveLoop([c1,c11,-c[0][1],-c10])
+ll3 = gmsh.model.geo.addCurveLoop([c1, c11, -c[0][1], -c10])
 s3 = gmsh.model.geo.addPlaneSurface([ll3])
-ll4 = gmsh.model.geo.addCurveLoop([c2,c12,-c[1][1],-c11])
+ll4 = gmsh.model.geo.addCurveLoop([c2, c12, -c[1][1], -c11])
 s4 = gmsh.model.geo.addPlaneSurface([ll4])
-ll5 = gmsh.model.geo.addCurveLoop([c3,c13,-c[2][1],-c12])
+ll5 = gmsh.model.geo.addCurveLoop([c3, c13, -c[2][1], -c12])
 s5 = gmsh.model.geo.addPlaneSurface([ll5])
-ll6 = gmsh.model.geo.addCurveLoop([c4,c10,-c[3][1],-c13])
+ll6 = gmsh.model.geo.addCurveLoop([c4, c10, -c[3][1], -c13])
 s6 = gmsh.model.geo.addPlaneSurface([ll6])
 sl1 = gmsh.model.geo.addSurfaceLoop([s1, s3, s4, s5, s6, s[0][1]])
 v1 = gmsh.model.geo.addVolume([sl1])
@@ -84,7 +84,6 @@ else:
     gmsh.option.setNumber('Mesh.CharacteristicLengthMin', 100)
     gmsh.option.setNumber('Mesh.CharacteristicLengthMax', 100)
 
-
 gmsh.fltk.run()
 
 gmsh.finalize()
diff --git a/demos/api/test.py b/demos/api/test.py
index a64951821d409f3f1aa87fa2cab7df0d4aaeffc6..89767d3fad2cbc8867be8c1ed49afefda37988e5 100644
--- a/demos/api/test.py
+++ b/demos/api/test.py
@@ -7,52 +7,50 @@ print(gmsh.option.getNumber("Mesh.Algorithm"))
 
 gmsh.open("square.msh")
 
-model = gmsh.model
-
-model.add("square")
-factory = model.geo
-factory.addPoint(0,0,0,0.1,1)
-factory.addPoint(1,0,0,0.1,2)
-factory.addPoint(1,1,0,0.1,3)
-factory.addPoint(0,1,0,0.1,4)
-factory.addLine(1,2,1)
-factory.addLine(2,3,2)
-factory.addLine(3,4,3)
-line4 = factory.addLine(4,1)
+gmsh.model.add("square")
+gmsh.model.geo.addPoint(0, 0, 0, 0.1, 1)
+gmsh.model.geo.addPoint(1, 0, 0, 0.1, 2)
+gmsh.model.geo.addPoint(1, 1, 0, 0.1, 3)
+gmsh.model.geo.addPoint(0, 1, 0, 0.1, 4)
+gmsh.model.geo.addLine(1, 2, 1)
+gmsh.model.geo.addLine(2, 3, 2)
+gmsh.model.geo.addLine(3, 4, 3)
+line4 = gmsh.model.geo.addLine(4, 1)
 print("line4 received tag ", line4)
-factory.addCurveLoop([1,2,3,line4],1)
-factory.addPlaneSurface([1],6)
-factory.synchronize()
+gmsh.model.geo.addCurveLoop([1, 2, 3, line4], 1)
+gmsh.model.geo.addPlaneSurface([1], 6)
+gmsh.model.geo.synchronize()
 
-ptag = model.addPhysicalGroup(1,[1,2,3,4])
-ent = model.getEntitiesForPhysicalGroup(1,ptag)
-print("new physical group ",ptag,":",ent, type(ent))
+ptag = gmsh.model.addPhysicalGroup(1, [1, 2, 3, 4])
+ent = gmsh.model.getEntitiesForPhysicalGroup(1, ptag)
+print("new physical group ", ptag, ":", ent, type(ent))
 
-model.addPhysicalGroup(2,[6])
+gmsh.model.addPhysicalGroup(2, [6])
 
 print(gmsh.option.getString("General.BuildOptions"))
 print(gmsh.option.getNumber("Mesh.Algorithm"))
 gmsh.option.setNumber("Mesh.Algorithm", 3.0)
 print(gmsh.option.getNumber("Mesh.Algorithm"))
-model.mesh.generate(2)
+gmsh.model.mesh.generate(2)
 
 gmsh.write("square.msh")
 
 print("Entities")
-entities = model.getEntities()
-for e in entities :
-    print("entity ",e)
-    types,tags,nodes = model.mesh.getElements(e[0],e[1])
+entities = gmsh.model.getEntities()
+for e in entities:
+    print("entity ", e)
+    types, tags, nodes = gmsh.model.mesh.getElements(e[0], e[1])
     for i in range(len(types)):
         print("type ", types[i])
         print("tags : ", list(tags[i]))
         print("nodes : ", list(nodes[i]))
-    if e[0] == [2] and e[1] == 6 :
-        model.mesh.addElements(e[0],e[1],types,[tags[0][:10]],[nodes[0][:30]])
+    if e[0] == [2] and e[1] == 6:
+        gmsh.model.mesh.addElements(e[0], e[1], types, [tags[0][:10]],
+                                    [nodes[0][:30]])
 
 gmsh.write("mesh_truncated.msh")
 print("Nodes")
-tags, coord, _ = model.mesh.getNodes(2,6)
+tags, coord, _ = gmsh.model.mesh.getNodes(2, 6)
 print(tags)
 print(coord)
 gmsh.finalize()
diff --git a/demos/api/view.py b/demos/api/view.py
index 6bb7a6b9c33eaf36459d96435a964a797fa77e6f..e7c75d34174148ac9c99f6b72d32af7af736da6d 100644
--- a/demos/api/view.py
+++ b/demos/api/view.py
@@ -5,16 +5,11 @@ gmsh.initialize(sys.argv)
 gmsh.option.setNumber("General.Terminal", 1)
 
 # Copied from discrete.py...
-gmsh.model.add("test");
+gmsh.model.add("test")
 gmsh.model.addDiscreteEntity(2, 1)
 gmsh.model.mesh.addNodes(2, 1, [1, 2, 3, 4],
-                         [0., 0., 0.,
-                          1., 0., 0.,
-                          1., 1., 0.,
-                          0., 1., 0.])
-gmsh.model.mesh.addElements(2, 1, [2], [[1, 2]],
-                            [[1, 2, 3,
-                              1, 3, 4]])
+                         [0., 0., 0., 1., 0., 0., 1., 1., 0., 0., 1., 0.])
+gmsh.model.mesh.addElements(2, 1, [2], [[1, 2]], [[1, 2, 3, 1, 3, 4]])
 # ... end of copy
 
 # Create a new post-processing view
@@ -22,9 +17,13 @@ t = gmsh.view.add("some data")
 
 # add 10 steps of model-based data, on the nodes of the mesh
 for step in range(0, 10):
-    gmsh.view.addModelData(t, step, "test", "NodeData",
-                           [1, 2, 3, 4], # tags of nodes
-                           [[10.],[10.],[12.+step],[13.+step]]) # data, per node
+    gmsh.view.addModelData(
+        t,
+        step,
+        "test",
+        "NodeData",
+        [1, 2, 3, 4],  # tags of nodes
+        [[10.], [10.], [12. + step], [13. + step]])  # data, per node
 
 gmsh.view.write(t, "data.msh")
 
diff --git a/demos/api/view_combine.py b/demos/api/view_combine.py
index 2825b606ad43394ae60b18d47d394396fcc8801e..c4b72db2a72ef15b57bb74fcfc1c63d276795a15 100644
--- a/demos/api/view_combine.py
+++ b/demos/api/view_combine.py
@@ -4,16 +4,16 @@ import sys
 gmsh.initialize(sys.argv)
 gmsh.option.setNumber("General.Terminal", 1)
 
-tri1 = [0., 1., 1.,
-        0., 0., 1.,
-        0., 0., 0.];
-tri2 = [0., 1., 0.,
-        0., 1., 1.,
-        0., 0., 0.];
+tri1 = [0., 1., 1., 0., 0., 1., 0., 0., 0.]
+tri2 = [0., 1., 0., 0., 1., 1., 0., 0., 0.]
 
 for step in range(0, 10):
-    tri1.append(10.); tri1.append(10.); tri1.append(12. + step)
-    tri2.append(10.); tri2.append(12. + step); tri2.append(13. + step)
+    tri1.append(10.)
+    tri1.append(10.)
+    tri1.append(12. + step)
+    tri2.append(10.)
+    tri2.append(12. + step)
+    tri2.append(13. + step)
 
 t = gmsh.view.add("some data")
 gmsh.view.addListData(t, "ST", 1, tri1)
diff --git a/demos/api/viewlist.py b/demos/api/viewlist.py
index 88881c722b15c7d827a3a56d9112f100dbd87ba4..6ec6bbc1f065045dd76dea69d829fb0f1e64981d 100644
--- a/demos/api/viewlist.py
+++ b/demos/api/viewlist.py
@@ -4,16 +4,16 @@ import sys
 gmsh.initialize(sys.argv)
 gmsh.option.setNumber("General.Terminal", 1)
 
-tri1 = [0., 1., 1.,
-        0., 0., 1.,
-        0., 0., 0.];
-tri2 = [0., 1., 0.,
-        0., 1., 1.,
-        0., 0., 0.];
+tri1 = [0., 1., 1., 0., 0., 1., 0., 0., 0.]
+tri2 = [0., 1., 0., 0., 1., 1., 0., 0., 0.]
 
 for step in range(0, 10):
-    tri1.append(10.); tri1.append(10.); tri1.append(12. + step)
-    tri2.append(10.); tri2.append(12. + step); tri2.append(13. + step)
+    tri1.append(10.)
+    tri1.append(10.)
+    tri1.append(12. + step)
+    tri2.append(10.)
+    tri2.append(12. + step)
+    tri2.append(13. + step)
 
 t = gmsh.view.add("some data")
 
diff --git a/demos/api/volume.py b/demos/api/volume.py
index d3f7cae9b710da98a583f17e75cca78c37b4a2a2..663e67f1aaabbdf5251ea7d54015fc6a804e55a0 100644
--- a/demos/api/volume.py
+++ b/demos/api/volume.py
@@ -4,7 +4,7 @@ import sys
 gmsh.initialize(sys.argv)
 gmsh.option.setNumber("General.Terminal", 1)
 
-s = gmsh.model.occ.addRectangle(0,0,0, 3,2)
+s = gmsh.model.occ.addRectangle(0, 0, 0, 3, 2)
 gmsh.model.occ.synchronize()
 
 m = gmsh.model.occ.getMass(2, s)
diff --git a/tutorial/python/t1.py b/tutorial/python/t1.py
index e900c14e116b3f3673b4a5a781e7af7e1670b338..ed0d9791a8311e877d54d2d44252eb02048ebdb0 100644
--- a/tutorial/python/t1.py
+++ b/tutorial/python/t1.py
@@ -46,7 +46,7 @@ gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
 #
 # We can then define some additional points. All points should have different
 # tags:
-gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
+gmsh.model.geo.addPoint(.1, 0, 0, lc, 2)
 gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
 
 # If the tag is not provided explicitly, a new tag is automatically created, and
diff --git a/tutorial/python/t10.py b/tutorial/python/t10.py
index e00b3f76c6680976fb62a3225a28fd82d5a07d94..6695ad4dfc3a35d6f9df97367c6bb92ebf9dd443 100644
--- a/tutorial/python/t10.py
+++ b/tutorial/python/t10.py
@@ -20,18 +20,18 @@ gmsh.model.add("t10")
 
 # Let's create a simple rectangular geometry:
 lc = .15
-gmsh.model.geo.addPoint(0.0,0.0,0, lc, 1)
-gmsh.model.geo.addPoint(1,0.0,0, lc, 2)
-gmsh.model.geo.addPoint(1,1,0, lc, 3)
-gmsh.model.geo.addPoint(0,1,0, lc, 4)
-gmsh.model.geo.addPoint(0.2,.5,0, lc, 5)
-
-gmsh.model.geo.addLine(1,2, 1);
-gmsh.model.geo.addLine(2,3, 2);
-gmsh.model.geo.addLine(3,4, 3);
-gmsh.model.geo.addLine(4,1, 4);
-
-gmsh.model.geo.addCurveLoop([1,2,3,4], 5)
+gmsh.model.geo.addPoint(0.0, 0.0, 0, lc, 1)
+gmsh.model.geo.addPoint(1, 0.0, 0, lc, 2)
+gmsh.model.geo.addPoint(1, 1, 0, lc, 3)
+gmsh.model.geo.addPoint(0, 1, 0, lc, 4)
+gmsh.model.geo.addPoint(0.2, .5, 0, lc, 5)
+
+gmsh.model.geo.addLine(1, 2, 1)
+gmsh.model.geo.addLine(2, 3, 2)
+gmsh.model.geo.addLine(3, 4, 3)
+gmsh.model.geo.addLine(4, 1, 4)
+
+gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 5)
 gmsh.model.geo.addPlaneSurface([5], 6)
 
 gmsh.model.geo.synchronize()
@@ -57,8 +57,8 @@ gmsh.model.mesh.field.setNumbers(1, "EdgesList", [2])
 # LcMin -o----------------/
 #        |                |       |
 #      Point           DistMin DistMax
-gmsh.model.mesh.field.add("Threshold", 2);
-gmsh.model.mesh.field.setNumber(2, "IField", 1);
+gmsh.model.mesh.field.add("Threshold", 2)
+gmsh.model.mesh.field.setNumber(2, "IField", 1)
 gmsh.model.mesh.field.setNumber(2, "LcMin", lc / 30)
 gmsh.model.mesh.field.setNumber(2, "LcMax", lc)
 gmsh.model.mesh.field.setNumber(2, "DistMin", 0.15)
@@ -67,7 +67,8 @@ gmsh.model.mesh.field.setNumber(2, "DistMax", 0.5)
 # Say we want to modulate the mesh element sizes using a mathematical function
 # of the spatial coordinates. We can do this with the MathEval field:
 gmsh.model.mesh.field.add("MathEval", 3)
-gmsh.model.mesh.field.setString(3, "F", "Cos(4*3.14*x) * Sin(4*3.14*y) / 10 + 0.101")
+gmsh.model.mesh.field.setString(3, "F",
+                                "Cos(4*3.14*x) * Sin(4*3.14*y) / 10 + 0.101")
 
 # We could also combine MathEval with values coming from other fields. For
 # example, let's define a `Distance' field around point 1
@@ -77,7 +78,7 @@ gmsh.model.mesh.field.setNumbers(4, "NodesList", [1])
 # We can then create a `MathEval' field with a function that depends on the
 # return value of the `Distance' field 4, i.e., depending on the distance to
 # point 1 (here using a cubic law, with minimum element size = lc / 100)
-gmsh.model.mesh.field.add("MathEval", 5);
+gmsh.model.mesh.field.add("MathEval", 5)
 gmsh.model.mesh.field.setString(5, "F", "F4^3 + " + str(lc / 100))
 
 # We could also use a `Box' field to impose a step change in element sizes
diff --git a/tutorial/python/t13.py b/tutorial/python/t13.py
index b0de023edd7eb92e1c568eec45dfd601c87e6d8d..e24fdd4c9165690bf5261e71b0bf4747fd31e92e 100644
--- a/tutorial/python/t13.py
+++ b/tutorial/python/t13.py
@@ -22,7 +22,7 @@ gmsh.merge(os.path.join(path, '..', 't13_data.stl'))
 # We first classify ("color") the surfaces by splitting the original surface
 # along sharp geometrical features. This will create new discrete surfaces,
 # curves and points.
-angle = 40 # Angle for surface detection
+angle = 40  # Angle for surface detection
 
 # For complex geometries, patches can be too complex, too elongated or too large
 # to be parametrized; setting the following option will force the creation of
@@ -33,12 +33,11 @@ forceParametrizablePatches = False
 includeBoundary = True
 
 # Force curves to be split on given angle:
-curveAngle = 180;
+curveAngle = 180
 
-gmsh.model.mesh.classifySurfaces(angle * math.pi/180.,
-                                 includeBoundary,
+gmsh.model.mesh.classifySurfaces(angle * math.pi / 180., includeBoundary,
                                  forceParametrizablePatches,
-                                 curveAngle * math.pi/180.)
+                                 curveAngle * math.pi / 180.)
 
 # Create a geometry for all the discrete curves and surfaces in the mesh, by
 # computing a parametrization for each one
@@ -53,11 +52,11 @@ gmsh.model.geo.synchronize()
 
 # We specify element sizes imposed by a size field, just because we can :-)
 funny = False
-f = gmsh.model.mesh.field.add("MathEval");
+f = gmsh.model.mesh.field.add("MathEval")
 if funny:
-  gmsh.model.mesh.field.setString(f, "F", "2*Sin((x+y)/5) + 3")
+    gmsh.model.mesh.field.setString(f, "F", "2*Sin((x+y)/5) + 3")
 else:
-  gmsh.model.mesh.field.setString(f, "F", "4")
+    gmsh.model.mesh.field.setString(f, "F", "4")
 gmsh.model.mesh.field.setAsBackgroundMesh(f)
 
 gmsh.model.mesh.generate(3)
diff --git a/tutorial/python/t14.py b/tutorial/python/t14.py
index 632df29c5df390afb4836f25f1a07f1c8b9d9f85..422da695f41d64b1e0c0192b53cfdcc16a97ab13 100644
--- a/tutorial/python/t14.py
+++ b/tutorial/python/t14.py
@@ -20,8 +20,8 @@ gmsh.option.setNumber("General.Terminal", 1)
 # Create an example geometry
 gmsh.model.add("t14")
 
-m = 0.5; # mesh characteristic length
-h = 2;   # geometry height in the z-direction
+m = 0.5  # mesh characteristic length
+h = 2  # geometry height in the z-direction
 
 gmsh.model.geo.addPoint(0, 0, 0, m, 1)
 gmsh.model.geo.addPoint(10, 0, 0, m, 2)
@@ -58,26 +58,29 @@ gmsh.model.geo.addCurveLoop([6, 7, 8, 1, 2, 3, 4, 5], 13)
 gmsh.model.geo.addCurveLoop([11, 12, 9, 10], 14)
 gmsh.model.geo.addPlaneSurface([13, 14], 15)
 
-ext_tags = gmsh.model.geo.extrude([(2,15)], 0, 0, h)
+ext_tags = gmsh.model.geo.extrude([(2, 15)], 0, 0, h)
 
 # Create physical groups, which are used to define the domain of the
 # (co)homology computation and the subdomain of the relative (co)homology
 # computation.
 
 # Whole domain
-domain_tag = 1;
-domain_physical_tag = 1001;
+domain_tag = 1
+domain_physical_tag = 1001
 gmsh.model.addPhysicalGroup(dim=3, tags=[domain_tag], tag=domain_physical_tag)
 gmsh.model.setPhysicalName(dim=3, tag=domain_physical_tag, name="Whole domain")
 
 # Four "terminals" of the model
-terminal_tags = [36, 44, 52, 60];
+terminal_tags = [36, 44, 52, 60]
 terminals_physical_tag = 2001
-gmsh.model.addPhysicalGroup(dim=2, tags=terminal_tags, tag=terminals_physical_tag)
+gmsh.model.addPhysicalGroup(dim=2,
+                            tags=terminal_tags,
+                            tag=terminals_physical_tag)
 gmsh.model.setPhysicalName(dim=2, tag=terminals_physical_tag, name="Terminals")
 
 # Find domain boundary tags
-boundary_dimtags = gmsh.model.getBoundary(dimTags=[(3, domain_tag)], oriented=False)
+boundary_dimtags = gmsh.model.getBoundary(dimTags=[(3, domain_tag)],
+                                          oriented=False)
 boundary_tags = []
 complement_tags = []
 for tag in boundary_dimtags:
@@ -88,16 +91,18 @@ for tag in terminal_tags:
 
 # Whole domain surface
 boundary_physical_tag = 2002
-gmsh.model.addPhysicalGroup(dim=2, tags=boundary_tags,
+gmsh.model.addPhysicalGroup(dim=2,
+                            tags=boundary_tags,
                             tag=boundary_physical_tag)
-gmsh.model.setPhysicalName(dim=2, tag=boundary_physical_tag,
-                           name="Boundary")
+gmsh.model.setPhysicalName(dim=2, tag=boundary_physical_tag, name="Boundary")
 
 # Complement of the domain surface respect to the four terminals
 complement_physical_tag = 2003
-gmsh.model.addPhysicalGroup(dim=2, tags=complement_tags,
+gmsh.model.addPhysicalGroup(dim=2,
+                            tags=complement_tags,
                             tag=complement_physical_tag)
-gmsh.model.setPhysicalName(dim=2, tag=complement_physical_tag,
+gmsh.model.setPhysicalName(dim=2,
+                           tag=complement_physical_tag,
                            name="Complement")
 
 gmsh.model.geo.synchronize()
@@ -106,19 +111,19 @@ gmsh.model.geo.synchronize()
 # terminals.
 gmsh.model.mesh.computeHomology(domainTags=[domain_physical_tag],
                                 subdomainTags=[terminals_physical_tag],
-                                dims=[0,1,2,3])
+                                dims=[0, 1, 2, 3])
 
 # Find homology space bases isomorphic to the previous bases: homology spaces
 # modulo the non-terminal domain surface, a.k.a the thin cuts.
 gmsh.model.mesh.computeHomology(domainTags=[domain_physical_tag],
                                 subdomainTags=[complement_physical_tag],
-                                dims=[0,1,2,3])
+                                dims=[0, 1, 2, 3])
 
 # Find cohomology space bases isomorphic to the previous bases: cohomology
 # spaces of the domain modulo the four terminals, a.k.a the thick cuts.
 gmsh.model.mesh.computeCohomology(domainTags=[domain_physical_tag],
                                   subdomainTags=[terminals_physical_tag],
-                                  dims=[0,1,2,3])
+                                  dims=[0, 1, 2, 3])
 
 # more examples
 # gmsh.model.mesh.computeHomology()
diff --git a/tutorial/python/t15.py b/tutorial/python/t15.py
index 76648972fe1f4658fd086aae4afbd3ffc2fdec9e..27099078b23be900b893f0bf6a124205cdecfef5 100644
--- a/tutorial/python/t15.py
+++ b/tutorial/python/t15.py
@@ -22,7 +22,7 @@ gmsh.option.setNumber("General.Terminal", 1)
 # Copied from t1.py...
 lc = 1e-2
 gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
-gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
+gmsh.model.geo.addPoint(.1, 0, 0, lc, 2)
 gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
 gmsh.model.geo.addPoint(0, .3, 0, lc, 4)
 gmsh.model.geo.addLine(1, 2, 1)
@@ -63,24 +63,24 @@ p = gmsh.model.geo.addPoint(0.07, 0.15, 0.025, lc)
 gmsh.model.geo.synchronize()
 gmsh.model.mesh.embed(0, [p], 3, 1)
 
-gmsh.model.geo.addPoint(0.025, 0.15, 0.025, lc, p+1)
-l = gmsh.model.geo.addLine(7, p+1)
+gmsh.model.geo.addPoint(0.025, 0.15, 0.025, lc, p + 1)
+l = gmsh.model.geo.addLine(7, p + 1)
 
 gmsh.model.geo.synchronize()
 gmsh.model.mesh.embed(1, [l], 3, 1)
 
 # Finally, we can also embed a surface in a volume:
-gmsh.model.geo.addPoint(0.02, 0.12, 0.05, lc, p+2)
-gmsh.model.geo.addPoint(0.04, 0.12, 0.05, lc, p+3)
-gmsh.model.geo.addPoint(0.04, 0.18, 0.05, lc, p+4)
-gmsh.model.geo.addPoint(0.02, 0.18, 0.05, lc, p+5)
+gmsh.model.geo.addPoint(0.02, 0.12, 0.05, lc, p + 2)
+gmsh.model.geo.addPoint(0.04, 0.12, 0.05, lc, p + 3)
+gmsh.model.geo.addPoint(0.04, 0.18, 0.05, lc, p + 4)
+gmsh.model.geo.addPoint(0.02, 0.18, 0.05, lc, p + 5)
 
-gmsh.model.geo.addLine(p+2, p+3, l+1)
-gmsh.model.geo.addLine(p+3, p+4, l+2)
-gmsh.model.geo.addLine(p+4, p+5, l+3)
-gmsh.model.geo.addLine(p+5, p+2, l+4)
+gmsh.model.geo.addLine(p + 2, p + 3, l + 1)
+gmsh.model.geo.addLine(p + 3, p + 4, l + 2)
+gmsh.model.geo.addLine(p + 4, p + 5, l + 3)
+gmsh.model.geo.addLine(p + 5, p + 2, l + 4)
 
-ll = gmsh.model.geo.addCurveLoop([l+1, l+2, l+3, l+4])
+ll = gmsh.model.geo.addCurveLoop([l + 1, l + 2, l + 3, l + 4])
 s = gmsh.model.geo.addPlaneSurface([ll])
 
 gmsh.model.geo.synchronize()
diff --git a/tutorial/python/t16.py b/tutorial/python/t16.py
index f551cb2df8af5ba5058359682a1ae742c89245d2..921239b1fb418e655bb04f9a4ffa4a3acae22195 100644
--- a/tutorial/python/t16.py
+++ b/tutorial/python/t16.py
@@ -25,30 +25,33 @@ gmsh.model.add("t16")
 gmsh.logger.start()
 
 # We first create two cubes:
-gmsh.model.occ.addBox(0,0,0, 1,1,1, 1)
-gmsh.model.occ.addBox(0,0,0, 0.5,0.5,0.5, 2)
+gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1, 1)
+gmsh.model.occ.addBox(0, 0, 0, 0.5, 0.5, 0.5, 2)
 
 # We apply a boolean difference to create the "cube minus one eigth" shape:
-gmsh.model.occ.cut([(3,1)], [(3,2)], 3)
+gmsh.model.occ.cut([(3, 1)], [(3, 2)], 3)
 
 # Boolean operations with OpenCASCADE always create new entities. By default the
 # extra arguments `removeObject' and `removeTool' in `cut()' are set to `True',
 # which will delete the original entities.
 
 # We then create the five spheres:
-x = 0; y = 0.75; z = 0; r = 0.09
+x = 0
+y = 0.75
+z = 0
+r = 0.09
 holes = []
 for t in range(1, 6):
     x += 0.166
     z += 0.166
-    gmsh.model.occ.addSphere(x,y,z,r, 3 + t)
+    gmsh.model.occ.addSphere(x, y, z, r, 3 + t)
     holes.append((3, 3 + t))
 
 # If we had wanted five empty holes we would have used `cut()' again. Here we
 # want five spherical inclusions, whose mesh should be conformal with the mesh
 # of the cube: we thus use `fragment()', which intersects all volumes in a
 # conformal manner (without creating duplicate interfaces):
-ov, ovv = gmsh.model.occ.fragment([(3,3)], holes)
+ov, ovv = gmsh.model.occ.fragment([(3, 3)], holes)
 
 # ov contains all the generated entities of the same dimension as the input
 # entities:
@@ -58,7 +61,7 @@ for e in ov:
 
 # ovv contains the parent-child relationships for all the input entities:
 print("before/after fragment relations:")
-for e in zip([(3,3)] + holes, ovv):
+for e in zip([(3, 3)] + holes, ovv):
     print("parent " + str(e[0]) + " -> child " + str(e[1]))
 
 gmsh.model.occ.synchronize()
@@ -94,12 +97,13 @@ lcar3 = .055
 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lcar1)
 
 # Override this constraint on the points of the five spheres:
-gmsh.model.mesh.setSize(gmsh.model.getBoundary(holes, False, False, True), lcar3)
+gmsh.model.mesh.setSize(gmsh.model.getBoundary(holes, False, False, True),
+                        lcar3)
 
 # Select the corner point by searching for it geometrically:
 eps = 1e-3
-ov = gmsh.model.getEntitiesInBoundingBox(0.5-eps, 0.5-eps, 0.5-eps,
-                                         0.5+eps, 0.5+eps, 0.5+eps, 0)
+ov = gmsh.model.getEntitiesInBoundingBox(0.5 - eps, 0.5 - eps, 0.5 - eps,
+                                         0.5 + eps, 0.5 + eps, 0.5 + eps, 0)
 gmsh.model.mesh.setSize(ov, lcar2)
 
 gmsh.model.mesh.generate(3)
diff --git a/tutorial/python/t18.py b/tutorial/python/t18.py
index a0d16e8cb4fb07c012176457fc21090776e10172..42c3e67a1f07228c7e74b4f6603a6b96da39e868 100644
--- a/tutorial/python/t18.py
+++ b/tutorial/python/t18.py
@@ -27,12 +27,12 @@ gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1, 1)
 gmsh.model.occ.synchronize()
 
 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), 0.1)
-gmsh.model.mesh.setSize([(0,1)], 0.02)
+gmsh.model.mesh.setSize([(0, 1)], 0.02)
 
 # To impose that the mesh on surface 2 (the right side of the cube) should
 # match the mesh from surface 1 (the left side), the following periodicity
 # constraint is set:
-translation = [1,0,0,1, 0,1,0,0, 0,0,1,0, 0,0,0,1]
+translation = [1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1]
 gmsh.model.mesh.setPeriodic(2, [2], [1], translation)
 
 # The periodicity transform is provided as a 4x4 affine transformation matrix,
@@ -42,8 +42,10 @@ gmsh.model.mesh.setPeriodic(2, [2], [1], translation)
 # the mesh from surface 1.
 
 # Multiple periodicities can be imposed in the same way:
-gmsh.model.mesh.setPeriodic(2, [6], [5], [1,0,0,0, 0,1,0,0, 0,0,1,1, 0,0,0,1])
-gmsh.model.mesh.setPeriodic(2, [4], [3], [1,0,0,0, 0,1,0,1, 0,0,1,0, 0,0,0,1])
+gmsh.model.mesh.setPeriodic(2, [6], [5],
+                            [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1])
+gmsh.model.mesh.setPeriodic(2, [4], [3],
+                            [1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1])
 
 # For more complicated cases, finding the corresponding surfaces by hand can
 # be tedious, especially when geometries are created through solid
@@ -51,19 +53,21 @@ gmsh.model.mesh.setPeriodic(2, [4], [3], [1,0,0,0, 0,1,0,1, 0,0,1,0, 0,0,0,1])
 
 # We start with a cube and some spheres:
 gmsh.model.occ.addBox(2, 0, 0, 1, 1, 1, 10)
-x = 2-0.3; y = 0; z = 0
+x = 2 - 0.3
+y = 0
+z = 0
 gmsh.model.occ.addSphere(x, y, z, 0.35, 11)
-gmsh.model.occ.addSphere(x+1, y, z, 0.35, 12)
-gmsh.model.occ.addSphere(x, y+1, z, 0.35, 13)
-gmsh.model.occ.addSphere(x, y, z+1, 0.35, 14)
-gmsh.model.occ.addSphere(x+1, y+1, z, 0.35, 15)
-gmsh.model.occ.addSphere(x, y+1, z+1, 0.35, 16)
-gmsh.model.occ.addSphere(x+1, y, z+1, 0.35, 17)
-gmsh.model.occ.addSphere(x+1, y+1, z+1, 0.35, 18)
+gmsh.model.occ.addSphere(x + 1, y, z, 0.35, 12)
+gmsh.model.occ.addSphere(x, y + 1, z, 0.35, 13)
+gmsh.model.occ.addSphere(x, y, z + 1, 0.35, 14)
+gmsh.model.occ.addSphere(x + 1, y + 1, z, 0.35, 15)
+gmsh.model.occ.addSphere(x, y + 1, z + 1, 0.35, 16)
+gmsh.model.occ.addSphere(x + 1, y, z + 1, 0.35, 17)
+gmsh.model.occ.addSphere(x + 1, y + 1, z + 1, 0.35, 18)
 
 # We first fragment all the volumes, which will leave parts of spheres
 # protruding outside the cube:
-out, _ = gmsh.model.occ.fragment([(3,10)], [(3,i) for i in range(11,19)])
+out, _ = gmsh.model.occ.fragment([(3, 10)], [(3, i) for i in range(11, 19)])
 gmsh.model.occ.synchronize()
 
 # Ask OpenCASCADE to compute more accurate bounding boxes of entities using
@@ -73,18 +77,17 @@ gmsh.option.setNumber("Geometry.OCCBoundsUseStl", 1)
 # We then retrieve all the volumes in the bounding box of the original cube,
 # and delete all the parts outside it:
 eps = 1e-3
-vin = gmsh.model.getEntitiesInBoundingBox(2-eps,-eps,-eps,
-                                          2+1+eps,1+eps,1+eps,
-                                          3)
-for v in vin: out.remove(v)
-gmsh.model.removeEntities(out, True) # Delete outside parts recursively
+vin = gmsh.model.getEntitiesInBoundingBox(2 - eps, -eps, -eps, 2 + 1 + eps,
+                                          1 + eps, 1 + eps, 3)
+for v in vin:
+    out.remove(v)
+gmsh.model.removeEntities(out, True)  # Delete outside parts recursively
 
 # We now set a non-uniform mesh size constraint (again to check results
 # visually):
-p = gmsh.model.getBoundary(vin, False, False, True) # Get all points
+p = gmsh.model.getBoundary(vin, False, False, True)  # Get all points
 gmsh.model.mesh.setSize(p, 0.1)
-p = gmsh.model.getEntitiesInBoundingBox(2-eps, -eps, -eps,
-                                        2+eps, eps, eps,
+p = gmsh.model.getEntitiesInBoundingBox(2 - eps, -eps, -eps, 2 + eps, eps, eps,
                                         0)
 gmsh.model.mesh.setSize(p, 0.001)
 
@@ -92,18 +95,17 @@ gmsh.model.mesh.setSize(p, 0.001)
 # geometry automatically.
 
 # First we get all surfaces on the left:
-sxmin = gmsh.model.getEntitiesInBoundingBox(2-eps, -eps, -eps,
-                                            2+eps, 1+eps, 1+eps,
-                                            2)
+sxmin = gmsh.model.getEntitiesInBoundingBox(2 - eps, -eps, -eps, 2 + eps,
+                                            1 + eps, 1 + eps, 2)
 
 for i in sxmin:
     # Then we get the bounding box of each left surface
     xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(i[0], i[1])
     # We translate the bounding box to the right and look for surfaces inside
     # it:
-    sxmax = gmsh.model.getEntitiesInBoundingBox(xmin-eps+1, ymin-eps, zmin-eps,
-                                                xmax+eps+1, ymax+eps, zmax+eps,
-                                                2)
+    sxmax = gmsh.model.getEntitiesInBoundingBox(xmin - eps + 1, ymin - eps,
+                                                zmin - eps, xmax + eps + 1,
+                                                ymax + eps, zmax + eps, 2)
     # For all the matches, we compare the corresponding bounding boxes...
     for j in sxmax:
         xmin2, ymin2, zmin2, xmax2, ymax2, zmax2 = gmsh.model.getBoundingBox(
@@ -111,9 +113,9 @@ for i in sxmin:
         xmin2 -= 1
         xmax2 -= 1
         # ...and if they match, we apply the periodicity constraint
-        if(abs(xmin2 - xmin) < eps and abs(xmax2 - xmax) < eps and
-           abs(ymin2 - ymin) < eps and abs(ymax2 - ymax) < eps and
-           abs(zmin2 - zmin) < eps and abs(zmax2 - zmax) < eps):
+        if (abs(xmin2 - xmin) < eps and abs(xmax2 - xmax) < eps
+                and abs(ymin2 - ymin) < eps and abs(ymax2 - ymax) < eps
+                and abs(zmin2 - zmin) < eps and abs(zmax2 - zmax) < eps):
             gmsh.model.mesh.setPeriodic(2, [j[1]], [i[1]], translation)
 
 gmsh.model.mesh.generate(3)
diff --git a/tutorial/python/t19.py b/tutorial/python/t19.py
index a51db6dcbe6203d96ee71432f3f302773d09a3b6..4dee47800d17c4a59ecc841375e7c439e9bac7e4 100644
--- a/tutorial/python/t19.py
+++ b/tutorial/python/t19.py
@@ -20,23 +20,23 @@ gmsh.model.add("t19")
 
 # Volumes can be constructed from (closed) curve loops thanks to the
 # `addThruSections()' function
-gmsh.model.occ.addCircle(0,0,0, 0.5, 1)
+gmsh.model.occ.addCircle(0, 0, 0, 0.5, 1)
 gmsh.model.occ.addCurveLoop([1], 1)
-gmsh.model.occ.addCircle(0.1,0.05,1, 0.1, 2)
+gmsh.model.occ.addCircle(0.1, 0.05, 1, 0.1, 2)
 gmsh.model.occ.addCurveLoop([2], 2)
-gmsh.model.occ.addCircle(-0.1,-0.1,2, 0.3, 3)
+gmsh.model.occ.addCircle(-0.1, -0.1, 2, 0.3, 3)
 gmsh.model.occ.addCurveLoop([3], 3)
-gmsh.model.occ.addThruSections([1,2,3], 1)
+gmsh.model.occ.addThruSections([1, 2, 3], 1)
 gmsh.model.occ.synchronize()
 
 # We can also force the creation of ruled surfaces:
-gmsh.model.occ.addCircle(2+0,0,0, 0.5, 11)
+gmsh.model.occ.addCircle(2 + 0, 0, 0, 0.5, 11)
 gmsh.model.occ.addCurveLoop([11], 11)
-gmsh.model.occ.addCircle(2+0.1,0.05,1, 0.1, 12)
+gmsh.model.occ.addCircle(2 + 0.1, 0.05, 1, 0.1, 12)
 gmsh.model.occ.addCurveLoop([12], 12)
-gmsh.model.occ.addCircle(2-0.1,-0.1,2, 0.3, 13)
+gmsh.model.occ.addCircle(2 - 0.1, -0.1, 2, 0.3, 13)
 gmsh.model.occ.addCurveLoop([13], 13)
-gmsh.model.occ.addThruSections([11,12,13], 11, True, True)
+gmsh.model.occ.addThruSections([11, 12, 13], 11, True, True)
 gmsh.model.occ.synchronize()
 
 # We copy the first volume, and fillet all its edges:
@@ -55,9 +55,9 @@ r = 1.
 h = 1. * nturns
 p = []
 for i in range(0, npts):
-    theta = i * 2 * math.pi * nturns/npts
+    theta = i * 2 * math.pi * nturns / npts
     gmsh.model.occ.addPoint(r * math.cos(theta), r * math.sin(theta),
-                            i * h/npts, 1, 1000 + i)
+                            i * h / npts, 1, 1000 + i)
     p.append(1000 + i)
 gmsh.model.occ.addSpline(p, 1000)
 
@@ -65,11 +65,11 @@ gmsh.model.occ.addSpline(p, 1000)
 gmsh.model.occ.addWire([1000], 1000)
 
 # We define the shape we would like to extrude along the spline (a disk):
-gmsh.model.occ.addDisk(1,0,0, 0.2, 0.2, 1000)
-gmsh.model.occ.rotate([(2,1000)], 0, 0, 0, 1, 0, 0, math.pi/2)
+gmsh.model.occ.addDisk(1, 0, 0, 0.2, 0.2, 1000)
+gmsh.model.occ.rotate([(2, 1000)], 0, 0, 0, 1, 0, 0, math.pi / 2)
 
 # We extrude the disk along the spline to create a pipe:
-gmsh.model.occ.addPipe([(2,1000)], 1000)
+gmsh.model.occ.addPipe([(2, 1000)], 1000)
 
 # We delete the source surface, and increase the number of sub-edges for a
 # nicer display of the geometry:
diff --git a/tutorial/python/t2.py b/tutorial/python/t2.py
index b295dcfbed6e2ad8887be2fc773095d2829c6ae0..0c364592ef812e1557e2b807d7814578a3549a53 100644
--- a/tutorial/python/t2.py
+++ b/tutorial/python/t2.py
@@ -21,7 +21,7 @@ gmsh.model.add("t2")
 # Copied from t1.py...
 lc = 1e-2
 gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
-gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
+gmsh.model.geo.addPoint(.1, 0, 0, lc, 2)
 gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
 gmsh.model.geo.addPoint(0, .3, 0, lc, 4)
 gmsh.model.geo.addLine(1, 2, 1)
@@ -48,7 +48,7 @@ gmsh.model.geo.translate([(0, 5)], -0.02, 0, 0)
 
 # And it can be further rotated by -Pi/4 around (0, 0.3, 0) (with the rotation
 # along the z axis) with:
-gmsh.model.geo.rotate([(0, 5)], 0,0.3,0, 0,0,1, -math.pi/4)
+gmsh.model.geo.rotate([(0, 5)], 0, 0.3, 0, 0, 0, 1, -math.pi / 4)
 
 # Note that there are no units in Gmsh: coordinates are just numbers - it's
 # up to the user to associate a meaning to them.
@@ -63,7 +63,7 @@ gmsh.model.geo.translate(ov, 0, 0.05, 0)
 # lines:
 gmsh.model.geo.addLine(3, ov[0][1], 7)
 gmsh.model.geo.addLine(ov[0][1], 5, 8)
-gmsh.model.geo.addCurveLoop([5,-8,-7,3], 10)
+gmsh.model.geo.addCurveLoop([5, -8, -7, 3], 10)
 gmsh.model.geo.addPlaneSurface([10], 11)
 
 # In the same way, we can translate copies of the two surfaces 1 and 11 to the
@@ -121,12 +121,12 @@ ov2 = gmsh.model.geo.extrude([ov[1]], 0, 0, 0.12)
 
 # Mesh sizes associated to geometrical points can be set by passing a vector of
 # (dim, tag) pairs for the corresponding points:
-gmsh.model.geo.mesh.setSize([(0,103), (0,105), (0,109), (0,102), (0,28),
-                             (0, 24), (0,6), (0,5)], lc * 3)
+gmsh.model.geo.mesh.setSize([(0, 103), (0, 105), (0, 109), (0, 102), (0, 28),
+                             (0, 24), (0, 6), (0, 5)], lc * 3)
 
 # We finally group volumes 129 and 130 in a single physical group with tag `1'
 # and name "The volume":
-gmsh.model.addPhysicalGroup(3, [129,130], 1)
+gmsh.model.addPhysicalGroup(3, [129, 130], 1)
 gmsh.model.setPhysicalName(3, 1, "The volume")
 
 # We finish by synchronizing the data from the built-in geometry kernel with the
diff --git a/tutorial/python/t20.py b/tutorial/python/t20.py
index dc34c0996d457b57af29d516c4f61d40b06b35f5..7dd7c584fb56a87ac456940129b888f274a5d613 100644
--- a/tutorial/python/t20.py
+++ b/tutorial/python/t20.py
@@ -33,7 +33,8 @@ v = gmsh.model.occ.importShapes(os.path.join(path, '..', 't20_data.step'))
 
 # Get the bounding box of the volume:
 gmsh.model.occ.synchronize()
-xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(v[0][0], v[0][1])
+xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(
+    v[0][0], v[0][1])
 
 # Note that the synchronization step can be avoided in Gmsh >= 4.6 by using
 # `gmsh.model.occ.getBoundingBox()' instead of `gmsh.model.getBoundingBox()'.
@@ -41,19 +42,18 @@ xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(v[0][0], v[0][1])
 # We want to slice the model into N slices, and either keep the volume slices
 # or just the surfaces obtained by the cutting:
 
-N = 5 # Number of slices
-surf = False # Keep only surfaces?
+N = 5  # Number of slices
+surf = False  # Keep only surfaces?
 
-dz = (zmax - zmin) / N;
+dz = (zmax - zmin) / N
 
 # Define the cutting planes
 for i in range(1, N):
-    gmsh.model.occ.addRectangle(xmin, ymin, zmin + i * dz,
-                                xmax-xmin, ymax-ymin,
-                                1000 + i)
+    gmsh.model.occ.addRectangle(xmin, ymin, zmin + i * dz, xmax - xmin,
+                                ymax - ymin, 1000 + i)
 
 # Fragment (i.e. intersect) the volume with all the cutting planes:
-gmsh.model.occ.fragment(v, [(2, i) for i in range(1000+1, 1000+N)])
+gmsh.model.occ.fragment(v, [(2, i) for i in range(1000 + 1, 1000 + N)])
 
 # Now remove all the surfaces (and their bounding entities) that are not on the
 # boundary of a volume, i.e. the parts of the cutting planes that "stick out" of
@@ -73,9 +73,11 @@ if surf:
     eps = 1e-4
     s = []
     for i in range(1, N):
-        s.extend(gmsh.model.getEntitiesInBoundingBox(
-            xmin-eps,ymin-eps,zmin + i * dz - eps,
-            xmax+eps,ymax+eps,zmin + i * dz + eps, 2))
+        s.extend(
+            gmsh.model.getEntitiesInBoundingBox(xmin - eps, ymin - eps,
+                                                zmin + i * dz - eps,
+                                                xmax + eps, ymax + eps,
+                                                zmin + i * dz + eps, 2))
     # ...and remove all the other entities (here directly in the model, as we
     # won't modify any OpenCASCADE entities later on):
     dels = gmsh.model.getEntities(2)
diff --git a/tutorial/python/t21.py b/tutorial/python/t21.py
index c29c9519c3886f0de803825e7e314caa0a7d8379..b7012e3dc8e459e9b3fe031c2c0dcef2141fbef6 100644
--- a/tutorial/python/t21.py
+++ b/tutorial/python/t21.py
@@ -31,7 +31,7 @@ gmsh.option.setNumber("General.Terminal", 1)
 gmsh.model.add("t21")
 gmsh.model.occ.addRectangle(0, 0, 0, 1, 1, 1)
 gmsh.model.occ.addRectangle(1, 0, 0, 1, 1, 2)
-gmsh.model.occ.fragment([(2,1)], [(2,2)])
+gmsh.model.occ.fragment([(2, 1)], [(2, 2)])
 gmsh.model.occ.synchronize()
 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), 0.05)
 
@@ -44,13 +44,13 @@ gmsh.model.setPhysicalName(2, 200, "Right")
 gmsh.model.mesh.generate(2)
 
 # We now define several constants to fine-tune how the mesh will be partitioned
-partitioner = 0 # 0 for Metis, 1 for SimplePartition
-N = 3 # Number of partitions
-topology = 1 # Create partition topology (BRep)?
-ghosts = 0 # Create ghost cells?
-physicals = 0 # Create new physical groups?
-write = 1 # Write file to disk?
-split = 1 # Write one file per partition?
+partitioner = 0  # 0 for Metis, 1 for SimplePartition
+N = 3  # Number of partitions
+topology = 1  # Create partition topology (BRep)?
+ghosts = 0  # Create ghost cells?
+physicals = 0  # Create new physical groups?
+write = 1  # Write file to disk?
+split = 1  # Write one file per partition?
 
 # Should we create the boundary representation of the partition entities?
 gmsh.option.setNumber("Mesh.PartitionCreateTopology", topology)
@@ -92,7 +92,8 @@ entities = gmsh.model.getEntities()
 for e in entities:
     partitions = gmsh.model.getPartitions(e[0], e[1])
     if len(partitions):
-        print("Entity " + str(e) + " of type " + gmsh.model.getType(e[0], e[1]))
+        print("Entity " + str(e) + " of type " +
+              gmsh.model.getType(e[0], e[1]))
         print(" - Partition(s): " + str(partitions))
         print(" - Parent: " + str(gmsh.model.getParent(e[0], e[1])))
         print(" - Boundary: " + str(gmsh.model.getBoundary([e])))
diff --git a/tutorial/python/t3.py b/tutorial/python/t3.py
index 88960d2ea46044e0ca419a790afd39ef2d2262fb..a98e33851ea059908c6d4ac136ab1cae0c309a8f 100644
--- a/tutorial/python/t3.py
+++ b/tutorial/python/t3.py
@@ -17,7 +17,7 @@ gmsh.model.add("t3")
 # Copied from t1.py...
 lc = 1e-2
 gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
-gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
+gmsh.model.geo.addPoint(.1, 0, 0, lc, 2)
 gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
 gmsh.model.geo.addPoint(0, .3, 0, lc, 4)
 gmsh.model.geo.addLine(1, 2, 1)
@@ -40,13 +40,14 @@ gmsh.model.setPhysicalName(2, ps, "My surface")
 
 h = 0.1
 angle = 90.
-ov = gmsh.model.geo.extrude([(2,1)], 0, 0, h, [8,2], [0.5,1])
+ov = gmsh.model.geo.extrude([(2, 1)], 0, 0, h, [8, 2], [0.5, 1])
 
 # The extrusion can also be performed with a rotation instead of a translation,
 # and the resulting mesh can be recombined into prisms (we use only one layer
 # here, with 7 subdivisions). All rotations are specified by an an axis point
 # (-0.1, 0, 0.1), an axis direction (0, 1, 0), and a rotation angle (-Pi/2):
-ov = gmsh.model.geo.revolve([(2,28)], -0.1,0,0.1, 0,1,0, -math.pi/2, [7])
+ov = gmsh.model.geo.revolve([(2, 28)], -0.1, 0, 0.1, 0, 1, 0, -math.pi / 2,
+                            [7])
 
 # Using the built-in geometry kernel, only rotations with angles < Pi are
 # supported. To do a full turn, you will thus need to apply at least 3
@@ -56,8 +57,8 @@ ov = gmsh.model.geo.revolve([(2,28)], -0.1,0,0.1, 0,1,0, -math.pi/2, [7])
 # also be combined to form a "twist".  The last (optional) argument for the
 # extrude() and twist() functions specifies whether the extruded mesh should be
 # recombined or not.
-ov = gmsh.model.geo.twist([(2,50)], 0,0.15,0.25, -2*h,0,0, 1,0,0,
-                          angle*math.pi/180., [10], [], True)
+ov = gmsh.model.geo.twist([(2, 50)], 0, 0.15, 0.25, -2 * h, 0, 0, 1, 0, 0,
+                          angle * math.pi / 180., [10], [], True)
 
 gmsh.model.geo.synchronize()
 
@@ -76,16 +77,16 @@ gmsh.write("t3.msh")
 # through the API, we can for example make point tags visible or redefine some
 # colors:
 
-gmsh.option.setNumber("Geometry.PointNumbers", 1);
-gmsh.option.setColor("Geometry.Points", 255, 165, 0);
-gmsh.option.setColor("General.Text", 255, 255, 255);
-gmsh.option.setColor("Mesh.Points", 255, 0, 0);
+gmsh.option.setNumber("Geometry.PointNumbers", 1)
+gmsh.option.setColor("Geometry.Points", 255, 165, 0)
+gmsh.option.setColor("General.Text", 255, 255, 255)
+gmsh.option.setColor("Mesh.Points", 255, 0, 0)
 
 # Note that color options are special: setting a color option of "X.Y"
 # actually sets the option "X.Color.Y".
 
-r, g, b, a = gmsh.option.getColor("Geometry.Points");
-gmsh.option.setColor("Geometry.Surfaces", r, g, b, a);
+r, g, b, a = gmsh.option.getColor("Geometry.Points")
+gmsh.option.setColor("Geometry.Surfaces", r, g, b, a)
 
 # Launch the GUI to see the effects of the color changes:
 
diff --git a/tutorial/python/t4.py b/tutorial/python/t4.py
index 291717ca71903a6dd7f169fa6595f0c688d312b7..9ed4d68a17a2f21630854b8b468e794e220698c1 100644
--- a/tutorial/python/t4.py
+++ b/tutorial/python/t4.py
@@ -15,87 +15,97 @@ gmsh.option.setNumber("General.Terminal", 1)
 gmsh.model.add("t4")
 
 cm = 1e-02
-e1 = 4.5 * cm; e2 = 6 * cm / 2; e3 =  5 * cm / 2
-h1 = 5 * cm; h2 = 10 * cm; h3 = 5 * cm; h4 = 2 * cm; h5 = 4.5 * cm
-R1 = 1 * cm; R2 = 1.5 * cm; r = 1 * cm
+e1 = 4.5 * cm
+e2 = 6 * cm / 2
+e3 = 5 * cm / 2
+h1 = 5 * cm
+h2 = 10 * cm
+h3 = 5 * cm
+h4 = 2 * cm
+h5 = 4.5 * cm
+R1 = 1 * cm
+R2 = 1.5 * cm
+r = 1 * cm
 Lc1 = 0.01
 Lc2 = 0.003
 
+
 def hypot(a, b):
     return math.sqrt(a * a + b * b)
 
-ccos = (-h5*R1 + e2 * hypot(h5, hypot(e2, R1))) / (h5*h5 + e2*e2)
-ssin = math.sqrt(1 - ccos*ccos)
+
+ccos = (-h5 * R1 + e2 * hypot(h5, hypot(e2, R1))) / (h5 * h5 + e2 * e2)
+ssin = math.sqrt(1 - ccos * ccos)
 
 # We start by defining some points and some lines. To make the code shorter we
 # can redefine a namespace:
 factory = gmsh.model.geo
-factory.addPoint(-e1-e2, 0    , 0, Lc1, 1)
-factory.addPoint(-e1-e2, h1   , 0, Lc1, 2)
-factory.addPoint(-e3-r , h1   , 0, Lc2, 3)
-factory.addPoint(-e3-r , h1+r , 0, Lc2, 4)
-factory.addPoint(-e3   , h1+r , 0, Lc2, 5)
-factory.addPoint(-e3   , h1+h2, 0, Lc1, 6)
-factory.addPoint( e3   , h1+h2, 0, Lc1, 7)
-factory.addPoint( e3   , h1+r , 0, Lc2, 8)
-factory.addPoint( e3+r , h1+r , 0, Lc2, 9)
-factory.addPoint( e3+r , h1   , 0, Lc2, 10)
-factory.addPoint( e1+e2, h1   , 0, Lc1, 11)
-factory.addPoint( e1+e2, 0    , 0, Lc1, 12)
-factory.addPoint( e2   , 0    , 0, Lc1, 13)
-
-factory.addPoint( R1 / ssin, h5+R1*ccos, 0, Lc2, 14)
-factory.addPoint( 0        , h5        , 0, Lc2, 15)
-factory.addPoint(-R1 / ssin, h5+R1*ccos, 0, Lc2, 16)
-factory.addPoint(-e2       , 0.0       , 0, Lc1, 17)
-
-factory.addPoint(-R2 , h1+h3   , 0, Lc2, 18)
-factory.addPoint(-R2 , h1+h3+h4, 0, Lc2, 19)
-factory.addPoint( 0  , h1+h3+h4, 0, Lc2, 20)
-factory.addPoint( R2 , h1+h3+h4, 0, Lc2, 21)
-factory.addPoint( R2 , h1+h3   , 0, Lc2, 22)
-factory.addPoint( 0  , h1+h3   , 0, Lc2, 23)
-
-factory.addPoint( 0, h1+h3+h4+R2, 0, Lc2, 24)
-factory.addPoint( 0, h1+h3-R2,    0, Lc2, 25)
-
-factory.addLine(1 , 17, 1)
+factory.addPoint(-e1 - e2, 0, 0, Lc1, 1)
+factory.addPoint(-e1 - e2, h1, 0, Lc1, 2)
+factory.addPoint(-e3 - r, h1, 0, Lc2, 3)
+factory.addPoint(-e3 - r, h1 + r, 0, Lc2, 4)
+factory.addPoint(-e3, h1 + r, 0, Lc2, 5)
+factory.addPoint(-e3, h1 + h2, 0, Lc1, 6)
+factory.addPoint(e3, h1 + h2, 0, Lc1, 7)
+factory.addPoint(e3, h1 + r, 0, Lc2, 8)
+factory.addPoint(e3 + r, h1 + r, 0, Lc2, 9)
+factory.addPoint(e3 + r, h1, 0, Lc2, 10)
+factory.addPoint(e1 + e2, h1, 0, Lc1, 11)
+factory.addPoint(e1 + e2, 0, 0, Lc1, 12)
+factory.addPoint(e2, 0, 0, Lc1, 13)
+
+factory.addPoint(R1 / ssin, h5 + R1 * ccos, 0, Lc2, 14)
+factory.addPoint(0, h5, 0, Lc2, 15)
+factory.addPoint(-R1 / ssin, h5 + R1 * ccos, 0, Lc2, 16)
+factory.addPoint(-e2, 0.0, 0, Lc1, 17)
+
+factory.addPoint(-R2, h1 + h3, 0, Lc2, 18)
+factory.addPoint(-R2, h1 + h3 + h4, 0, Lc2, 19)
+factory.addPoint(0, h1 + h3 + h4, 0, Lc2, 20)
+factory.addPoint(R2, h1 + h3 + h4, 0, Lc2, 21)
+factory.addPoint(R2, h1 + h3, 0, Lc2, 22)
+factory.addPoint(0, h1 + h3, 0, Lc2, 23)
+
+factory.addPoint(0, h1 + h3 + h4 + R2, 0, Lc2, 24)
+factory.addPoint(0, h1 + h3 - R2, 0, Lc2, 25)
+
+factory.addLine(1, 17, 1)
 factory.addLine(17, 16, 2)
 
 # Gmsh provides other curve primitives than straight lines: splines, B-splines,
 # circle arcs, ellipse arcs, etc. Here we define a new circle arc, starting at
 # point 14 and ending at point 16, with the circle's center being the point 15:
-factory.addCircleArc(14,15,16, 3)
+factory.addCircleArc(14, 15, 16, 3)
 
 # Note that, in Gmsh, circle arcs should always be smaller than Pi. The
 # OpenCASCADE geometry kernel does not have this limitation.
 
 # We can then define additional lines and circles, as well as a new surface:
-factory.addLine(14,13, 4)
-factory.addLine(13,12, 5)
-factory.addLine(12,11, 6)
-factory.addLine(11,10, 7)
-factory.addCircleArc(8,9,10, 8)
-factory.addLine(8,7, 9)
-factory.addLine(7,6, 10)
-factory.addLine(6,5, 11)
-factory.addCircleArc(3,4,5, 12)
-factory.addLine(3,2, 13)
-factory.addLine(2,1, 14)
-factory.addLine(18,19, 15)
-factory.addCircleArc(21,20,24, 16)
-factory.addCircleArc(24,20,19, 17)
-factory.addCircleArc(18,23,25, 18)
-factory.addCircleArc(25,23,22, 19)
-factory.addLine(21,22, 20)
-
-factory.addCurveLoop([17,-15,18,19,-20,16], 21)
+factory.addLine(14, 13, 4)
+factory.addLine(13, 12, 5)
+factory.addLine(12, 11, 6)
+factory.addLine(11, 10, 7)
+factory.addCircleArc(8, 9, 10, 8)
+factory.addLine(8, 7, 9)
+factory.addLine(7, 6, 10)
+factory.addLine(6, 5, 11)
+factory.addCircleArc(3, 4, 5, 12)
+factory.addLine(3, 2, 13)
+factory.addLine(2, 1, 14)
+factory.addLine(18, 19, 15)
+factory.addCircleArc(21, 20, 24, 16)
+factory.addCircleArc(24, 20, 19, 17)
+factory.addCircleArc(18, 23, 25, 18)
+factory.addCircleArc(25, 23, 22, 19)
+factory.addLine(21, 22, 20)
+
+factory.addCurveLoop([17, -15, 18, 19, -20, 16], 21)
 factory.addPlaneSurface([21], 22)
 
 # But we still need to define the exterior surface. Since this surface has a
 # hole, its definition now requires two curves loops:
-factory.addCurveLoop([11,-12,13,14,1,2,-3,4,5,6,7,-8,9,10], 23)
-factory.addPlaneSurface([23,21], 24)
+factory.addCurveLoop([11, -12, 13, 14, 1, 2, -3, 4, 5, 6, 7, -8, 9, 10], 23)
+factory.addPlaneSurface([23, 21], 24)
 
 # As a general rule, if a surface has N holes, it is defined by N+1 curve loops:
 # the first loop defines the exterior boundary; the other loops define the
@@ -121,8 +131,7 @@ gmsh.view.addListDataString(v, [0, 0.11, 0], ["Hole"],
 # a `@' symbol in the form `widthxheight' (if one of `width' or `height' is
 # zero, natural scaling is used; if both are zero, original image dimensions in
 # pixels are used):
-gmsh.view.addListDataString(v, [0, 0.09, 0],
-                            ["file://../t4_image.png@0.01x0"],
+gmsh.view.addListDataString(v, [0, 0.09, 0], ["file://../t4_image.png@0.01x0"],
                             ["Align", "Center"])
 
 # The 3D orientation of the image can be specified by proving the direction
@@ -149,15 +158,16 @@ gmsh.view.addListDataString(v, [150, -7], ["file://../t4_image.png@20x0"])
 # events, here to print some messages to the console:
 gmsh.option.setString("View[0].DoubleClickedCommand",
                       "Printf('View[0] has been double-clicked!');")
-gmsh.option.setString("Geometry.DoubleClickedLineCommand",
-                      "Printf('Curve %g has been double-clicked!', "
-                      "Geometry.DoubleClickedEntityTag);")
+gmsh.option.setString(
+    "Geometry.DoubleClickedLineCommand",
+    "Printf('Curve %g has been double-clicked!', "
+    "Geometry.DoubleClickedEntityTag);")
 
 # We can also change the color of some entities:
-gmsh.model.setColor([(2, 22)], 127, 127, 127) # Gray50
-gmsh.model.setColor([(2, 24)], 160, 32, 240) # Purple
-gmsh.model.setColor([(1, i) for i in range(1, 15)], 255, 0, 0) # Red
-gmsh.model.setColor([(1, i) for i in range(15, 21)], 255, 255, 0) # Yellow
+gmsh.model.setColor([(2, 22)], 127, 127, 127)  # Gray50
+gmsh.model.setColor([(2, 24)], 160, 32, 240)  # Purple
+gmsh.model.setColor([(1, i) for i in range(1, 15)], 255, 0, 0)  # Red
+gmsh.model.setColor([(1, i) for i in range(15, 21)], 255, 255, 0)  # Yellow
 
 gmsh.model.mesh.generate(2)
 
diff --git a/tutorial/python/t5.py b/tutorial/python/t5.py
index 1d04488ee709ae4a8437a74c0e2fd7cfeb81b361..7a1ab2d364e9f3c2b651861ad6ea21625de9e5c0 100644
--- a/tutorial/python/t5.py
+++ b/tutorial/python/t5.py
@@ -40,93 +40,103 @@ lcar3 = .055
 #
 # See `t10.py' for more information about mesh sizes.
 
-
 # We proceed by defining some elementary entities describing a truncated cube:
 
-gmsh.model.geo.addPoint(0.5,0.5,0.5, lcar2, 1)
-gmsh.model.geo.addPoint(0.5,0.5,0, lcar1, 2)
-gmsh.model.geo.addPoint(0,0.5,0.5, lcar1, 3)
-gmsh.model.geo.addPoint(0,0,0.5, lcar1, 4)
-gmsh.model.geo.addPoint(0.5,0,0.5, lcar1, 5)
-gmsh.model.geo.addPoint(0.5,0,0, lcar1, 6)
-gmsh.model.geo.addPoint(0,0.5,0, lcar1, 7)
-gmsh.model.geo.addPoint(0,1,0, lcar1, 8)
-gmsh.model.geo.addPoint(1,1,0, lcar1, 9)
-gmsh.model.geo.addPoint(0,0,1, lcar1, 10)
-gmsh.model.geo.addPoint(0,1,1, lcar1, 11)
-gmsh.model.geo.addPoint(1,1,1, lcar1, 12)
-gmsh.model.geo.addPoint(1,0,1, lcar1, 13)
-gmsh.model.geo.addPoint(1,0,0, lcar1, 14)
-
-gmsh.model.geo.addLine(8,9, 1);   gmsh.model.geo.addLine(9,12, 2)
-gmsh.model.geo.addLine(12,11, 3); gmsh.model.geo.addLine(11,8, 4)
-gmsh.model.geo.addLine(9,14, 5);  gmsh.model.geo.addLine(14,13, 6)
-gmsh.model.geo.addLine(13,12, 7); gmsh.model.geo.addLine(11,10, 8)
-gmsh.model.geo.addLine(10,13, 9); gmsh.model.geo.addLine(10,4, 10)
-gmsh.model.geo.addLine(4,5, 11);  gmsh.model.geo.addLine(5,6, 12)
-gmsh.model.geo.addLine(6,2, 13);  gmsh.model.geo.addLine(2,1, 14)
-gmsh.model.geo.addLine(1,3, 15);  gmsh.model.geo.addLine(3,7, 16)
-gmsh.model.geo.addLine(7,2, 17);  gmsh.model.geo.addLine(3,4, 18)
-gmsh.model.geo.addLine(5,1, 19);  gmsh.model.geo.addLine(7,8, 20)
-gmsh.model.geo.addLine(6,14, 21)
-
-gmsh.model.geo.addCurveLoop([-11,-19,-15,-18], 22)
+gmsh.model.geo.addPoint(0.5, 0.5, 0.5, lcar2, 1)
+gmsh.model.geo.addPoint(0.5, 0.5, 0, lcar1, 2)
+gmsh.model.geo.addPoint(0, 0.5, 0.5, lcar1, 3)
+gmsh.model.geo.addPoint(0, 0, 0.5, lcar1, 4)
+gmsh.model.geo.addPoint(0.5, 0, 0.5, lcar1, 5)
+gmsh.model.geo.addPoint(0.5, 0, 0, lcar1, 6)
+gmsh.model.geo.addPoint(0, 0.5, 0, lcar1, 7)
+gmsh.model.geo.addPoint(0, 1, 0, lcar1, 8)
+gmsh.model.geo.addPoint(1, 1, 0, lcar1, 9)
+gmsh.model.geo.addPoint(0, 0, 1, lcar1, 10)
+gmsh.model.geo.addPoint(0, 1, 1, lcar1, 11)
+gmsh.model.geo.addPoint(1, 1, 1, lcar1, 12)
+gmsh.model.geo.addPoint(1, 0, 1, lcar1, 13)
+gmsh.model.geo.addPoint(1, 0, 0, lcar1, 14)
+
+gmsh.model.geo.addLine(8, 9, 1)
+gmsh.model.geo.addLine(9, 12, 2)
+gmsh.model.geo.addLine(12, 11, 3)
+gmsh.model.geo.addLine(11, 8, 4)
+gmsh.model.geo.addLine(9, 14, 5)
+gmsh.model.geo.addLine(14, 13, 6)
+gmsh.model.geo.addLine(13, 12, 7)
+gmsh.model.geo.addLine(11, 10, 8)
+gmsh.model.geo.addLine(10, 13, 9)
+gmsh.model.geo.addLine(10, 4, 10)
+gmsh.model.geo.addLine(4, 5, 11)
+gmsh.model.geo.addLine(5, 6, 12)
+gmsh.model.geo.addLine(6, 2, 13)
+gmsh.model.geo.addLine(2, 1, 14)
+gmsh.model.geo.addLine(1, 3, 15)
+gmsh.model.geo.addLine(3, 7, 16)
+gmsh.model.geo.addLine(7, 2, 17)
+gmsh.model.geo.addLine(3, 4, 18)
+gmsh.model.geo.addLine(5, 1, 19)
+gmsh.model.geo.addLine(7, 8, 20)
+gmsh.model.geo.addLine(6, 14, 21)
+
+gmsh.model.geo.addCurveLoop([-11, -19, -15, -18], 22)
 gmsh.model.geo.addPlaneSurface([22], 23)
-gmsh.model.geo.addCurveLoop([16,17,14,15], 24)
+gmsh.model.geo.addCurveLoop([16, 17, 14, 15], 24)
 gmsh.model.geo.addPlaneSurface([24], 25)
-gmsh.model.geo.addCurveLoop([-17,20,1,5,-21,13], 26)
+gmsh.model.geo.addCurveLoop([-17, 20, 1, 5, -21, 13], 26)
 gmsh.model.geo.addPlaneSurface([26], 27)
-gmsh.model.geo.addCurveLoop([-4,-1,-2,-3], 28)
+gmsh.model.geo.addCurveLoop([-4, -1, -2, -3], 28)
 gmsh.model.geo.addPlaneSurface([28], 29)
-gmsh.model.geo.addCurveLoop([-7,2,-5,-6], 30)
+gmsh.model.geo.addCurveLoop([-7, 2, -5, -6], 30)
 gmsh.model.geo.addPlaneSurface([30], 31)
-gmsh.model.geo.addCurveLoop([6,-9,10,11,12,21], 32)
+gmsh.model.geo.addCurveLoop([6, -9, 10, 11, 12, 21], 32)
 gmsh.model.geo.addPlaneSurface([32], 33)
-gmsh.model.geo.addCurveLoop([7,3,8,9], 34)
+gmsh.model.geo.addCurveLoop([7, 3, 8, 9], 34)
 gmsh.model.geo.addPlaneSurface([34], 35)
-gmsh.model.geo.addCurveLoop([-10,18,-16,-20,4,-8], 36)
+gmsh.model.geo.addCurveLoop([-10, 18, -16, -20, 4, -8], 36)
 gmsh.model.geo.addPlaneSurface([36], 37)
-gmsh.model.geo.addCurveLoop([-14,-13,-12,19], 38)
+gmsh.model.geo.addCurveLoop([-14, -13, -12, 19], 38)
 gmsh.model.geo.addPlaneSurface([38], 39)
 
 shells = []
 
-sl = gmsh.model.geo.addSurfaceLoop([35,31,29,37,33,23,39,25,27])
+sl = gmsh.model.geo.addSurfaceLoop([35, 31, 29, 37, 33, 23, 39, 25, 27])
 shells.append(sl)
 
+
 def cheeseHole(x, y, z, r, lc, shells):
     # This function will create a spherical hole in a volume. We don't specify
     # tags manually, and let the functions return them automatically:
 
-    p1 = gmsh.model.geo.addPoint(x,  y,  z,   lc)
-    p2 = gmsh.model.geo.addPoint(x+r,y,  z,   lc)
-    p3 = gmsh.model.geo.addPoint(x,  y+r,z,   lc)
-    p4 = gmsh.model.geo.addPoint(x,  y,  z+r, lc)
-    p5 = gmsh.model.geo.addPoint(x-r,y,  z,   lc)
-    p6 = gmsh.model.geo.addPoint(x,  y-r,z,   lc)
-    p7 = gmsh.model.geo.addPoint(x,  y,  z-r, lc)
-
-    c1 = gmsh.model.geo.addCircleArc(p2,p1,p7)
-    c2 = gmsh.model.geo.addCircleArc(p7,p1,p5)
-    c3 = gmsh.model.geo.addCircleArc(p5,p1,p4)
-    c4 = gmsh.model.geo.addCircleArc(p4,p1,p2)
-    c5 = gmsh.model.geo.addCircleArc(p2,p1,p3)
-    c6 = gmsh.model.geo.addCircleArc(p3,p1,p5)
-    c7 = gmsh.model.geo.addCircleArc(p5,p1,p6)
-    c8 = gmsh.model.geo.addCircleArc(p6,p1,p2)
-    c9 = gmsh.model.geo.addCircleArc(p7,p1,p3)
-    c10 = gmsh.model.geo.addCircleArc(p3,p1,p4)
-    c11 = gmsh.model.geo.addCircleArc(p4,p1,p6)
-    c12 = gmsh.model.geo.addCircleArc(p6,p1,p7)
-
-    l1 = gmsh.model.geo.addCurveLoop([c5,c10,c4])
-    l2 = gmsh.model.geo.addCurveLoop([c9,-c5,c1])
-    l3 = gmsh.model.geo.addCurveLoop([c12,-c8,-c1])
-    l4 = gmsh.model.geo.addCurveLoop([c8,-c4,c11])
-    l5 = gmsh.model.geo.addCurveLoop([-c10,c6,c3])
-    l6 = gmsh.model.geo.addCurveLoop([-c11,-c3,c7])
-    l7 = gmsh.model.geo.addCurveLoop([-c2,-c7,-c12])
-    l8 = gmsh.model.geo.addCurveLoop([-c6,-c9,c2])
+    p1 = gmsh.model.geo.addPoint(x, y, z, lc)
+    p2 = gmsh.model.geo.addPoint(x + r, y, z, lc)
+    p3 = gmsh.model.geo.addPoint(x, y + r, z, lc)
+    p4 = gmsh.model.geo.addPoint(x, y, z + r, lc)
+    p5 = gmsh.model.geo.addPoint(x - r, y, z, lc)
+    p6 = gmsh.model.geo.addPoint(x, y - r, z, lc)
+    p7 = gmsh.model.geo.addPoint(x, y, z - r, lc)
+
+    c1 = gmsh.model.geo.addCircleArc(p2, p1, p7)
+    c2 = gmsh.model.geo.addCircleArc(p7, p1, p5)
+    c3 = gmsh.model.geo.addCircleArc(p5, p1, p4)
+    c4 = gmsh.model.geo.addCircleArc(p4, p1, p2)
+    c5 = gmsh.model.geo.addCircleArc(p2, p1, p3)
+    c6 = gmsh.model.geo.addCircleArc(p3, p1, p5)
+    c7 = gmsh.model.geo.addCircleArc(p5, p1, p6)
+    c8 = gmsh.model.geo.addCircleArc(p6, p1, p2)
+    c9 = gmsh.model.geo.addCircleArc(p7, p1, p3)
+    c10 = gmsh.model.geo.addCircleArc(p3, p1, p4)
+    c11 = gmsh.model.geo.addCircleArc(p4, p1, p6)
+    c12 = gmsh.model.geo.addCircleArc(p6, p1, p7)
+
+    l1 = gmsh.model.geo.addCurveLoop([c5, c10, c4])
+    l2 = gmsh.model.geo.addCurveLoop([c9, -c5, c1])
+    l3 = gmsh.model.geo.addCurveLoop([c12, -c8, -c1])
+    l4 = gmsh.model.geo.addCurveLoop([c8, -c4, c11])
+    l5 = gmsh.model.geo.addCurveLoop([-c10, c6, c3])
+    l6 = gmsh.model.geo.addCurveLoop([-c11, -c3, c7])
+    l7 = gmsh.model.geo.addCurveLoop([-c2, -c7, -c12])
+    l8 = gmsh.model.geo.addCurveLoop([-c6, -c9, c2])
 
     # We need non-plane surfaces to define the spherical holes. Here we use the
     # `gmsh.model.geo.addSurfaceFilling()' function, which can be used for
@@ -154,8 +164,12 @@ def cheeseHole(x, y, z, r, lc, shells):
     shells.append(sl)
     return v
 
+
 # We create five holes in the cube:
-x = 0; y = 0.75; z = 0; r = 0.09
+x = 0
+y = 0.75
+z = 0
+r = 0.09
 for t in range(1, 6):
     x += 0.166
     z += 0.166
@@ -184,8 +198,9 @@ gmsh.model.geo.synchronize()
 # gmsh.option.setNumber("Mesh.MeshOnlyVisible", 1)
 
 # Meshing algorithms can changed globally using options:
-gmsh.option.setNumber("Mesh.Algorithm", 6) # 2D algorithm set to Frontal-Delaunay
-gmsh.option.setNumber("Mesh.Algorithm3D", 1) # 3D algorithm set to Delaunay
+gmsh.option.setNumber("Mesh.Algorithm",
+                      6)  # 2D algorithm set to Frontal-Delaunay
+gmsh.option.setNumber("Mesh.Algorithm3D", 1)  # 3D algorithm set to Delaunay
 
 # They can also be set for individual surfaces, e.g. for using `MeshAdapt' on
 # surface 1:
diff --git a/tutorial/python/t6.py b/tutorial/python/t6.py
index c12de9f0a4e6d1554e34bad7a24563367ac1ddac..c2b13d6af6194f0d147313984d73049269c36132 100644
--- a/tutorial/python/t6.py
+++ b/tutorial/python/t6.py
@@ -17,7 +17,7 @@ gmsh.model.add("t6")
 # Copied from t1.py...
 lc = 1e-2
 gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
-gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
+gmsh.model.geo.addPoint(.1, 0, 0, lc, 2)
 gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
 gmsh.model.geo.addPoint(0, .3, 0, lc, 4)
 gmsh.model.geo.addLine(1, 2, 1)
@@ -28,7 +28,7 @@ gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1)
 gmsh.model.geo.addPlaneSurface([1], 1)
 
 # Delete the surface and the left line, and replace the line with 3 new ones:
-gmsh.model.geo.remove([[2,1], [1,4]])
+gmsh.model.geo.remove([[2, 1], [1, 4]])
 
 p1 = gmsh.model.geo.addPoint(-0.05, 0.05, 0, lc)
 p2 = gmsh.model.geo.addPoint(-0.05, 0.1, 0, lc)
@@ -63,7 +63,7 @@ gmsh.model.geo.mesh.setTransfiniteCurve(3, 30, "Progression", 1.2)
 # nodes on the boundary using a structured grid. If the surface has more than 4
 # corner points, the corners of the transfinite interpolation have to be
 # specified by hand:
-gmsh.model.geo.mesh.setTransfiniteSurface(1, "Left", [1,2,3,4])
+gmsh.model.geo.mesh.setTransfiniteSurface(1, "Left", [1, 2, 3, 4])
 
 # To create quadrangles instead of triangles, one can use the `setRecombine'
 # constraint:
@@ -82,7 +82,7 @@ gmsh.model.geo.addLine(10, 7, 12)
 gmsh.model.geo.addLine(7, 8, 13)
 gmsh.model.geo.addCurveLoop([13, 10, 11, 12], 14)
 gmsh.model.geo.addPlaneSurface([14], 15)
-for i in range(10,14):
+for i in range(10, 14):
     gmsh.model.geo.mesh.setTransfiniteCurve(i, 10)
 gmsh.model.geo.mesh.setTransfiniteSurface(15)
 
diff --git a/tutorial/python/t7.py b/tutorial/python/t7.py
index 9939441a151a29e467e5ceb1eb4fbf66aa2ac4bc..05b3b7318adc0e47e68f826525623b3f48ed0979 100644
--- a/tutorial/python/t7.py
+++ b/tutorial/python/t7.py
@@ -18,7 +18,7 @@ gmsh.option.setNumber("General.Terminal", 1)
 # Create a simple rectangular geometry
 lc = 1e-2
 gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
-gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
+gmsh.model.geo.addPoint(.1, 0, 0, lc, 2)
 gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
 gmsh.model.geo.addPoint(0, .3, 0, lc, 4)
 gmsh.model.geo.addLine(1, 2, 1)
diff --git a/tutorial/python/t8.py b/tutorial/python/t8.py
index 6bb74933c31307ab4b9072945129d9bb8ddf1765..ce1d457345b5837d4430b95ad61ca3b6eb3e0b70 100644
--- a/tutorial/python/t8.py
+++ b/tutorial/python/t8.py
@@ -18,7 +18,7 @@ gmsh.option.setNumber("General.Terminal", 1)
 # We first create a simple geometry
 lc = 1e-2
 gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
-gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
+gmsh.model.geo.addPoint(.1, 0, 0, lc, 2)
 gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
 gmsh.model.geo.addPoint(0, .3, 0, lc, 4)
 gmsh.model.geo.addLine(1, 2, 1)
@@ -34,7 +34,7 @@ gmsh.model.geo.synchronize()
 path = os.path.dirname(os.path.abspath(__file__))
 gmsh.merge(os.path.join(path, '..', 'view1.pos'))
 gmsh.merge(os.path.join(path, '..', 'view1.pos'))
-gmsh.merge(os.path.join(path, '..', 'view4.pos')) # contains 2 views inside
+gmsh.merge(os.path.join(path, '..', 'view4.pos'))  # contains 2 views inside
 
 # Gmsh can read post-processing views in various formats. Here the `view1.pos'
 # and `view4.pos' files are in the Gmsh "parsed" format, which is interpreted by
@@ -63,7 +63,7 @@ set_color("General.Foreground", black)
 set_color("General.Text", black)
 
 gmsh.option.setNumber("General.Orthographic", 0)
-gmsh.option.setNumber("General.Axes", 0);
+gmsh.option.setNumber("General.Axes", 0)
 gmsh.option.setNumber("General.SmallAxes", 0)
 
 # Show the GUI
@@ -79,6 +79,7 @@ view_tags = [v0, v1, v2, v3] = gmsh.view.getTags()
 # View name format helper function: returns "View[index]." for a given view tag
 view_fmt = lambda v_tag: "View[" + str(gmsh.view.getIndex(v_tag)) + "]."
 
+
 # Option setter
 def set_opt(name, val):
     # if it's a string, call the string method
@@ -93,6 +94,7 @@ def set_opt(name, val):
         print("error: set_opt is only meant for numbers and strings, aborting")
         quit(1)
 
+
 # We'll use this helper function for our views
 set_view_opt = lambda v_tag, name, val: set_opt(view_fmt(v_tag) + name, val)
 
@@ -136,7 +138,7 @@ set_view_opt(v3, "Visible", 0)
 # saved to disk as an image, and multiple frames can be encoded to form a
 # movie. Below is an example of such a custom animation.
 
-t = 0 # Initial step
+t = 0  # Initial step
 
 for num in range(1, 4):
 
diff --git a/tutorial/python/x1.py b/tutorial/python/x1.py
index d87d76d04428763b0068927e8c2a96bbc44b8931..6e1c7e63cd741a883e5ffcff855b560869406259 100644
--- a/tutorial/python/x1.py
+++ b/tutorial/python/x1.py
@@ -104,12 +104,13 @@ for e in entities:
     # * Is the entity a partition entity? If so, what is its parent entity?
     partitions = gmsh.model.getPartitions(e[0], e[1])
     if len(partitions):
-           print(" - Partition tags: " + str(partitions) +
-                 " - parent entity " + str(gmsh.model.getParent(e[0], e[1])))
+        print(" - Partition tags: " + str(partitions) + " - parent entity " +
+              str(gmsh.model.getParent(e[0], e[1])))
 
     # * List all types of elements making up the mesh of the entity:
     for t in elemTypes:
-        name, dim, order, numv, parv, _ = gmsh.model.mesh.getElementProperties(t)
+        name, dim, order, numv, parv, _ = gmsh.model.mesh.getElementProperties(
+            t)
         print(" - Element type: " + name + ", order " + str(order) + " (" +
               str(numv) + " nodes in param coord: " + str(parv) + ")")
 
diff --git a/tutorial/python/x2.py b/tutorial/python/x2.py
index 8fe0e11a89a138ba935997b98de66fce6254b666..69370b01e252c0e8a5ecd6dc9a2cec4748777c97 100644
--- a/tutorial/python/x2.py
+++ b/tutorial/python/x2.py
@@ -30,9 +30,11 @@ gmsh.model.add("x2")
 # We will create the terrain surface mesh from N x N input data points:
 N = 100
 
+
 # Helper function to return a node tag given two indices i and j:
 def tag(i, j):
-    return (N+1) * i + j + 1
+    return (N + 1) * i + j + 1
+
 
 # The x, y, z coordinates of all the nodes:
 coords = []
@@ -46,7 +48,7 @@ tris = []
 
 # The connectivities of the line elements on the 4 boundaries (2 node tags
 # for each line element):
-lin = [[],[],[],[]]
+lin = [[], [], [], []]
 
 # The connectivities of the point elements on the 4 corners (1 node tag for each
 # point element):
@@ -55,26 +57,29 @@ pnt = [tag(0, 0), tag(N, 0), tag(N, N), tag(0, N)]
 for i in range(N + 1):
     for j in range(N + 1):
         nodes.append(tag(i, j))
-        coords.extend([float(i)/N, float(j)/N, 0.05 * math.sin(10*float(i+j) / N)])
+        coords.extend([
+            float(i) / N,
+            float(j) / N, 0.05 * math.sin(10 * float(i + j) / N)
+        ])
         if i > 0 and j > 0:
-            tris.extend([tag(i-1, j-1), tag(i, j-1), tag(i-1, j)])
-            tris.extend([tag(i, j-1), tag(i, j), tag(i-1, j)])
+            tris.extend([tag(i - 1, j - 1), tag(i, j - 1), tag(i - 1, j)])
+            tris.extend([tag(i, j - 1), tag(i, j), tag(i - 1, j)])
         if (i == 0 or i == N) and j > 0:
-            lin[3 if i == 0 else 1].extend([tag(i, j-1), tag(i, j)])
+            lin[3 if i == 0 else 1].extend([tag(i, j - 1), tag(i, j)])
         if (j == 0 or j == N) and i > 0:
-            lin[0 if j == 0 else 2].extend([tag(i-1, j), tag(i, j)])
+            lin[0 if j == 0 else 2].extend([tag(i - 1, j), tag(i, j)])
 
 # Create 4 discrete points for the 4 corners of the terrain surface:
 for i in range(4):
-    gmsh.model.addDiscreteEntity(0, i+1)
-gmsh.model.setCoordinates(1, 0, 0, coords[3*tag(0,0)-1])
-gmsh.model.setCoordinates(2, 1, 0, coords[3*tag(N,0)-1])
-gmsh.model.setCoordinates(3, 1, 1, coords[3*tag(N,N)-1])
-gmsh.model.setCoordinates(4, 0, 1, coords[3*tag(0,N)-1])
+    gmsh.model.addDiscreteEntity(0, i + 1)
+gmsh.model.setCoordinates(1, 0, 0, coords[3 * tag(0, 0) - 1])
+gmsh.model.setCoordinates(2, 1, 0, coords[3 * tag(N, 0) - 1])
+gmsh.model.setCoordinates(3, 1, 1, coords[3 * tag(N, N) - 1])
+gmsh.model.setCoordinates(4, 0, 1, coords[3 * tag(0, N) - 1])
 
 # Create 4 discrete bounding curves, with their (CAD) boundary points:
 for i in range(4):
-    gmsh.model.addDiscreteEntity(1, i+1, [i+1 , i+2 if i < 3 else 1])
+    gmsh.model.addDiscreteEntity(1, i + 1, [i + 1, i + 2 if i < 3 else 1])
 
 # Create one discrete surface, with its bounding curves:
 gmsh.model.addDiscreteEntity(2, 1, [1, 2, -3, -4])
@@ -118,15 +123,15 @@ c10 = gmsh.model.geo.addLine(p1, 1)
 c11 = gmsh.model.geo.addLine(p2, 2)
 c12 = gmsh.model.geo.addLine(p3, 3)
 c13 = gmsh.model.geo.addLine(p4, 4)
-ll1 = gmsh.model.geo.addCurveLoop([c1,c2,c3,c4])
+ll1 = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4])
 s1 = gmsh.model.geo.addPlaneSurface([ll1])
-ll3 = gmsh.model.geo.addCurveLoop([c1,c11,-1,-c10])
+ll3 = gmsh.model.geo.addCurveLoop([c1, c11, -1, -c10])
 s3 = gmsh.model.geo.addPlaneSurface([ll3])
-ll4 = gmsh.model.geo.addCurveLoop([c2,c12,-2,-c11])
+ll4 = gmsh.model.geo.addCurveLoop([c2, c12, -2, -c11])
 s4 = gmsh.model.geo.addPlaneSurface([ll4])
-ll5 = gmsh.model.geo.addCurveLoop([c3,c13,3,-c12])
+ll5 = gmsh.model.geo.addCurveLoop([c3, c13, 3, -c12])
 s5 = gmsh.model.geo.addPlaneSurface([ll5])
-ll6 = gmsh.model.geo.addCurveLoop([c4,c10,4,-c13])
+ll6 = gmsh.model.geo.addCurveLoop([c4, c10, 4, -c13])
 s6 = gmsh.model.geo.addPlaneSurface([ll6])
 sl1 = gmsh.model.geo.addSurfaceLoop([s1, s3, s4, s5, s6, 1])
 v1 = gmsh.model.geo.addVolume([sl1])