diff --git a/hysop/fields/continuous_field.py b/hysop/fields/continuous_field.py
index 7f1ba42b639d6afbd859c2bc59a994db2d620541..a3603b10e9e56f0e427332b0cf274bb9fe25e806 100644
--- a/hysop/fields/continuous_field.py
+++ b/hysop/fields/continuous_field.py
@@ -10,7 +10,7 @@ import textwrap
 import sympy as sm
 from abc import ABCMeta, abstractmethod
 
-from hysop.constants         import HYSOP_REAL, HYSOP_BOOL, BoundaryCondition
+from hysop.constants         import HYSOP_REAL, HYSOP_BOOL, BoundaryCondition, DirectionLabels
 from hysop.tools.decorators  import debug
 from hysop.tools.types       import check_instance, first_not_None, to_tuple
 from hysop.tools.warning    import HysopWarning
@@ -30,7 +30,7 @@ class FieldContainerI(TaggedObject):
     """Common abstract interface for scalar and tensor-like fields."""
 
     @debug
-    def __new__(cls, domain, 
+    def __new__(cls, domain,
             name=None, nb_components=None, shape=None, is_vector=None, **kwds):
         """
         Create a FieldContainer on a specific domain.
@@ -52,7 +52,7 @@ class FieldContainerI(TaggedObject):
                   ((nb_components is not None) and (nb_components > 1))):
                 nb_components = first_not_None(nb_components, domain.dim)
                 assert (is_vector is not True) or (nb_components == domain.dim)
-                return VectorField(domain=domain, name=name, 
+                return VectorField(domain=domain, name=name,
                                         nb_components=nb_components, **kwds)
             else:
                 return ScalarField(domain=domain, name=name, **kwds)
@@ -66,12 +66,12 @@ class FieldContainerI(TaggedObject):
     @property
     def is_scalar(self):
         return (not self.is_tensor)
-    
+
     @abstractmethod
     def field_like(self, name, **kwds):
         """Create a ScalarField like this object, possibly altered."""
         pass
-    
+
     @abstractmethod
     def tmp_like(self, name, **kwds):
         """Create a temporary field like self, possibly altered."""
@@ -103,17 +103,17 @@ class FieldContainerI(TaggedObject):
             The topology state on which to discretize this ScalarField.
         """
         pass
-    
+
     @classmethod
-    def from_sympy_expressions(cls, name, exprs, space_symbols, 
+    def from_sympy_expressions(cls, name, exprs, space_symbols,
                                     scalar_name_prefix=None, scalar_pretty_name_prefix=None,
                                     pretty_name=None,  **kwds):
         """
         Create a field wich has the same shape as exprs, with optional names.
         Expressions should be of kind sympy.Expr and are converted to FieldExpression: this
         means they all have to contain at least one FieldExpression.
-        Note that field.symbol is always a SymbolicField which is a FieldExpression. 
-        FieldExpression make sure boundary conditions match between fields for derivatives 
+        Note that field.symbol is always a SymbolicField which is a FieldExpression.
+        FieldExpression make sure boundary conditions match between fields for derivatives
         and integrations.
         """
         if isinstance(exprs, sm.Expr):
@@ -144,7 +144,7 @@ class FieldContainerI(TaggedObject):
                 sname  = None
                 spname = None
 
-            fields[idx] = cls.from_sympy_expression(expr=exprs[idx], 
+            fields[idx] = cls.from_sympy_expression(expr=exprs[idx],
                             space_symbols=space_symbols,
                             name=sname, pretty_name=spname, **kwds)
         return TensorField.from_field_array(name=name, pretty_name=pretty_name,
@@ -178,15 +178,15 @@ class FieldContainerI(TaggedObject):
 
         # finally return create and return the ScalarField
         return ScalarField(**kwds)
-        
 
-    def gradient(self, name=None, pretty_name=None, 
+
+    def gradient(self, name=None, pretty_name=None,
                        scalar_name_prefix=None, scalar_pretty_name_prefix=None,
-                       directions=None, axis=-1, 
+                       directions=None, axis=-1,
                        space_symbols=None,
                        dtype=None, **kwds):
         """
-        Create a field capable of storing the gradient of self, 
+        Create a field capable of storing the gradient of self,
         possibly altered.
         """
         dim    = self.dim  # dimension of the domain
@@ -195,7 +195,7 @@ class FieldContainerI(TaggedObject):
 
         directions = to_tuple(first_not_None(directions, range(dim)))
         space_symbols = to_tuple(first_not_None(space_symbols, frame.coords))
-        check_instance(directions, tuple, minval=0, maxval=self.dim-1, minsize=1, unique=True) 
+        check_instance(directions, tuple, minval=0, maxval=self.dim-1, minsize=1, unique=True)
         check_instance(axis, int, minval=-ndim, maxval=ndim-1)
         check_instance(space_symbols, tuple, values=SpaceSymbol, size=dim, unique=True)
 
@@ -221,7 +221,7 @@ class FieldContainerI(TaggedObject):
                 i = idx[:axis+1] + idx[axis+2:]
                 d = directions[idx[axis+1]]
                 if self.is_tensor:
-                    exprs[idx] = self[i].symbol(frame.time, 
+                    exprs[idx] = self[i].symbol(frame.time,
                             *space_symbols).diff(space_symbols[d])
                 else:
                     assert i==(), i
@@ -233,7 +233,7 @@ class FieldContainerI(TaggedObject):
                     scalar_pretty_name_prefix=scalar_pretty_name_prefix,
                     dtype=dtype, **kwds)
 
-    def laplacian(self, name=None, pretty_name=None, 
+    def laplacian(self, name=None, pretty_name=None,
                   scalar_name_prefix=None, scalar_pretty_name_prefix=None,
                   dtype=None, **kwds):
         from hysop.symbolic.field import laplacian
@@ -243,7 +243,7 @@ class FieldContainerI(TaggedObject):
         name = first_not_None(name, 'laplacian_{}'.format(self.name))
         pretty_name = first_not_None(pretty_name, u'\u0394{}'.format(
                                      self.pretty_name.decode('utf-8')))
-        
+
         if isinstance(exprs, npw.ndarray):
             if (exprs.size == 1):
                 expr = exprs.item()
@@ -264,11 +264,11 @@ class FieldContainerI(TaggedObject):
                                               dtype=dtype, **kwds)
 
 
-    def div(self, name=None, pretty_name=None, 
+    def div(self, name=None, pretty_name=None,
                   scalar_name_prefix=None, scalar_pretty_name_prefix=None,
                   axis=-1, dtype=None, **kwds):
         """
-        Create a field capable of storing the divergence of self, 
+        Create a field capable of storing the divergence of self,
         on chosen axis.
         """
         from hysop.symbolic.field import div
@@ -278,7 +278,7 @@ class FieldContainerI(TaggedObject):
         name = first_not_None(name, 'div_{}'.format(self.name))
         pretty_name = first_not_None(pretty_name, u'{}\u22c5{}'.format(nabla,
                                      self.pretty_name.decode('utf-8')))
-        
+
         if exprs.size in (0,1):
             expr = npw.asscalar(exprs)
             return self.from_sympy_expression(expr=expr, space_symbols=frame.coords,
@@ -291,23 +291,23 @@ class FieldContainerI(TaggedObject):
                                                scalar_pretty_name_prefix=scalar_pretty_name_prefix,
                                                dtype=dtype, **kwds)
 
-    def curl(self, name=None, pretty_name=None, 
+    def curl(self, name=None, pretty_name=None,
                    scalar_name_prefix=None, scalar_pretty_name_prefix=None,
                    dtype=None, **kwds):
         """
-        Create a field capable of storing the curl of self, 
+        Create a field capable of storing the curl of self,
 
-        Only 2D and 3D fields are supported as the curl brings 
+        Only 2D and 3D fields are supported as the curl brings
         a 1-vector to a 2-vector:
-        
-        - A vector to a pseudoscalar or a pseudoscalar to a vector in 2D 
+
+        - A vector to a pseudoscalar or a pseudoscalar to a vector in 2D
         - A vector to a pseudovector or a pseudovector to a vector in 3D
 
         In 1D the curl is 0, and in 4D the curl would be a 6D 'field'.
         """
         from hysop.symbolic.field import curl
-        
-        
+
+
         if (self.dim==2):
             msg='Can only take curl for a 2D field with one or two components.'
             assert self.nb_components in (1,2), msg
@@ -329,7 +329,7 @@ class FieldContainerI(TaggedObject):
             return self.from_sympy_expressions(
                     exprs=exprs, space_symbols=frame.coords,
                     name=name, pretty_name=pretty_name,
-                    scalar_name_prefix=scalar_name_prefix,            
+                    scalar_name_prefix=scalar_name_prefix,
                     scalar_pretty_name_prefix=scalar_pretty_name_prefix,
                     dtype=dtype, **kwds)
         else:
@@ -341,7 +341,7 @@ class FieldContainerI(TaggedObject):
         """See curl."""
         return self.curl(*args, **kwds)
 
-    
+
     def get_attributes(self, *attrs):
         """
         Return all matching attributes contained in self.fields,
@@ -391,7 +391,7 @@ class FieldContainerI(TaggedObject):
     def get_unique_attribute(self, *attr):
         """
         Try to return the unique attribute common to all contained fields.
-        Raise an AttributeError if a attribute is not unique accross contained 
+        Raise an AttributeError if a attribute is not unique accross contained
         field views.
         """
         if self.has_unique_attribute(*attr):
@@ -458,7 +458,7 @@ class FieldContainerI(TaggedObject):
         return (self is not other)
     def __hash__(self):
         return id(self)
-    
+
 
     def _get_domain(self):
         """Return the physical domain where this field is defined."""
@@ -476,7 +476,7 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
 
     This object handles a dictionnary of discrete fields
     (from 0 to any number).
-    
+
     Each discrete field is uniquely defined by the topology used to
     discretize it.
 
@@ -523,9 +523,9 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
             Underlying data type of this field
         initial_values: numeric value, or tuple of numeric values, optional
             Fields are initialized to specified initial value everywhere in the domain
-            on first discretization. 
+            on first discretization.
             The input values are cast to given dtype.
-            
+
             If None, leaves the memory uninitialized.
             If a single value is given, the whole field is initialized to this value,
              the default being None (ie. no initialization at all).
@@ -539,7 +539,7 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
             Specify that this field is a temporary continuous field.
             Basically a ScalarField that yields a temporary discrete field upon discretization.
 
-            /!\ ** WARNING *******************************************/!\ 
+            /!\ ** WARNING *******************************************/!\
             TemporaryDiscreteFields are allocated during setup using
             temporary work buffers. Those work buffers are only
             available withing the scope of operators thats use
@@ -547,7 +547,7 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
             /!\ ***************************************************** /!\
         kwds: dict
             Base class keyword arguments.
-        
+
         Attributes
         ----------
         boundaries: tuple of numpy.ndarray of BoundaryCondition
@@ -572,44 +572,44 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
             warnings.warn(msg, HysopWarning)
             dtype = HYSOP_BOOL
         dtype = npw.dtype(dtype)
-        
+
         # Name and pretty name
         pretty_name = first_not_None(pretty_name, name)
         if isinstance(pretty_name, unicode):
             pretty_name = pretty_name.encode('utf-8')
         check_instance(pretty_name, str)
-        
+
         # Initial values
         if not isinstance(initial_values,(list,tuple)):
             initial_values = (initial_values, initial_values)
         assert len(initial_values)==2
         initial_values = tuple(initial_values)
         check_instance(initial_values, tuple, size=2)
-        
+
         # Field boundary conditions
-        lboundaries = npw.asarray(first_not_None(lboundaries, 
+        lboundaries = npw.asarray(first_not_None(lboundaries,
             cls.default_boundaries_from_domain(domain.lboundaries)))
-        rboundaries = npw.asarray(first_not_None(rboundaries, 
+        rboundaries = npw.asarray(first_not_None(rboundaries,
             cls.default_boundaries_from_domain(domain.rboundaries)))
-        check_instance(lboundaries, npw.ndarray, values=BoundaryCondition, 
+        check_instance(lboundaries, npw.ndarray, values=BoundaryCondition,
                 ndim=1, size=domain.dim, dtype=object, allow_none=True)
-        check_instance(rboundaries, npw.ndarray, values=BoundaryCondition, 
+        check_instance(rboundaries, npw.ndarray, values=BoundaryCondition,
                 ndim=1, size=domain.dim, dtype=object, allow_none=True)
         assert lboundaries.size == rboundaries.size == domain.dim
         for i,(lb,rb) in enumerate(zip(lboundaries,rboundaries)):
             if (lb==BoundaryCondition.PERIODIC) ^ (rb==BoundaryCondition.PERIODIC):
                 msg='Periodic BoundaryCondition mismatch on axis {}.'.format(i)
                 raise ValueError(msg)
-        check_instance(lboundaries, npw.ndarray, values=BoundaryCondition, 
+        check_instance(lboundaries, npw.ndarray, values=BoundaryCondition,
                 ndim=1, size=domain.dim, dtype=object)
-        check_instance(rboundaries, npw.ndarray, values=BoundaryCondition, 
+        check_instance(rboundaries, npw.ndarray, values=BoundaryCondition,
                 ndim=1, size=domain.dim, dtype=object)
 
         periodic    = BoundaryCondition.PERIODIC
         periodicity = (lboundaries==periodic)
-        
+
         obj = super(ScalarField, cls).__new__(cls, domain=domain,
-                name=name, pretty_name=pretty_name, 
+                name=name, pretty_name=pretty_name,
                 var_name=var_name, latex_name=latex_name,
                 tag_prefix='f', tagged_cls=ScalarField, **kwds)
         obj._dtype  = dtype
@@ -619,7 +619,7 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
         obj._lboundaries = lboundaries
         obj._rboundaries = rboundaries
         obj._periodicity = periodicity
-    
+
         # Symbolic representation of this field
         from hysop.symbolic.field import SymbolicField
         obj._symbol = SymbolicField(field=obj)
@@ -651,7 +651,7 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
                 raise NotImplementedError(msg)
             field_boundaries[i] = fbd
         return field_boundaries
-    
+
     @classmethod
     def __check_vars(cls, obj):
         """Check properties and types."""
@@ -663,17 +663,17 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
         check_instance(obj.nb_components, int, minval=1)
         check_instance(obj.discrete_fields, dict)
         check_instance(obj.initial_values, tuple, size=2)
-        check_instance(obj.lboundaries, npw.ndarray, values=BoundaryCondition, 
+        check_instance(obj.lboundaries, npw.ndarray, values=BoundaryCondition,
                 ndim=1, size=obj.domain.dim, dtype=object)
-        check_instance(obj.rboundaries, npw.ndarray, values=BoundaryCondition, 
+        check_instance(obj.rboundaries, npw.ndarray, values=BoundaryCondition,
                 ndim=1, size=obj.domain.dim, dtype=object)
-        check_instance(obj.periodicity, npw.ndarray, dtype=bool, 
+        check_instance(obj.periodicity, npw.ndarray, dtype=bool,
                 ndim=1, size=obj.domain.dim)
         check_instance(obj.is_tmp, bool)
 
     def field_like(self, name, pretty_name=None,
             latex_name=None, var_name=None,
-            domain=None, dtype=None, is_tmp=None, 
+            domain=None, dtype=None, is_tmp=None,
             lboundaries=None, rboundaries=None,
             initial_values=None, **kwds):
         """Create a ScalarField like this object, possibly altered."""
@@ -686,23 +686,23 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
         initial_values = first_not_None(initial_values, self.initial_values)
         return ScalarField(name=name, pretty_name=pretty_name,
                 var_name=var_name, latex_name=latex_name,
-                domain=domain, dtype=dtype, is_tmp=is_tmp, 
-                lboundaries=lboundaries, rboundaries=rboundaries, 
+                domain=domain, dtype=dtype, is_tmp=is_tmp,
+                lboundaries=lboundaries, rboundaries=rboundaries,
                 initial_values=initial_values, **kwds)
 
     def tmp_like(self, name, **kwds):
         """Create a TemporaryField like self, possibly altered."""
         return self.field_like(name=name, is_tmp=True, **kwds)
-    
+
     def short_description(self):
         """Short description of this field."""
         s = '{}[pname={}, dim={}, dtype={}, bc=[{}], iv={}]'
-        s = s.format(self.full_tag, self.name, self.dim, 
-                     self.dtype, 
+        s = s.format(self.full_tag, self.name, self.dim,
+                     self.dtype,
                      self.format_boundaries(),
                      self.initial_values)
         return s
-    
+
     def format_boundaries(self):
         from hysop.constants import format_boundaries as fb
         return fb(*self.boundaries)
@@ -723,15 +723,15 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
           *initial values: {}
           *topology tags:  [{}]
         ''').format(self.full_tag,
-                self.name, self.pretty_name, 
+                self.name, self.pretty_name,
                 self.var_name, self.latex_name,
                 self.dim, self.dtype,
-                self.lboundaries.tolist(), self.rboundaries.tolist(), 
+                self.lboundaries.tolist(), self.rboundaries.tolist(),
                 self.initial_values,
                 ','.join([k.full_tag for k in self.discrete_fields.keys()]))
         return s[1:]
-                
-    
+
+
     @debug
     def discretize(self, topology, topology_state=None):
         """
@@ -797,7 +797,7 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
         return self._is_tmp
     def _get_mem_tag(self):
         return self._mem_tag
-    
+
     dtype = property(_get_dtype)
     initial_values = property(_get_initial_values)
     discrete_fields = property(_get_discrete_fields)
@@ -807,15 +807,15 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
     periodicity = property(_get_periodicity)
     is_tmp = property(_get_is_tmp)
     mem_tag = property(_get_mem_tag)
-    
+
     @property
     def is_tensor(self):
         return False
-    
+
     @property
     def fields(self):
         return (self,)
-    
+
     @property
     def nb_components(self):
         return 1
@@ -832,13 +832,13 @@ class ScalarField(NamedScalarContainerI, FieldContainerI):
 
 class TensorField(NamedTensorContainerI, FieldContainerI):
     """
-    A continuous tensor field is a collection of scalar fields 
+    A continuous tensor field is a collection of scalar fields
     defined on a physical domain, organized as a multi-dimensional array..
 
     This object handles a numpy.ndarray of continuous scalar fields,
-    which may have different attributes (different data types for 
+    which may have different attributes (different data types for
     example). It is mainly a view on scalar fields with the FieldContainerI interface.
-    
+
     A tensor field garanties that each different field objects have
     unique names and pretty names within the tensor field. A given
     continuous scalar field may appear in at multiple indices (to allow
@@ -846,13 +846,13 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
     (to allow upper triangular matrices for example). Is also garanties
     that all fields shares the same domain.
     """
-    
+
     @property
     def is_tensor(self):
         return True
-    
-    def __new__(cls, domain, name, shape, 
-                    pretty_name=None, 
+
+    def __new__(cls, domain, name, shape,
+                    pretty_name=None,
                     name_formatter=None,
                     pretty_name_formatter=None,
                     skip_field=None,
@@ -868,28 +868,28 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
         check_instance(shape, tuple, values=int)
         if (len(shape)==1) and not issubclass(cls, VectorField):
             obj = VectorField(domain=domain, shape=shape,
-                                name=name, 
+                                name=name,
                                 name_formatter=name_formatter,
-                                pretty_name=pretty_name, 
+                                pretty_name=pretty_name,
                                 pretty_name_formatter=pretty_name_formatter,
                                 skip_field=skip_field, make_field=make_field,
                                 fields=fields, base_kwds=base_kwds, **kwds)
             return obj
-        
+
         name_formatter = first_not_None(name_formatter, cls.default_name_formatter)
-        pretty_name_formatter = first_not_None(pretty_name_formatter, 
+        pretty_name_formatter = first_not_None(pretty_name_formatter,
                                                cls.default_pretty_name_formatter)
         skip_field = first_not_None(skip_field, cls.default_skip_field)
         make_field = first_not_None(make_field, cls.default_make_field)
         base_kwds = first_not_None(base_kwds, {})
-        
+
         check_instance(domain, Domain)
 
         if npw.prod(shape)<=0:
             msg='Invalid shape for a tensor-like field, got {}.'
             msg=msg.format(shape)
             raise ValueError(msg)
-        
+
         if (fields is None):
             fields = ()
             for idx in npw.ndindex(*shape):
@@ -898,29 +898,29 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
                 else:
                     fname  = name_formatter(basename=name, idx=idx)
                     pfname = pretty_name_formatter(basename=pretty_name, idx=idx)
-                    field  = make_field(idx, domain=domain, name=fname, pretty_name=pfname, 
+                    field  = make_field(idx, domain=domain, name=fname, pretty_name=pfname,
                                         **kwds)
                 fields += (field,)
             cls._check_fields(*fields)
             fields = npw.asarray(fields, dtype=object).reshape(shape)
         else:
             assert (not kwds)
-        
+
         check_instance(fields, npw.ndarray, dtype=object)
         assert npw.array_equal(fields.shape, shape)
 
-        obj = super(TensorField, cls).__new__(cls, domain=domain, 
-                name=name, pretty_name=pretty_name, 
-                tag_prefix='tf', tagged_cls=TensorField, 
+        obj = super(TensorField, cls).__new__(cls, domain=domain,
+                name=name, pretty_name=pretty_name,
+                tag_prefix='tf', tagged_cls=TensorField,
                 contained_objects=fields, **base_kwds)
         obj._fields = fields
         obj._name_formatter = name_formatter
         obj._pretty_name_formatter = pretty_name_formatter
         obj._skip_field = skip_field
-        
+
         from hysop.symbolic.field import SymbolicFieldTensor
         obj._symbol = SymbolicFieldTensor(field=obj)
-        
+
         obj._check_domains(domain)
         obj._check_names()
         return obj
@@ -929,7 +929,7 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
         from hysop.fields.discrete_field import DiscreteTensorField
         dfields = npw.empty(shape=self.shape, dtype=object)
         for (idx, field) in self.nd_iter():
-            dfields[idx] = field.discretize(topology=topology, 
+            dfields[idx] = field.discretize(topology=topology,
                                             topology_state=topology_state)
         return DiscreteTensorField(field=self, dfields=dfields)
 
@@ -953,7 +953,7 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
 
         return Field(domain=domain, name=name, shape=shape, pretty_name=pretty_name,
                     fields=fields, **kwds)
-    
+
     @classmethod
     def from_field_array(cls, name, fields, pretty_name=None, **kwds):
         """Create a TensorField from numpy.ndarray of fields."""
@@ -962,7 +962,7 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
         check_instance(pretty_name, (str, unicode), allow_none=True)
         check_instance(fields, npw.ndarray, dtype=object, values=ScalarField)
         shape = fields.shape
-        
+
         _fields = tuple(fields.ravel().tolist())
         cls._check_fields(*_fields)
 
@@ -971,7 +971,7 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
 
         return Field(domain=domain, name=name, pretty_name=pretty_name,
                 shape=shape, fields=fields, **kwds)
-        
+
     @classmethod
     def _check_fields(cls, *fields):
         """Check that at least one field is specified."""
@@ -980,9 +980,10 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
             msg='Tensor field {} should at least contain a valid ScalarField.'
             msg=msg.format(name)
             raise ValueError(msg)
-    
+
     @classmethod
     def default_name_formatter(cls, basename, idx):
+
         assert len(basename)>0
         if basename[-1] in '0123456789':
             sep = '_'
@@ -1001,11 +1002,11 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
     @classmethod
     def default_make_field(cls, idx, **kwds):
         return ScalarField(**kwds)
-    
+
     @classmethod
     def default_skip_field(cls, idx):
         return False
-    
+
     def _check_domains(self, domain):
         """Check that fields share a unique domain."""
         for field in self:
@@ -1013,7 +1014,7 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
             if (field.domain.domain is not domain):
                 msg='Domain mismatch for field {}.'.format(field.name)
                 raise ValueError(msg)
-    
+
     def _check_names(self):
         """Check that fields names are unique."""
         names = {}
@@ -1031,7 +1032,7 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
                 raise ValueError(msg)
             names[name]  = field
             pnames[name] = field
-    
+
     @property
     def fields(self):
         """Return all unique scalar fields contained in this field-like interface."""
@@ -1045,18 +1046,18 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
             if (field in fields) and (field not in unique_fields):
                 unique_fields += (field,)
         return unique_fields
-    
+
     @property
     def nb_components(self):
         return len(self.fields)
-    
+
     def short_description(self):
         """Short description of this tensor field."""
         s = '{}[name={}, pname={}, dim={}, shape={}]'
         s = s.format(self.full_tag, self.name, self.pretty_name, self.dim,
                         self.shape)
         return s
-    
+
     def long_description(self):
         """Long description of this tensor field as a string."""
         s=textwrap.dedent(
@@ -1073,8 +1074,8 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
         s+='      '+'\n      '.join(str(self.symbol).split('\n'))
         return s
 
-    def field_like(self, name, pretty_name=None, 
-                        shape=None, nb_components=None, 
+    def field_like(self, name, pretty_name=None,
+                        shape=None, nb_components=None,
                         fn='field_like', **kwds):
         """Create a TensorField like this object, possibly altered."""
         if (shape is None) and (nb_components is not None):
@@ -1082,13 +1083,13 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
         del nb_components
         shape = first_not_None(shape, self.shape)
         nb_components = npw.prod(shape, dtype=npw.int64)
-        
+
         pretty_name = first_not_None(pretty_name, name)
         check_instance(name, str)
         check_instance(pretty_name, (str,unicode))
         if not isinstance(pretty_name, str):
             pretty_name = pretty_name.encode('utf-8')
-        
+
         if (nb_components == 1):
             return getattr(self.fields[0], fn)(name=name, pretty_name=pretty_name, **kwds)
         else:
@@ -1105,7 +1106,7 @@ class TensorField(NamedTensorContainerI, FieldContainerI):
                     pfname = self._pretty_name_formatter(basename=pretty_name, idx=idx)
                     fields[idx] = getattr(field, fn)(name=fname, pretty_name=pfname, **kwds)
             return self.from_field_array(name=name, pretty_name=pretty_name, fields=fields)
-    
+
     def tmp_like(self, name, **kwds):
         """Create a temporary field like self, possibly altered."""
         return self.field_like(name=name, fn='tmp_like', **kwds)
@@ -1130,11 +1131,23 @@ class VectorField(TensorField):
         check_instance(shape, tuple, values=int, size=1)
         if shape[0]==1:
             return ScalarField(domain=domain, name=name, **kwds)
-        
-        obj = super(VectorField, cls).__new__(cls, domain=domain, name=name, 
+        obj = super(VectorField, cls).__new__(cls, domain=domain, name=name,
                                                    shape=shape, **kwds)
         return obj
 
+    @classmethod
+    def default_name_formatter(cls, basename, idx):
+        assert len(basename)>0
+        if basename[-1] in '0123456789':
+            sep = '_'
+        else:
+            sep = ''
+        if len(idx)==1:
+            name = basename + sep + '_'.join(DirectionLabels[i] for i in idx)
+        else:
+            name = basename + sep + '_'.join(str(i) for i in idx)
+        return name
+
 
 Field = FieldContainerI
 """A Field is just a alias of FieldContainerI"""
diff --git a/hysop/operator/hdf_io.py b/hysop/operator/hdf_io.py
index cadbcb95938c2ef3fb2012d6f7f5671916b7dfef..c913ae8c702b10b6786ca03a88feb7dd195b9590 100755
--- a/hysop/operator/hdf_io.py
+++ b/hysop/operator/hdf_io.py
@@ -12,7 +12,7 @@ from abc import ABCMeta, abstractmethod
 from hysop import __H5PY_PARALLEL_COMPRESSION_ENABLED__, vprint
 from hysop.deps import h5py, sys
 from hysop.core.graph.graph import discretized
-from hysop.constants import DirectionLabels, HYSOP_REAL, Backend, TranspositionState
+from hysop.constants import HYSOP_REAL, Backend, TranspositionState
 from hysop.tools.decorators import debug
 from hysop.tools.types import check_instance, first_not_None
 from hysop.tools.numpywrappers import npw
@@ -196,7 +196,7 @@ class HDF_IO(ComputationalGraphOperator):
         self._global_grid_resolution = refmesh.grid_resolution
         # Local resolution for hdf5 output
         self._local_grid_resolution = refmesh.compute_resolution
-        
+
         assert self.io_params.io_leader<topo.cart_comm.size
         self._all_local_grid_resolution = topo.cart_comm.gather(
             self._local_grid_resolution, root=self.io_params.io_leader)
@@ -224,7 +224,10 @@ class HDF_IO(ComputationalGraphOperator):
             # Get field names and initialize dataset dict.
             for df in self.discrete_fields:
                 for d in xrange(df.nb_components):
-                    name = name_prefix + df.name + '_' + DirectionLabels[d] + name_postfix
+                    if df.nb_components == 1:
+                        name = name_prefix + df.name + name_postfix
+                    else:
+                        name = name_prefix + df.name + "_{}".format(d) + name_postfix
                     self.dataset[name] = df.data[d]
                     var_names[df.field] = name
             self.var_names = var_names
@@ -516,7 +519,7 @@ class HDF_Writer(HDF_IO):
             raise RuntimeError(msg)
         self._xdmf_data_files.append((self._count, simu.t()))
         self._last_written_time = simu.t()
-        
+
         self._hdf_file.close()
         self.updateXMFFile()
 
@@ -525,18 +528,18 @@ class HDF_Writer(HDF_IO):
             op_name         = self.name
             actual_filepath = self.io_params.filepath
             disk_filepath   = self.io_params.disk_filepath
-            xmf_file   = self._xmf_file.name 
+            xmf_file   = self._xmf_file.name
             hdf_file   = filename
             hdf_is_tmp = self.io_params.dump_is_temporary
             iteration = self._count
             time = self._last_written_time
-            
+
             vprint('>Executing postprocessing script: {}'.format(postprocess_cmd))
-            
+
             # execute command OP_NAME  ACTUAL_FILEPATH  DISK_FILEPATH  XMF_FILE  HDF5_FILE  IS_TMP
-            command = [str(postprocess_cmd), 
+            command = [str(postprocess_cmd),
                        str(op_name), str(actual_filepath), str(disk_filepath),
-                       str(xmf_file), str(hdf_file), 
+                       str(xmf_file), str(hdf_file),
                        '1' if hdf_is_tmp else '0',
                        str(iteration), str(time)]
             try:
@@ -557,7 +560,7 @@ class HDF_Writer(HDF_IO):
                     print(msg)
                     print
                     raise
-        
+
         if self.io_params.dump_is_temporary:
             vprint('>Deleting HDF5 data {}...'.format(filename))
             os.remove(filename)