Fix remeshing + Horner algorithm for sequential code.

parmepy/new_ParMePy/Utils/M6Prime.py

@@ -32,35 +32,36 @@ class M6Prime(RemeshingMethod):
         ind[...] = (np.round(np.abs(a0))) * np.sign(a0)
         ind[..., dir] = np.floor(a0[..., dir])
         x0 = (ppos[..., dir] - ind[..., dir] * self.grid.elementSize[dir] - self.grid.min[dir]) / self.grid.elementSize[dir]
-        am2 = (-(x0) * (5. * (x0 + 2.) - 8.) * (x0 - 1.) ** 3) / 24.
-        am = x0 * (x0 - 1.) * (25. * (x0 + 1.) ** 3 - 114. * (x0 + 1.) ** 2 + 153. * (x0 + 1.) - 48.) / 24.
-        a0 = -(x0 - 1.) * (25. * x0 ** 4 - 38. * x0 ** 3 - 3. * x0 ** 2 + 12. * x0 + 12.) / 12.
-        ap = x0 * (25. * (1. - x0) ** 4 - 38. * (1. - x0) ** 3 - 3. * (1. - x0) ** 2 + 12. * (1. - x0) + 12.) / 12.
-        ap2 = (1. - x0) * (-x0) * (25. * (2. - x0) ** 3 - 114. * (2. - x0) ** 2 + 153. * (2. - x0) - 48.) / 24.
-        ap3 = -(1. - x0) * ((5. * (3. - x0) - 8.) * (-x0) ** 3) / 24.
+        am2 = x0 * (x0 * (x0 * (x0 * (-5. * x0 + 13.) - 9.) - 1.) + 2.) / 24.
+        am = x0 * (x0 * (x0 * (x0 * (25. * x0 - 64.) + 39.) + 16.) - 16.) / 24.
+        a0 = x0 * x0 * (x0 * (x0 * (-50. * x0 + 126.) - 70.) - 30.) / 24. + 1.
+        ap = x0 * (x0 * (x0 * (x0 * (50. * x0 - 124.) + 66.) + 16.) + 16.) / 24.
+        ap2 = x0 * (x0 * (x0 * (x0 * (-25. * x0 + 61.) - 33.) - 1.) - 2.) / 24.
+        ap3 = x0 * x0 * x0 * (x0 * (5. * x0 - 12.) + 7.) / 24.
 
         ind[..., dir] = (ind[..., dir] - 2 + self.grid.elementNumber[dir]) % self.grid.elementNumber[dir]
-        for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), am2.reshape((-1))):
+
+        for i_index, ps, w in zip(ind.reshape((-1,self.grid.dimension)), pscal.reshape((-1)), am2.reshape((-1))):
             self.gvalues.values[tuple(i_index)] += ps * w
 
         ind[..., dir] = (ind[..., dir] + 1) % self.grid.elementNumber[dir]
-        for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), am.reshape((-1))):
+        for i_index, ps, w in zip(ind.reshape((-1,self.grid.dimension)), pscal.reshape((-1)), am.reshape((-1))):
             self.gvalues.values[tuple(i_index)] += ps * w
 
         ind[..., dir] = (ind[..., dir] + 1) % self.grid.elementNumber[dir]
-        for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), a0.reshape((-1))):
+        for i_index, ps, w in zip(ind.reshape((-1,self.grid.dimension)), pscal.reshape((-1)), a0.reshape((-1))):
             self.gvalues.values[tuple(i_index)] += ps * w
 
         ind[..., dir] = (ind[..., dir] + 1) % self.grid.elementNumber[dir]
-        for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), ap.reshape((-1))):
+        for i_index, ps, w in zip(ind.reshape((-1,self.grid.dimension)), pscal.reshape((-1)), ap.reshape((-1))):
             self.gvalues.values[tuple(i_index)] += ps * w
 
         ind[..., dir] = (ind[..., dir] + 1) % self.grid.elementNumber[dir]
-        for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), ap2.reshape((-1))):
+        for i_index, ps, w in zip(ind.reshape((-1,self.grid.dimension)), pscal.reshape((-1)), ap2.reshape((-1))):
             self.gvalues.values[tuple(i_index)] += ps * w
 
         ind[..., dir] = (ind[..., dir] + 1) % self.grid.elementNumber[dir]
-        for i_index, ps, w in zip(ind.reshape((-1,2)), pscal.reshape((-1)), ap3.reshape((-1))):
+        for i_index, ps, w in zip(ind.reshape((-1,self.grid.dimension)), pscal.reshape((-1)), ap3.reshape((-1))):
             self.gvalues.values[tuple(i_index)] += ps * w
 
     def __str__(self):