From f391f89029a50ed5c67a292831bd7e5360fc81c1 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Mon, 10 Dec 2018 22:14:11 +0100
Subject: [PATCH] ext force

---
 examples/sediment_deposit/sediment_deposit.py |  30 +--
 .../device/opencl/operator/external_force.py  | 103 ++++++++
 hysop/fields/continuous_field.py              | 174 +------------
 hysop/operator/external_force.py              | 233 ++++++++++++------
 hysop/operators.py                            |   2 +-
 hysop/tools/interface.py                      | 176 +++++++++++++
 6 files changed, 448 insertions(+), 270 deletions(-)
 create mode 100644 hysop/backend/device/opencl/operator/external_force.py
 create mode 100644 hysop/tools/interface.py

diff --git a/examples/sediment_deposit/sediment_deposit.py b/examples/sediment_deposit/sediment_deposit.py
index ce1eb0855..e5f827dc5 100644
--- a/examples/sediment_deposit/sediment_deposit.py
+++ b/examples/sediment_deposit/sediment_deposit.py
@@ -121,7 +121,7 @@ def compute(args):
                                 Enstrophy, MinMaxFieldStatistics, StrangSplitting, \
                                 ParameterPlotter, Integrate, HDF_Writer,           \
                                 CustomSymbolicOperator, DirectionalSymbolic,       \
-                                SymbolicSpectralExternalForce
+                                SpectralExternalForce, SymbolicExternalForce
 
     from hysop.methods import SpaceDiscretization, Remesh, TimeIntegrator, \
                               ComputeGranularity, Interpolation
@@ -231,17 +231,6 @@ def compute(args):
         msg='Unsupported dimension {}.'.format(dim)
         raise RuntimeError(msg)
     
-    #> Diffusion of vorticity and S
-    diffuse_S = Diffusion(implementation=impl,
-             enforce_implementation=enforce_implementation,
-             name='diffuse_S',
-             pretty_name='diffS',
-             nu = nu_S,
-             Fin = S,
-             variables = {S: npts},
-             dt=dt, 
-             **extra_op_kwds)
-    
     splitting = StrangSplitting(splitting_dim=dim, 
                     order=args.strang_order)
     splitting.push_operators(advec, stretch)
@@ -255,14 +244,14 @@ def compute(args):
                             enforce_implementation=enforce_implementation,
                             **extra_op_kwds)
 
-    #> External force rot(-rho*g) = rot(-S)
-    external_force = SymbolicExternalForce(vorticity=Ws, 
-                                   compute_statistics=True,
-                                   Fext=[0,-Ss],
-                                   dt=dt,
+    #> External force rot(-rho*g) = rot(-(1+S)) = rot(-S)
+    g = 1.0
+    Fext = SymbolicExternalForce(name='S', Fext=(0,-g*Ss),
+                                   diffusion = {S: nu_S})
+    external_force = SpectralExternalForce(vorticity=Ws, dt=dt,
+                                   Fext=Fext, Finf=True,
                                    implementation=impl,
-                                   variables={vorti: npts,
-                                              S: npts},
+                                   variables={vorti: npts, S: npts},
                                    **extra_op_kwds)
     
     #> Operator to compute the infinite norm of the velocity
@@ -306,7 +295,7 @@ def compute(args):
                                                  name='dt_lcfl', pretty_name='LCFL')
     dt_force = adapt_dt.push_cst_criteria(cst=1,
                                         Finf=external_force.Finf,
-                                        name='dt_force', pretty_name='CST')
+                                        name='dt_force', pretty_name='FEXT')
 
     
     ## Create the problem we want to solve and insert our 
@@ -327,7 +316,6 @@ def compute(args):
     problem.insert(poisson, 
                    min_max_U, min_max_W, adapt_dt,
                    splitting, 
-                   #diffuse_S,
                    external_force,
                    dump_fields,
                    compute_mean_fields)
diff --git a/hysop/backend/device/opencl/operator/external_force.py b/hysop/backend/device/opencl/operator/external_force.py
new file mode 100644
index 000000000..d5e10fb66
--- /dev/null
+++ b/hysop/backend/device/opencl/operator/external_force.py
@@ -0,0 +1,103 @@
+
+import primefac, functools
+from hysop.tools.numerics import float_to_complex_dtype
+from hysop.tools.numpywrappers import npw
+from hysop.tools.decorators  import debug
+from hysop.tools.warning import HysopWarning
+from hysop.operator.base.curl import SpectralExternalForceOperatorBase
+from hysop.backend.device.opencl.opencl_symbolic import OpenClSymbolic
+from hysop.core.graph.graph import op_apply
+from hysop.core.memory.memory_request import MemoryRequest
+from hysop.backend.device.opencl.opencl_fft import OpenClFFT
+from hysop.backend.device.codegen.base.variables import dtype_to_ctype
+from hysop.symbolic import local_indices_symbols
+from hysop.symbolic.relational import Assignment
+from hysop.symbolic.complex import ComplexMul
+from hysop.symbolic.misc import Select
+
+class OpenClSpectralExternalForce(SpectralExternalForceOperatorBase, OpenClSymbolic):
+    """
+    Operator to compute the curl of a symbolic expression.
+    """
+    
+    @debug
+    def __init__(self, **kwds):
+        """
+        What is computed in 2D:
+            FEXT_X = fft(Fext_X)
+            tmp    = dFEXT/dy
+            FEXT_Y = fft(Fext_Y)
+            tmp   -= dFEXT/dx
+            Fmin[0], Fmax[0], Finf[0] = min(tmp), max(tmp), max(abs(tmp))
+            Wz  += dt*tmp
+        What is computed in 3D:
+            FEXT_X = fft(Fext_X)
+            FEXT_Y = fft(Fext_Y)
+            FEXT_Z = fft(Fext_Z)
+            
+            tmp = 
+            Fmin[0], Fmax[0], Finf[0] = min(tmp), max(tmp), max(abs(tmp))
+            Wx += dt*tmp
+
+            tmp = 
+            Fmin[1], Fmax[1], Finf[1] = min(tmp), max(tmp), max(abs(tmp))
+            Wy += dt*tmp
+
+            tmp = 
+            Fmin[1], Fmax[1], Finf[1] = min(tmp), max(tmp), max(abs(tmp))
+            Wz += dt*tmp
+        """
+        super(OpenClSpectralExternalForce, self).__init__(**kwds)
+
+        assert (len(self.forward_transforms) % 2 == 0)
+        N = len(self.forward_transforms)//2
+        assert len(self.K)==2*N
+
+        kernel_names = ()
+        for (i,(Ft,(tg,Ki))) in enumerate(zip(self.forward_transforms, self.K)):
+            Fhs = Ft.output_symbolic_array('F{}_hat'.format(i))
+            
+            kname = 'filter_curl_{}d_{}'.format(Fhs.dim, i)
+            kernel_names += (kname,)
+            
+            is_complex = Ki.is_complex
+            Ki = tg._indexed_wave_numbers[Ki]
+
+            if is_complex:
+                expr = ComplexMul(Ki, Fhs)
+            else:
+                expr = Ki*Fhs
+
+            if (i<N):
+                expr = Assignment(Fhs, +expr)
+            else:
+                expr = Assignment(Fhs, -expr)
+
+            self.require_symbolic_kernel(kname, expr)
+
+        self._kernel_names = kernel_names
+    
+    @debug
+    def setup(self, work):
+        super(OpenClSpectralExternalForce, self).setup(work)
+        curl_filters = ()
+        for kname in self._kernel_names:
+            kernel, _ = self.symbolic_kernels[kname]
+            kernel = functools.partial(kernel, queue=self.cl_env.default_queue)
+            curl_filters += (kernel,)
+        self.curl_filters = curl_filters
+        self.exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True)
+        assert len(self.forward_transforms)==len(self.backward_transforms)==len(curl_filters)
+
+    @op_apply
+    def apply(self, **kwds):
+        """Solve the ExternalForce equation."""
+        super(OpenClSpectralExternalForce, self).apply(**kwds)
+        for (Ft,filter_curl,Bt) in zip(self.forward_transforms,
+                                       self.curl_filters,
+                                       self.backward_transforms):
+            evt = Ft()
+            evt = filter_curl()
+            evt = Bt()
+        if (self.exchange_ghosts is not None):
+            evt = self.exchange_ghosts()
diff --git a/hysop/fields/continuous_field.py b/hysop/fields/continuous_field.py
index 2815c76a2..d60f6776f 100644
--- a/hysop/fields/continuous_field.py
+++ b/hysop/fields/continuous_field.py
@@ -22,178 +22,8 @@ from hysop.topology.topology import Topology, TopologyState
 from hysop.tools.sympy_utils import nabla, partial, subscript, subscripts, \
                                     exponent, exponents, xsymbol
 from hysop.symbolic import SpaceSymbol
-
-class SymbolContainerI(object):
-    __metaclass__ = ABCMeta
-    
-    def _get_symbol(self):
-        """
-        Return a Symbol that can be used to compute symbolic expressions
-        referring to this continuous field.
-        """
-        assert hasattr(self, '_symbol'), 'Symbol has not been defined.'
-        return self._symbol
-    
-    symbol = property(_get_symbol)
-    s = property(_get_symbol)
-
-
-class NamedObjectI(object):
-    __metaclass__ = ABCMeta
-    
-    def __new__(cls, name, pretty_name=None, latex_name=None, var_name=None, **kwds):
-        """
-        Create an abstract named object that contains a symbolic value.
-        name : string
-            A name for the field.
-        pretty_name: string or unicode, optional.
-            A pretty name used for display whenever possible (unicode supported).
-            Defaults to name.
-        kwds: dict
-            Keywords arguments for base class.
-        """
-        
-        obj = super(NamedObjectI, cls).__new__(cls, **kwds)
-        obj.rename(name=name, pretty_name=pretty_name, 
-                    latex_name=latex_name, var_name=var_name)
-        return obj
-    
-    def rename(self, name, pretty_name=None, latex_name=None):
-        """Change the names of this object."""
-        check_instance(name, str)
-        check_instance(pretty_name, (str,unicode), allow_none=True)
-        check_instance(latex_name, str, allow_none=True)
-        if isinstance(pretty_name, unicode):
-            pretty_name = pretty_name.encode('utf-8')
-        check_instance(pretty_name, str)
-        
-        pretty_name = first_not_None(pretty_name, name)
-        latex_name  = first_not_None(latex_name, name)
-
-        self._name   = name
-        self._pretty_name = pretty_name
-        self._latex_name = latex_name
-        
-    def _get_name(self):
-        """Return the name of this field."""
-        return self._name
-    def _get_pretty_name(self):
-        """Return the pretty name of this field."""
-        return self._pretty_name
-    def _get_latex_name(self):
-        """Return the latex name of this field."""
-        return self._latex_name
-    
-    def __str__(self):
-        return self.long_description()
-    
-    @abstractmethod
-    def short_description(self):
-        """Short description of this field as a string."""
-        pass
-
-    @abstractmethod
-    def long_description(self):
-        """Long description of this field as a string."""
-        pass
-
-    name = property(_get_name)
-    pretty_name = property(_get_pretty_name)
-    latex_name = property(_get_latex_name)
-
-
-class NamedScalarContainerI(NamedObjectI, SymbolContainerI):
-
-    @property
-    def ndim(self):
-        """Number of dimensions of this this tensor."""
-        return 0
-    
-    def _get_var_name(self):
-        """Return the variable name of this field."""
-        return self._var_name
-    
-    def rename(self, name, pretty_name=None, 
-                latex_name=None, var_name=None):
-        """Change the names of this object."""
-        super(NamedScalarContainerI, self).rename(name=name, 
-                pretty_name=pretty_name, latex_name=latex_name)
-        self.check_and_set_varname(first_not_None(var_name, self._name))
-    
-    def check_and_set_varname(self, var_name):
-        check_instance(var_name, str, allow_none=True)
-
-        msg='Invalid variable name {}.'.format(var_name)
-        if var_name[0] in tuple(str(x) for x in range(10)):
-            raise RuntimeError(msg)
-        for c in '/*+-=|&()[]{}-!?:;,\'"#$^%<>@':
-            if c in var_name:
-                raise RuntimeError(msg)
-        self._var_name = var_name
-    
-    var_name = property(_get_var_name)
-
-
-class NamedTensorContainerI(NamedObjectI, SymbolContainerI):
-    @debug
-    def __new__(cls, contained_objects, **kwds):
-        check_instance(contained_objects, npw.ndarray)
-        obj = super(NamedTensorContainerI, cls).__new__(cls, **kwds)
-        obj._contained_objects = contained_objects
-        return obj
-    
-    def rename(self, name, pretty_name=None, 
-                latex_name=None, var_name=None):
-        """Change the names of this object."""
-        assert (var_name is None), 'Tensor do not have variable names.'
-        super(NamedTensorContainerI, self).rename(name=name, 
-                pretty_name=pretty_name, latex_name=latex_name)
-    
-    @property
-    def size(self):
-        """Full size of this container as if it was a 1D tensor."""
-        return self._contained_objects.size
-
-    @property
-    def shape(self):
-        """Shape of this tensor."""
-        return self._contained_objects.shape
-    
-    @property
-    def ndim(self):
-        """Number of dimensions of this this tensor."""
-        return self._contained_objects.ndim
-
-    def new_empty_array(self, dtype=object):
-        """Return a new empty array of the same shape as self."""
-        if (dtype is object):
-            array = npw.empty(shape=self.shape, dtype=dtype)
-            array[...] = None
-        else:
-            array = npw.zeros(shape=self.shape, dtype=dtype)
-        return array
-
-    def iter_fields(self):
-        """Return an iterator on unique scalar object along with 1d index."""
-        for (i,obj) in enumerate(self._contained_objects.ravel()):
-            yield (i,obj)
-
-    def nd_iter(self):
-        """Return an nd-indexed iterator of contained objects."""
-        for idx in npw.ndindex(*self._contained_objects.shape):
-            yield (idx, self._contained_objects[idx])
-    
-    def __iter__(self):
-        """Return an iterator on unique scalar objects."""
-        return self._contained_objects.ravel().__iter__()
-
-    def __contains__(self, obj):
-        """Check if a scalar object is contained in self."""
-        return obj in self._contained_objects
-    
-    @abstractmethod
-    def __getitem__(self, slc):
-        pass
+from hysop.tools.interface import NamedObjectI, SymbolContainerI, \
+        NamedScalarContainerI, NamedTensorContainerI
 
 
 class FieldContainerI(TaggedObject):
diff --git a/hysop/operator/external_force.py b/hysop/operator/external_force.py
index e94cad86b..564a113a5 100644
--- a/hysop/operator/external_force.py
+++ b/hysop/operator/external_force.py
@@ -1,66 +1,160 @@
 
-from hysop.b
+from abc import ABCMeta, abstractmethod
+from hysop.fields
+from hysop.symbolic.field import AppliedSymbolicField, SymbolicDiscreteField
+from hysop.symbolic.parameter import SymbolicScalarParameter
+from hysop.parameters.scalar_parameter import ScalarParameter
+from hysop.operator.min_max import MinMaxFieldStatisticsBase
 
-class SymbolicSpectralExternalForce(SpectralCurl):
-    """Generate an operator to compute the curl of a symbolic expression."""
+class ExternalForce(NamedObjectI):
+    __metaclass__ = ABCMeta
+
+    def __init__(self, name, dim, Fext, **kwds):
+        super(ExternalForce, self).__init__(name=name, **kwds)
+
+        check_instance(dim, int)
+        self._dim  = dim
+        self._Fext = Fext
+
+        for f in self.input_fields:
+            msg='Expected field dimension to be {} but got {} from field {}.'
+            msg=msg.format(dim, field.dim, field.short_description())
+            assert f.dim == dim, msg
+
+    
+    @property
+    def name(self):
+        return self._diffusion
+    @property
+    def pretty_name(self):
+        return self._diffusion
+    @property
+    def Fext(self):
+        return self._Fext
+    @property
+    def dim(self):
+        return self._dim
+    @property
+    def diffusion(self):
+        return self._diffusion
+
+    @abstractmethod
+    def input_fields(self):
+        pass
+    @abstractmethod
+    def output_fields(self):
+        pass
+    @abstractmethod
+    def input_params(self):
+        pass
+    @abstractmethod
+    def output_params(self):
+        pass
+
+
+class SymbolicExternalForce(object):
+    def __init__(self, name, Fext, diffusion=None, **kwds):
+        """
+        Specify an external force as a tuple of symbolic expressions.
+
+        2D ExternalForce example:  
+            1) Fext = -rho*g*e_y where rho is a field and g a constant
+                Fext = (0, -rho.s*g)
+            2) Fext = (Rs*S+C)*e_y
+                Fext = (0, -Rs*S.s+C.s)
+        """
+        check_instance(Fext, tuple)
+        dim = len(Fext)
+
+        Fext = list(Fext)
+        for i,e in enumerate(Fext):
+            if isinstance(e, (type(None), int, long, float)):
+                Fext[i] = 0 # curl(const) = 0
+            msg='Expression "{}" contains a SymbolicField, did you forget to apply it ?'
+            msg=msg.format(e)
+            assert not e.atoms(SymbolicField), msg
+        Fext = tuple(Fext)
+
+        super(SymbolicExternalForce, self).__init__(name=name, dim=dim, Fext=Fext, **kwds)
+        
+        diffusion = first_not_None(diffusion, {})
+        for (k,v) in diffusion.iteritems():
+            assert k in self.input_fields
+            assert isinstance(v, ScalarParameter)
+        self._diffusion = diffusion
+
+    def input_fields(self):
+        return set(map(lambda f: f.field, self._extract_objects(AppliedSymbolicField)))
+    def output_fields(self):
+        return set(self._diffusion.keys())
+    def input_params(self):
+        p0 = set(map(lambda p: p.parameter, self._extract_objects(SymbolicScalarParameter)))
+        p1 = set(self._diffusion.values())
+        return p0.union(p1)
+    def output_params(self):
+        return set()
+
+    def _extract_objects(self, obj_type):
+        objs = ()
+        for e in self.Fext:
+            objs += e.atoms(obj_type)
+        return objs
+
+
+class SpectralExternalForce(ComputationalGraphNodeFrontend):
+    """
+    Generate an operator to compute the curl of a symbolic expression.
+    """
 
     @classmethod
     def implementations(cls):
-        raise NotImplementedError
-
+        from hysop.backend.device.opencl.operator.external_force import OpenClSpectralExternalForce
+        __implementations = {
+                Implementation.OPENCL: OpenClSpectralExternalForce,
+        }
+        return __implementations
+    
     @classmethod
     def default_implementation(cls):
-        raise NotImplementedError
+        return Implementation.OPENCL
 
     @debug
-    def __init__(self, vorticity, dt, variables,
-                    Fext=None, Fmin=None, Fmax=None, Finf=None,
+    def __init__(self, vorticity, Fext, dt, variables,
+                    Fmin=None, Fmax=None, Finf=None,
                     implementation=None, base_kwds=None, **kwds):
         """
-        Create an operator that computes the curl of an input field Fin.
+        Create an operator that computes the curl of a given input field Fext or a set of symbolic expressions.
 
-        Given Fin, a 2D ScalarField or VectorField or a 3D VectorField, compute Fout = curl(Fin).
-        
         Only the following configurations are supported:
                      dim   nb_components  |   dim   nb_components
-        Fin:          2          1        |    3          3
-
-        What is computed in 2D:
-            FEXT_X = fft(Fext_X)
-            tmp    = dFEXT/dy
-            FEXT_Y = fft(Fext_Y)
-            tmp   -= dFEXT/dx
-            Fmin[0], Fmax[0], Finf[0] = min(tmp), max(tmp), max(abs(tmp))
-            Wz  += dt*tmp
-        What is computed in 3D:
-            FEXT_X = fft(Fext_X)
-            FEXT_Y = fft(Fext_Y)
-            FEXT_Z = fft(Fext_Z)
-            
-            tmp = 
-            Fmin[0], Fmax[0], Finf[0] = min(tmp), max(tmp), max(abs(tmp))
-            Wx += dt*tmp
-
-            tmp = 
-            Fmin[1], Fmax[1], Finf[1] = min(tmp), max(tmp), max(abs(tmp))
-            Wy += dt*tmp
-
-            tmp = 
-            Fmin[1], Fmax[1], Finf[1] = min(tmp), max(tmp), max(abs(tmp))
-            Wz += dt*tmp
+        vorticity:    2          1        |    3          3
+
+        What is computed:
+            tmp  = curl(Fext) by using a spectral backend
+            Fmin = min(tmp)
+            Fmax = max(tmp)
+            Finf = max(abs(Fmin), abs(Fmax))
+            W   += dt*tmp
         
-        where Fext can be:
-            - an external tensor field of (dim 2 and nb_components=2) or (dim 3 and nb_components 3) 
-            - an array-like of 2 or 3 symbolic expressions to compute Fext components
+        where Fext is computed from user given ExternalForce.
 
         Parameters
         ----------
         vorticity: hysop.field.continuous_field.Field
             Continuous field as input ScalarField or VectorField.
             All contained field have to live on the same domain.
-        Fext: hysop.field.continuous_field.Field or array-like of expressions
-            Continuous field as input VectorField.
-            All contained field have to live on the same domain.
+        Fext: hysop.operator.external_force.ExternalForce
+            Expression of the external force.
+        F...: TensorParameter or boolean, optional
+            TensorParameters should match the shape of tmp (see Notes).
+            If set to True, the TensorParameter will be generated automatically.
+            Autogenerated TensorParameters that are not required by the user are set to be quiet.
+        pbasename: str, optional
+            Parameters basename for created parameters.
+            Defaults to 'curl_{}'.format(Fext.name).
+        ppbasename: str, optional
+            Parameters pretty basename for created parameters.
+            Defaults to '|{} x {}|'.format(nabla, Fext.pretty_name).
         variables: dict
             dictionary of fields as keys and topologies as values.
         implementation: Implementation, optional, defaults to None
@@ -68,11 +162,27 @@ class SymbolicSpectralExternalForce(SpectralCurl):
             If None, implementation will be set to default_implementation().
         kwds: dict, optional
             Extra parameters passed towards base class (MultiSpaceDerivatives).
+
+        Notes
+        -----
+        If Fext.dim == 2, Fext
         """
         base_kwds = first_not_None(base_kwds, {})
-        check_instance(Fin, Field)
-        check_instance(Fout, Field)
+        check_instance(vorticity, Field)
+        check_instance(Fext, ExternalForce)
         check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
+        
+        # Pregenerate parameters so that we can directly store them in self.
+        default_pbasename  = 'curl_{}'.format(Fext.name)
+        default_ppbasename = '{}_{}'.format(nabla, Fext.pretty_name)
+        pbasename  = first_not_None(pbasename,  default_pbasename)
+        ppbasename = first_not_None(ppbasename, default_ppbasename)
+        parameters = MinMaxFieldStatisticsBase.build_parameters(field=field, 
+                components=components, all_quiet=all_quiet,
+                Fmin=Fmin, Fmax=Fmax, Finf=Finf, 
+                pbasename=pbasename, ppbasename=ppbasename)
+
+        (Fmin, Fmax, Finf) = tuple(parameters[k] for k in ('Fmin', 'Fmax', 'Finf'))
 
         super(Curl, self).__init__(Fin=Fin, Fout=Fout, variables=variables,
                 candidate_input_tensors=(Fin,),
@@ -80,35 +190,6 @@ class SymbolicSpectralExternalForce(SpectralCurl):
                 implementation=implementation, 
                 base_kwds=base_kwds, **kwds)
 
-
-class SpectralCurl(Curl):
-    @classmethod
-    def implementations(cls):
-        from hysop.backend.host.python.operator.curl import PythonSpectralCurl
-        from hysop.backend.device.opencl.operator.curl import OpenClSpectralCurl
-        __implementations = {
-                Implementation.PYTHON: PythonSpectralCurl,
-                Implementation.OPENCL: OpenClSpectralCurl,
-        }
-        return __implementations
-    
-    @classmethod
-    def default_implementation(cls):
-        return Implementation.PYTHON
-
-
-class FiniteDifferencesCurl(Curl):
-    @classmethod
-    def implementations(cls):
-        raise NotImplementedError
-        from hysop.backend.host.python.operator.curl   import PythonFiniteDifferencesCurl
-        from hysop.backend.device.opencl.operator.curl import OpenClFiniteDifferencesCurl
-        __implementations = {
-                Implementation.PYTHON: PythonFiniteDifferencesCurl,
-                Implementation.OPENCL: OpenClFiniteDifferencesCurl,
-        }
-        return __implementations
-    
-    @classmethod
-    def default_implementation(cls):
-        return Implementation.PYTHON
+        self.Fmin = Fmin
+        self.Fmax = Fmax
+        self.Finf = Finf
diff --git a/hysop/operators.py b/hysop/operators.py
index 3fb77a513..ba66000fe 100644
--- a/hysop/operators.py
+++ b/hysop/operators.py
@@ -32,7 +32,7 @@ from hysop.operator.min_max import MinMaxFieldStatistics,
 
 from hysop.operator.gradient import Gradient, MinMaxGradientStatistics
 from hysop.operator.curl     import Curl, SpectralCurl
-from hysop.operator.external_force import SymbolicSpectralExternalForce
+from hysop.operator.external_force import SymbolicExternalForce, SpectralExternalForce
 
 from hysop.numerics.splitting.strang import StrangSplitting
 from hysop.operator.directional.symbolic_dir   import DirectionalSymbolic
diff --git a/hysop/tools/interface.py b/hysop/tools/interface.py
new file mode 100644
index 000000000..be755838b
--- /dev/null
+++ b/hysop/tools/interface.py
@@ -0,0 +1,176 @@
+
+from abc import ABCMeta, abstractmethod
+from hysop.tools.types import check_instance, first_not_None, to_tuple
+from hysop.tools.numpywrappers import npw
+
+
+class SymbolContainerI(object):
+    __metaclass__ = ABCMeta
+    
+    def _get_symbol(self):
+        """
+        Return a Symbol that can be used to compute symbolic expressions
+        referring to this continuous field.
+        """
+        assert hasattr(self, '_symbol'), 'Symbol has not been defined.'
+        return self._symbol
+    
+    symbol = property(_get_symbol)
+    s = property(_get_symbol)
+
+
+class NamedObjectI(object):
+    __metaclass__ = ABCMeta
+    
+    def __new__(cls, name, pretty_name=None, latex_name=None, var_name=None, **kwds):
+        """
+        Create an abstract named object that contains a symbolic value.
+        name : string
+            A name for the field.
+        pretty_name: string or unicode, optional.
+            A pretty name used for display whenever possible (unicode supported).
+            Defaults to name.
+        kwds: dict
+            Keywords arguments for base class.
+        """
+        
+        obj = super(NamedObjectI, cls).__new__(cls, **kwds)
+        obj.rename(name=name, pretty_name=pretty_name, 
+                    latex_name=latex_name, var_name=var_name)
+        return obj
+    
+    def rename(self, name, pretty_name=None, latex_name=None):
+        """Change the names of this object."""
+        check_instance(name, str)
+        check_instance(pretty_name, (str,unicode), allow_none=True)
+        check_instance(latex_name, str, allow_none=True)
+        if isinstance(pretty_name, unicode):
+            pretty_name = pretty_name.encode('utf-8')
+        check_instance(pretty_name, str)
+        
+        pretty_name = first_not_None(pretty_name, name)
+        latex_name  = first_not_None(latex_name, name)
+
+        self._name   = name
+        self._pretty_name = pretty_name
+        self._latex_name = latex_name
+        
+    def _get_name(self):
+        """Return the name of this field."""
+        return self._name
+    def _get_pretty_name(self):
+        """Return the pretty name of this field."""
+        return self._pretty_name
+    def _get_latex_name(self):
+        """Return the latex name of this field."""
+        return self._latex_name
+    
+    def __str__(self):
+        return self.long_description()
+    
+    @abstractmethod
+    def short_description(self):
+        """Short description of this field as a string."""
+        pass
+
+    @abstractmethod
+    def long_description(self):
+        """Long description of this field as a string."""
+        pass
+
+    name = property(_get_name)
+    pretty_name = property(_get_pretty_name)
+    latex_name = property(_get_latex_name)
+
+
+class NamedScalarContainerI(NamedObjectI, SymbolContainerI):
+    @property
+    def ndim(self):
+        """Number of dimensions of this this tensor."""
+        return 0
+    
+    def _get_var_name(self):
+        """Return the variable name of this field."""
+        return self._var_name
+    
+    def rename(self, name, pretty_name=None, 
+                latex_name=None, var_name=None):
+        """Change the names of this object."""
+        super(NamedScalarContainerI, self).rename(name=name, 
+                pretty_name=pretty_name, latex_name=latex_name)
+        self.check_and_set_varname(first_not_None(var_name, self._name))
+    
+    def check_and_set_varname(self, var_name):
+        check_instance(var_name, str, allow_none=True)
+
+        msg='Invalid variable name {}.'.format(var_name)
+        if var_name[0] in tuple(str(x) for x in range(10)):
+            raise RuntimeError(msg)
+        for c in '/*+-=|&()[]{}-!?:;,\'"#$^%<>@':
+            if c in var_name:
+                raise RuntimeError(msg)
+        self._var_name = var_name
+    
+    var_name = property(_get_var_name)
+
+
+class NamedTensorContainerI(NamedObjectI, SymbolContainerI):
+    def __new__(cls, contained_objects, **kwds):
+        check_instance(contained_objects, npw.ndarray)
+        obj = super(NamedTensorContainerI, cls).__new__(cls, **kwds)
+        obj._contained_objects = contained_objects
+        return obj
+    
+    def rename(self, name, pretty_name=None, 
+                latex_name=None, var_name=None):
+        """Change the names of this object."""
+        assert (var_name is None), 'Tensor do not have variable names.'
+        super(NamedTensorContainerI, self).rename(name=name, 
+                pretty_name=pretty_name, latex_name=latex_name)
+    
+    @property
+    def size(self):
+        """Full size of this container as if it was a 1D tensor."""
+        return self._contained_objects.size
+
+    @property
+    def shape(self):
+        """Shape of this tensor."""
+        return self._contained_objects.shape
+    
+    @property
+    def ndim(self):
+        """Number of dimensions of this this tensor."""
+        return self._contained_objects.ndim
+
+    def new_empty_array(self, dtype=object):
+        """Return a new empty array of the same shape as self."""
+        if (dtype is object):
+            array = npw.empty(shape=self.shape, dtype=dtype)
+            array[...] = None
+        else:
+            array = npw.zeros(shape=self.shape, dtype=dtype)
+        return array
+
+    def iter_fields(self):
+        """Return an iterator on unique scalar object along with 1d index."""
+        for (i,obj) in enumerate(self._contained_objects.ravel()):
+            yield (i,obj)
+
+    def nd_iter(self):
+        """Return an nd-indexed iterator of contained objects."""
+        for idx in npw.ndindex(*self._contained_objects.shape):
+            yield (idx, self._contained_objects[idx])
+    
+    def __iter__(self):
+        """Return an iterator on unique scalar objects."""
+        return self._contained_objects.ravel().__iter__()
+
+    def __contains__(self, obj):
+        """Check if a scalar object is contained in self."""
+        return obj in self._contained_objects
+    
+    @abstractmethod
+    def __getitem__(self, slc):
+        pass
+
-- 
GitLab