Add M4 and M8 remeshing formula

hysop/numerics/remesh/remesh.py

 from hysop.constants import __VERBOSE__, __DEBUG__
-from hysop.tools.enum  import EnumFactory
+from hysop.tools.enum import EnumFactory
 from hysop.tools.types import check_instance
 from hysop.numerics.remesh.kernel_generator import Kernel, SymmetricKernelGenerator
 
-Remesh = EnumFactory.create('Remesh',
-    ['L1_0', 'L2_1','L2_2','L4_2','L4_4','L6_4','L6_6','L8_4',        # lambda remesh kernels
-     'L2_1s','L2_2s','L4_2s','L4_4s','L6_4s','L6_6s','L8_4s', # splitted lambda remesh kernels
-     'Mp4', 'Mp6', 'Mp8', # Mprimes kernels: Mp4 = M'4 = L2_1 and Mp6 = M'6 = L4_2
+Remesh = EnumFactory.create(
+    'Remesh',
+    ['L1_0', 'L2_1', 'L2_2', 'L4_2', 'L4_4', 'L6_4', 'L6_6', 'L8_4',        # lambda remesh kernels
+     'L2_1s', 'L2_2s', 'L4_2s', 'L4_4s', 'L6_4s', 'L6_6s', 'L8_4s',  # splitted lambda remesh kernels
+     'Mp4', 'Mp6', 'Mp8',  # Mprimes kernels: Mp4 = M'4 = L2_1 and Mp6 = M'6 = L4_2
+     'M4', 'M8',  # M kernels
      'O2', 'O4',          # Corrected kernels, allow a large CFL number
      'L2',                # Corrected and limited lambda 2
      ])
 
+
 class RemeshKernelGenerator(SymmetricKernelGenerator):
-    def configure(self,n):
-        return super(RemeshKernelGenerator,self).configure(n=n, H=None)
+    def configure(self, n):
+        return super(RemeshKernelGenerator, self).configure(n=n, H=None)
+
 
 class RemeshKernel(Kernel):
 
     def __init__(self, moments, regularity,
-            verbose = __DEBUG__,
-            split_polys=False,
-            override_cache=False):
+                 verbose=__DEBUG__,
+                 split_polys=False,
+                 override_cache=False):
 
         generator = RemeshKernelGenerator(verbose=verbose)
         generator.configure(n=moments)
 
         kargs = generator.solve(r=regularity, override_cache=override_cache,
-                        split_polys=split_polys, no_wrap=True)
+                                split_polys=split_polys, no_wrap=True)
 
-        super(RemeshKernel,self).__init__(**kargs)
+        super(RemeshKernel, self).__init__(**kargs)
 
     @staticmethod
     def from_enum(remesh):
         check_instance(remesh, Remesh)
         remesh = str(remesh)
-        assert remesh[0]=='L' and (remesh!='L2'), \
-                'Only lambda remesh kernels are supported.'
-        remesh=remesh[1:]
-        if remesh[-1] == 's':
-            remesh = remesh[:-1]
-            split_polys=True
+        assert remesh[0] == 'L' and (remesh != 'L2') or (remesh in ('M4', 'M8')), \
+            'Only lambda remesh kernels are supported.'
+        if remesh in ('M4', 'M8'):
+            # given M4 or M8 kernels
+            from hysop.deps import sm
+            x = sm.abc.x
+            if remesh == 'M4':
+                M4 = (sm.Poly((1/sm.Rational(6))*((2-x)**3-4*(1-x)**3), x),
+                      sm.Poly((1/sm.Rational(6))*((2-x)**3), x))
+                return Kernel(n=2, r=4, deg=3, Ms=2, Mh=None, H=None, remesh=True,
+                              P=(M4[1].subs(x, -x), M4[0].subs(x, -x), M4[0], M4[1]))
         else:
-            split_polys=False
-        remesh = [int(x) for x in remesh.split('_')]
-        assert len(remesh) == 2
-        assert remesh[0] >= 1
-        assert remesh[1] >= 0
-        return RemeshKernel(remesh[0], remesh[1], split_polys=split_polys)
+            remesh = remesh[1:]
+            if remesh[-1] == 's':
+                remesh = remesh[:-1]
+                split_polys = True
+            else:
+                split_polys = False
+            remesh = [int(x) for x in remesh.split('_')]
+            assert len(remesh) == 2
+            assert remesh[0] >= 1
+            assert remesh[1] >= 0
+            return RemeshKernel(remesh[0], remesh[1], split_polys=split_polys)
 
     def __str__(self):
         return 'RemeshKernel(n={}, r={}, split={})'.format(self.n, self.r, self.poly_splitted)
 
 
-if __name__=='__main__':
+if __name__ == '__main__':
     import numpy as np
     from matplotlib import pyplot as plt
 
-    for i in xrange(1,5):
+    for i in xrange(1, 5):
         p = 2*i
         kernels = []
-        for r in [1,2,4,8]:
+        for r in [1, 2, 4, 8]:
             try:
-                kernel = RemeshKernel(p,r)
+                kernel = RemeshKernel(p, r)
                 kernels.append(kernel)
             except RuntimeError:
-                print 'Solver failed for p={} and r={}.'.format(p,r)
+                print 'Solver failed for p={} and r={}.'.format(p, r)
 
-        if len(kernels)==0:
+        if len(kernels) == 0:
             continue
         k0 = kernels[0]
 
         fig = plt.figure()
         plt.xlabel(r'$x$')
-        plt.ylabel(r'$\Lambda_{'+'{},{}'.format(p,'r')+'}$')
-        X = np.linspace(-k0.Ms-1,+k0.Ms+1,1000)
-        s = plt.subplot(1,1,1)
-        for i,k in enumerate(kernels):
-            s.plot(X,k(X),label=r'$\Lambda_{'+'{},{}'.format(p,k.r)+'}$')
-        s.plot(k0.I,k0.H,'or')
+        plt.ylabel(r'$\Lambda_{'+'{},{}'.format(p, 'r')+'}$')
+        X = np.linspace(-k0.Ms-1, +k0.Ms+1, 1000)
+        s = plt.subplot(1, 1, 1)
+        for i, k in enumerate(kernels):
+            s.plot(X, k(X), label=r'$\Lambda_{'+'{},{}'.format(p, k.r)+'}$')
+        s.plot(k0.I, k0.H, 'or')
         axe_scaling = 0.10
         ylim = s.get_ylim()
         Ly = ylim[1] - ylim[0]
-        s.set_ylim(ylim[0]-axe_scaling*Ly,ylim[1]+axe_scaling*Ly)
+        s.set_ylim(ylim[0]-axe_scaling*Ly, ylim[1]+axe_scaling*Ly)
         s.legend()
 
         plt.show(block=True)