diff --git a/HySoP/hysop/fields/vector.py b/HySoP/hysop/fields/vector.py
index cce238fc58f12183ea3e06fbba6d3ea033fbf7a4..aa0087f73681254d1bc1f9ec80e6c7dfd4851ae6 100644
--- a/HySoP/hysop/fields/vector.py
+++ b/HySoP/hysop/fields/vector.py
@@ -90,7 +90,10 @@ class VectorField(DiscreteField):
         """
         if formula is not None:
             print 'initializing', self.name, '...'
-            v_formula = np.vectorize(formula)
+            if isinstance(formula, np.lib.function_base.vectorize):
+                v_formula = formula
+            else:
+                v_formula = np.vectorize(formula)
             if self.dimension == 3:
                 self.data[0][...], self.data[1][...], self.data[2][...] = \
                     v_formula(*(self.topology.mesh.coords + args))
diff --git a/HySoP/hysop/gpu/gpu_discrete.py b/HySoP/hysop/gpu/gpu_discrete.py
index a0155d9bee89bb3851de6a131dc630b6dbe29dcb..f87c87a82feb9746b14d227591013c70905c1c1e 100644
--- a/HySoP/hysop/gpu/gpu_discrete.py
+++ b/HySoP/hysop/gpu/gpu_discrete.py
@@ -103,7 +103,7 @@ class GPUVectorField(VectorField, GPUDiscreteField):
             print vfield.mem_size / (1024 ** 2), "MB)"
 
     @debug
-    def initialize(self, formula=None):
+    def initialize(self, formula=None, *args):
         """
         GPU data initialization.
         Performs the initialization from different ways if device not already
@@ -111,59 +111,58 @@ class GPUVectorField(VectorField, GPUDiscreteField):
           - with an OpenCL kernel,
           - with a python formula (as VectorField) and the copy data to device.
         """
-        if not self.data_on_device:
-            if self.initialization_kernel is None:
-                VectorField.initialize(self, formula)
-                if self.contains_data:
-                    if self.dimension == 2:
-                        self.data[0] = np.asarray(
-                            self.data[0],
-                            dtype=self.precision, order=ORDER)
-                        self.data[1] = np.asarray(
-                            self.data[1],
-                            dtype=self.precision, order=ORDER)
-                    else:
-                        self.data[0] = np.asarray(
-                            self.data[0],
-                            dtype=self.precision, order=ORDER)
-                        self.data[1] = np.asarray(
-                            self.data[1],
-                            dtype=self.precision, order=ORDER)
-                        self.data[2] = np.asarray(
-                            self.data[2],
-                            dtype=self.precision, order=ORDER)
-                    data, time = self.toDevice()
-                    if CL_PROFILE:
-                        print "Transfers : ", self.mem_size,
-                        print "Bytes transfered at ",
-                        print "{0:.3f} GBytes/sec".format(
-                            (self.mem_size * 1e-9) / time)
-            else:
-                coord_min = np.ones(4, dtype=self.precision)
-                mesh_size = np.ones(4, dtype=self.precision)
-                coord_min[0:self.dimension] = np.asarray(
-                    self.topology.mesh.origin,
-                    dtype=self.precision)
-                mesh_size[0:self.dimension] = np.asarray(
-                    self.topology.mesh.space_step,
-                    dtype=self.precision)
+        if self.initialization_kernel is None:
+            VectorField.initialize(self, formula, *args)
+            if self.contains_data:
                 if self.dimension == 2:
-                    evt = self.initialization_kernel(self.gpu_data[0],
-                                                     self.gpu_data[1],
-                                                     coord_min,
-                                                     mesh_size)
+                    self.data[0] = np.asarray(
+                        self.data[0],
+                        dtype=self.precision, order=ORDER)
+                    self.data[1] = np.asarray(
+                        self.data[1],
+                        dtype=self.precision, order=ORDER)
                 else:
-                    evt = self.initialization_kernel(self.gpu_data[0],
-                                                     self.gpu_data[1],
-                                                     self.gpu_data[2],
-                                                     coord_min,
-                                                     mesh_size)
+                    self.data[0] = np.asarray(
+                        self.data[0],
+                        dtype=self.precision, order=ORDER)
+                    self.data[1] = np.asarray(
+                        self.data[1],
+                        dtype=self.precision, order=ORDER)
+                    self.data[2] = np.asarray(
+                        self.data[2],
+                        dtype=self.precision, order=ORDER)
+                data, time = self.toDevice()
                 if CL_PROFILE:
-                    self.queue.finish()
-                    cTime = (evt.profile.end - evt.profile.start) * 1e-9
-                    print "Computing  : ", cTime, "sec"
-            self.contains_data = False
-            self.data_on_device = True
+                    print "Transfers : ", self.mem_size,
+                    print "Bytes transfered at ",
+                    print "{0:.3f} GBytes/sec".format(
+                        (self.mem_size * 1e-9) / time)
+        else:
+            coord_min = np.ones(4, dtype=self.precision)
+            mesh_size = np.ones(4, dtype=self.precision)
+            coord_min[0:self.dimension] = np.asarray(
+                self.topology.mesh.origin,
+                dtype=self.precision)
+            mesh_size[0:self.dimension] = np.asarray(
+                self.topology.mesh.space_step,
+                dtype=self.precision)
+            if self.dimension == 2:
+                evt = self.initialization_kernel(self.gpu_data[0],
+                                                 self.gpu_data[1],
+                                                 coord_min,
+                                                 mesh_size)
+            else:
+                evt = self.initialization_kernel(self.gpu_data[0],
+                                                 self.gpu_data[1],
+                                                 self.gpu_data[2],
+                                                 coord_min,
+                                                 mesh_size)
+            if CL_PROFILE:
+                self.queue.finish()
+                cTime = (evt.profile.end - evt.profile.start) * 1e-9
+                print "Computing  : ", cTime, "sec"
+        self.contains_data = False
+        self.data_on_device = True
 
     def toDevice(self):
         """
@@ -315,7 +314,7 @@ class GPUScalarField(ScalarField, GPUDiscreteField):
             print sfield.mem_size / (1024 ** 2), "MB)"
 
     @debug
-    def initialize(self, formula=None):
+    def initialize(self, formula=None, *args):
         """
         GPU data initialization.
         Performs the initialization from different ways if device not already
@@ -323,36 +322,35 @@ class GPUScalarField(ScalarField, GPUDiscreteField):
           - with an OpenCL kernel,
           - with a python formula (as VectorField) and the copy data to device.
         """
-        if not self.data_on_device:
-            if self.initialization_kernel is None:
-                ScalarField.initialize(self, formula)
-                if self.contains_data:
-                    self.data = np.asarray(self.data,
-                                           dtype=self.precision, order=ORDER)
-                    data, time = self.toDevice()
-                    if CL_PROFILE:
-                        print "Transfers : ", self.mem_size,
-                        print "Bytes transfered at ",
-                        print "{0:.3f} GBytes/sec".format(
-                            (self.mem_size * 1e-9) / time)
-            else:
-                coord_min = np.ones(4, dtype=self.precision)
-                mesh_size = np.ones(4, dtype=self.precision)
-                coord_min[0:self.dimension] = np.asarray(
-                    self.topology.mesh.origin,
-                    dtype=self.precision)
-                mesh_size[0:self.dimension] = np.asarray(
-                    self.topology.mesh.space_step,
-                    dtype=self.precision)
-                evt = self.initialization_kernel(self.gpu_data,
-                                                 coord_min,
-                                                 mesh_size)
+        if self.initialization_kernel is None:
+            ScalarField.initialize(self, formula, *args)
+            if self.contains_data:
+                self.data = np.asarray(self.data,
+                                       dtype=self.precision, order=ORDER)
+                data, time = self.toDevice()
                 if CL_PROFILE:
-                    self.queue.finish()
-                    cTime = (evt.profile.end - evt.profile.start) * 1e-9
-                    print "Computing  : ", cTime, "sec"
-            self.contains_data = False
-            self.data_on_device = True
+                    print "Transfers : ", self.mem_size,
+                    print "Bytes transfered at ",
+                    print "{0:.3f} GBytes/sec".format(
+                        (self.mem_size * 1e-9) / time)
+        else:
+            coord_min = np.ones(4, dtype=self.precision)
+            mesh_size = np.ones(4, dtype=self.precision)
+            coord_min[0:self.dimension] = np.asarray(
+                self.topology.mesh.origin,
+                dtype=self.precision)
+            mesh_size[0:self.dimension] = np.asarray(
+                self.topology.mesh.space_step,
+                dtype=self.precision)
+            evt = self.initialization_kernel(self.gpu_data,
+                                             coord_min,
+                                             mesh_size)
+            if CL_PROFILE:
+                self.queue.finish()
+                cTime = (evt.profile.end - evt.profile.start) * 1e-9
+                print "Computing  : ", cTime, "sec"
+        self.contains_data = False
+        self.data_on_device = True
 
     def toDevice(self):
         """
diff --git a/HySoP/hysop/operator/discrete/analytic.py b/HySoP/hysop/operator/discrete/analytic.py
index 31aeef6e63baa2ca76fb5be64e1eb9a6b4138065..575c9fd094603ebb1bd16798c06ec39d10b73783 100644
--- a/HySoP/hysop/operator/discrete/analytic.py
+++ b/HySoP/hysop/operator/discrete/analytic.py
@@ -3,7 +3,7 @@
 
 Apply user-defined formula to a set of discrete variables.
 """
-from parmepy.constants import np, debug, ORDER
+from parmepy.constants import np, debug
 from discrete import DiscreteOperator
 from parmepy.mpi import MPI
 
@@ -39,31 +39,7 @@ class Analytic_D(DiscreteOperator):
     def apply(self, t, dt, ite):
         self.compute_time = MPI.Wtime()
         for df in self.variables:
-            datatype = df.data[0].dtype
-            if df.isVector:
-                if df.topology.dim == 3:
-                    df.data[0][...], df.data[1][...], df.data[2][...] = \
-                        self.vformula(*(
-                            df.topology.mesh.coords + (t, dt, ite)))
-                    df.data[0][...] = np.asarray(df.data[0],
-                                                 dtype=datatype, order=ORDER)
-                    df.data[1][...] = np.asarray(df.data[1],
-                                                 dtype=datatype, order=ORDER)
-                    df.data[2][...] = np.asarray(df.data[2],
-                                                 dtype=datatype, order=ORDER)
-                elif df.topology.dim == 2:
-                    df.data[0][...], df.data[1][...] = \
-                        self.vformula(*(
-                            df.topology.mesh.coords + (t, dt, ite)))
-                    df.data[0][...] = np.asarray(df.data[0],
-                                                 dtype=datatype, order=ORDER)
-                    df.data[1][...] = np.asarray(df.data[1],
-                                                 dtype=datatype, order=ORDER)
-            else:
-                df.data[...] = \
-                    self.vformula(*(df.topology.mesh.coords + (t, dt, ite)))
-                df.data[...] = np.asarray(df.data,
-                                          dtype=datatype, order=ORDER)
+            df.initialize(self.vformula, t, dt, ite)
         for df in self.output:
             df.contains_data = True