polynomial

hysop/numerics/interpolation/interpolation.py

@@ -2,8 +2,17 @@
 from hysop.tools.enum import EnumFactory
 
 Interpolation = EnumFactory.create('Interpolation',
-        ['LINEAR', 'CUBIC', 'CHEBYSHEV',
-         'L4_4', 'M4', 'Mp4'])
+        ['LINEAR',       # requires 0 ghosts
+         'CUBIC_FDC2',   # requires 1 ghosts (estimate derivatives with 2nd order centered fd)
+         'CUBIC_FDC4',   # requires 2 ghosts (estimate derivatives with 4th order centered fd)
+         'QUINTIC_FDC2', # requires 2 ghosts
+         'QUINTIC_FDC4'  # requires 3 ghosts
+         'SEPTIC_FDC2',  # requires 3 ghosts  
+         'SEPTIC_FDC4'   # requires 4 ghosts
+         'CUBIC',        # derivatives order is specified with SpaceDiscretization instead
+         'QUINTIC',      # derivatives order is specified with SpaceDiscretization instead
+         'SEPTIC',       # derivatives order is specified with SpaceDiscretization instead
+         'L4_4', 'M4', 'Mp4']) # scales interpolators
 
 class MultiScaleInterpolation(object):
     pass

hysop/numerics/interpolation/polynomial.py

+import itertools as it
+import numpy as np
+import scipy as sp
+import sympy as sm
+#import sympy.abc, sympy.polys, sympy.solvers, sympy.functions
+
+from hysop.numerics.stencil.stencil_generator import StencilGenerator, MPQ
+from hysop.tools.sympy_utils import tensor_xreplace
+
+class PolynomialInterpolator(object):
+
+    def __init__(self, dim, deg, fd_order, verbose=False):
+        """
+        Create a PolynomialInterpolator.
+
+        Parameters
+        ----------
+        dim: int
+            Number of dimensions to interpolate.
+        deg: int
+            Polynomial degree (1=linear, 3=cubic, 5=quintic, 7=septic, ...)
+            Degree should be odd.
+        fd_order: int
+            Order of centered finite differences used to compute derivatives.
+            Will affect the number of ghosts of the method.
+            Should be even because this interpolator only use centered dinite differences.
+        verbose: bool
+            Enable or disabled verbosity, default to False.
+
+        Attributes
+        ----------
+        dim: int
+            Number of dimensions to interpolate.
+        deg: int
+            Polynomial degree (1=linear, 3=cubic, 5=quintic, 7=septic, ...)
+            Degree should be odd: deg=2k+1 
+        fd: int
+            Order of centered finite differences used to compute derivatives.
+        p: int
+            Corresponds to deg+1.
+            The number of polynomial coefficients corresponds to p**dim.
+        k: int
+            Max derivative order required to compute the polynomial interpolator coefficients.
+            Corresponds to (deg-1)/2.
+                linear:  0
+                cubic:   1
+                quintic: 2
+                septic:  3
+                nonic:   4
+        ghosts: int
+            Return the number of ghosts required by the interpolator on each axis.
+            Corresponds to (k>0)*[fd_order//2 - 1 + (k+1)//2]
+                      deg    k   (k+1)/2  | FDC2  FDC4  FDC6
+            linear:    1     0      0     |  0     0     0
+            cubic:     3     1      1     |  1     2     3
+            quintic:   5     2      1     |  1     2     3
+            septic:    7     3      2     |  2     3     4
+            nonic:     9     4      2     |  2     3     4
+        n: int
+            Corresponds to 2*ghosts+1, the number of required nodes to generate the
+            polynomial coefficients (in each direction).
+
+        See Also
+        --------
+        :class:`PolynomialPatchInterpolator`: Precompute weights for fixed subgrid interpolation.
+        """
+        
+        assert dim>0, 'dim<=0'
+        assert deg%2==1, 'deg % 2 != 1'
+        assert fd_order%2==0, 'fd_order % 2 != 0'
+
+        dim = dim
+        k = (deg-1)/2
+
+        if (k>0):
+            required_ghosts = (k+1)/2
+            ghosts = (fd_order/2) - 1 + required_ghosts
+        else:
+            ghosts = 0
+
+        p = deg+1
+        n = 2*ghosts+1
+
+        self.dim             = dim
+        self.deg             = deg
+        self.p               = p
+        self.k               = k
+        self.n               = n
+        self.fd              = fd_order
+        self.ghosts          = ghosts
+        self.verbose         = verbose
+        
+        self._collect_stencils()
+        self._build_system()
+
+    def _collect_stencils(self):
+        dim     = self.dim
+        ghosts  = self.ghosts
+        verbose = self.verbose
+        order   = self.fd
+
+        k = self.k
+        n = self.n
+        ghosts = self.ghosts
+        
+        if verbose:
+            print '\nCollecting 1D stencils:'
+        
+        S = []
+        SG = StencilGenerator()
+        for i in xrange(k+1):
+            if (k==0):
+                Si = (0,)*ghosts + (1,) + (0,)*ghosts
+            else:
+                Si = SG.generate_exact_stencil(dx=1, dim=1, dtype=MPQ,
+                        derivative=i, origin=ghosts, order=order, shape=n)
+            print Si
+            S.append(Si.stencil)
+            if verbose:
+                print '  {}-th derivative: {}'.format(i,Si)
+        self.S = S
+
+    def _build_stencil(self, dvec):
+        dvec = np.asarray(dvec)    
+        assert dvec.size == self.dim, 'dvec.size != dim'
+        assert (dvec>=0).all(), 'dvec < 0    => {}'.format(dvec)
+        assert (dvec<=self.k).all(), 'dvec > dmax => {}'.format(dvec)
+
+        Sd = self.S[dvec[0]].copy()
+        for di in dvec[1:]:
+            Sdi = self.S[di]
+            Sd = np.tensordot(Sd, Sdi, axes=0)
+        for i in xrange(0,self.dim-1):
+            Sd = np.flip(Sd, axis=i)
+        return Sd
+
+    def _build_system(self):
+        dim    = self.dim
+        deg      = self.deg
+        maxd   = self.k
+        ghosts = self.ghosts
+        verbose = self.verbose
+
+        forigin = (ghosts,)*dim
+        fshape  = (2*(ghosts+1),)*dim
+        pshape  = (deg+1,)*dim
+        
+        xvals,  xvars  = tensor_symbol('x',shape=(dim,))
+        fvals,  fvars  = tensor_symbol('F',fshape,forigin)
+        pvals,  pvars  = tensor_symbol('C',pshape)
+
+        self.xvals,  self.xvars  = xvals, xvars
+        self.fvals,  self.fvars  = fvals, fvars
+        self.pvals,  self.pvars  = pvals, pvars
+
+        P0 = 0
+        for idx in it.product(range(0,deg+1), repeat=dim):
+            val = pvals[idx]
+            for i in xrange(dim):
+                val*= pow(xvals[i],idx[i])
+            P0 += val
+        self.P0 = P0
+        
+        if verbose:
+            print '\nGenerating variables:'
+            print '  *space vars: '
+            print xvals
+            print '  *grid values:'
+            print fvals
+            print '  *polynomial coefficients:'
+            print pvals
+            print '  *polynomial patch:'
+            print P0
+            print '\nBuilding system...'
+        
+        eqs = []
+        for dvec in it.product( range(0,maxd+1), repeat=dim ):
+            if verbose:
+                print '  => derivative {}'.format(dvec)
+            
+            dP0 = P0
+            for i,deg in enumerate(dvec):
+                dP0 = sm.diff(dP0, xvals[i], deg)
+            
+            stencil = self._build_stencil(dvec)
+
+            for idx in it.product( range(ghosts,ghosts+2), repeat=dim):
+                if verbose:
+                    print '    -> point {}'.format(idx)
+                
+                pos = np.asarray(idx)-ghosts
+                pos = dict(zip(xvals, pos))
+                eq = dP0.xreplace(pos)
+
+                for offset in it.product( range(-ghosts, ghosts+1), repeat=dim):
+                    fidx = tuple(np.asarray(idx)+np.asarray(offset))
+                    sidx = tuple(np.asarray(offset)+ghosts)
+                    eq -= fvals[fidx]*stencil[sidx]
+                
+                eqs.append(eq)
+                if verbose:
+                    print '        {}'.format(eq)
+
+        if verbose:
+            print '\nSolving system...'
+
+        sol = sm.solve(eqs, pvars, simplify=False)
+        for pvar in pvars:
+            psol = sol[pvar]
+            if verbose:
+                print '{} = {}'.format(pvar,psol)
+        self.sol = sol
+        
+        if verbose:
+            print '\nBuilding matrix...'
+        M = np.empty(shape=(pvals.size,fvals.size), dtype=object)
+        for i,pvar in enumerate(pvars):
+            psol = sol[pvar]
+            for j,fvar in enumerate(fvars):
+                M[i,j] = psol.coeff(fvar)
+        self.M = M
+        self.Mflt = M.astype(float)
+        if verbose:
+            original = np.get_printoptions()
+            np.set_printoptions(precision=2, suppress=True, linewidth=1000, threshold=10000)
+            np.set_printoptions(**original)
+            #print self.Mflt
+
+    def interpolate(self, fvals):
+        fvals = np.asarray(fvals)
+        assert fvals.shape == self.fvals.shape
+        #fvars = dict(zip(self.fvars,fvals.ravel()))
+        #cvars = self.sol.copy()
+        #for k,v in cvars.iteritems():
+            #cvars[k] = v.xreplace(fvars)
+        pvals = self.Mflt.dot(fvals.ravel())
+        P0 = self.P0.xreplace(dict(zip(self.pvars, pvals)))
+        return sm.utilities.lambdify(self.xvars[::-1], P0)
+
+
+
+if __name__ == '__main__':
+    P1 = PolynomialInterpolator(dim=2, deg=1, fd_order=2, verbose=True)