diff --git a/pystencils/fd/derivation.py b/pystencils/fd/derivation.py
index bfb4f7393a20f1423ea0c3da829a09b8e89873a8..5a9d79e8fe4038a090c8da84507747ed97633b86 100644
--- a/pystencils/fd/derivation.py
+++ b/pystencils/fd/derivation.py
@@ -1,10 +1,12 @@
import warnings
from collections import defaultdict
+import itertools
import numpy as np
import sympy as sp
from pystencils.field import Field
+from pystencils.stencil import direction_string_to_offset
from pystencils.sympyextensions import multidimensional_sum, prod
from pystencils.utils import LinearEquationSystem, fully_contains
@@ -228,3 +230,120 @@ class FiniteDifferenceStencilDerivation:
def __repr__(self):
return "Finite difference stencil of accuracy {}, isotropic error: {}".format(self.accuracy,
self.is_isotropic)
+
+
+class FiniteDifferenceStaggeredStencilDerivation:
+ """Derives a finite difference stencil for application at a staggered position
+
+ Args:
+ neighbor: the neighbor direction string or vector at whose staggered position to calculate the derivative
+ dim: how many dimensions (2 or 3)
+ derivative: a tuple of directions over which to perform derivatives
+ """
+
+ def __init__(self, neighbor, dim, derivative=tuple()):
+ if type(neighbor) is str:
+ neighbor = direction_string_to_offset(neighbor)
+ if dim == 2:
+ assert neighbor[dim:] == 0
+ assert derivative is tuple() or max(derivative) < dim
+ neighbor = sp.Matrix(neighbor[:dim])
+ pos = neighbor / 2
+
+ def unitvec(i):
+ """return the `i`-th unit vector in three dimensions"""
+ a = np.zeros(dim, dtype=int)
+ a[i] = 1
+ return a
+
+ def flipped(a, i):
+ """return `a` with its `i`-th element's sign flipped"""
+ a = a.copy()
+ a[i] *= -1
+ return a
+
+ # determine the points to use, coordinates are relative to position
+ points = []
+ if np.linalg.norm(neighbor, 1) == 1:
+ main_points = [neighbor / 2, neighbor / -2]
+ elif np.linalg.norm(neighbor, 1) == 2:
+ nonzero_indices = [i for i, v in enumerate(neighbor) if v != 0 and i < dim]
+ main_points = [neighbor / 2, neighbor / -2, flipped(neighbor / 2, nonzero_indices[0]),
+ flipped(neighbor / -2, nonzero_indices[0])]
+ else:
+ main_points = [neighbor.multiply_elementwise(sp.Matrix(c) / 2)
+ for c in itertools.product([-1, 1], repeat=3)]
+ points += main_points
+ zero_indices = [i for i, v in enumerate(neighbor) if v == 0 and i < dim]
+ for i in zero_indices:
+ points += [point + sp.Matrix(unitvec(i)) for point in main_points]
+ points += [point - sp.Matrix(unitvec(i)) for point in main_points]
+ points_tuple = tuple([tuple(p) for p in points])
+ self._stencil = points_tuple
+
+ # determine the stencil weights
+ if len(derivative) == 0:
+ weights = None
+ else:
+ derivation = FiniteDifferenceStencilDerivation(derivative, points_tuple).get_stencil()
+ if not derivation.accuracy:
+ raise Exception('the requested derivative cannot be performed with the available neighbors')
+ weights = derivation.weights
+
+ # if the weights are underdefined, we can choose the free symbols to find the sparsest stencil
+ free_weights = set(itertools.chain(*[w.free_symbols for w in weights]))
+ if len(free_weights) > 0:
+ zero_counts = defaultdict(list)
+ for values in itertools.product([-1, -sp.Rational(1, 2), 0, 1, sp.Rational(1, 2)],
+ repeat=len(free_weights)):
+ subs = {free_weight: value for free_weight, value in zip(free_weights, values)}
+ weights = [w.subs(subs) for w in derivation.weights]
+ if not all(a == 0 for a in weights):
+ zero_count = sum([1 for w in weights if w == 0])
+ zero_counts[zero_count].append(weights)
+ best = zero_counts[max(zero_counts.keys())]
+ if len(best) > 1: # if there are multiple, pick the one that contains a nonzero center weight
+ center = [tuple(p + pos) for p in points].index((0, 0, 0))
+ best = [b for b in best if b[center] != 0]
+ if len(best) > 1:
+ raise NotImplementedError("more than one suitable set of weights found, don't know how to proceed")
+ weights = best[0]
+ assert weights
+
+ points_tuple = tuple([tuple(p + pos) for p in points])
+ self._points = points_tuple
+ self._weights = weights
+
+ @property
+ def points(self):
+ """return the points of the stencil"""
+ return self._points
+
+ @property
+ def stencil(self):
+ """return the points of the stencil relative to the staggered position specified by neighbor"""
+ return self._stencil
+
+ @property
+ def weights(self):
+ """return the weights of the stencil"""
+ assert self._weights is not None
+ return self._weights
+
+ def visualize(self):
+ if self._weights is None:
+ ws = None
+ else:
+ ws = np.array([w for w in self.weights if w != 0], dtype=float)
+ pts = np.array([p for i, p in enumerate(self.points) if self.weights[i] != 0], dtype=int)
+ from pystencils.stencil import plot
+ plot(pts, data=ws)
+
+ def apply(self, field):
+ if field.index_dimensions == 0:
+ return sum([field.__getitem__(point) * weight for point, weight in zip(self.points, self.weights)])
+ else:
+ total = field.neighbor_vector(self.points[0]) * self.weights[0]
+ for point, weight in zip(self.points[1:], self.weights[1:]):
+ total += field.neighbor_vector(point) * weight
+ return total
diff --git a/pystencils/field.py b/pystencils/field.py
index 6b3fa9f8c4c5fa21cd799fc9ff37d090224794af..dfe71c692ecf307753459fb219431c977546f3c6 100644
--- a/pystencils/field.py
+++ b/pystencils/field.py
@@ -441,6 +441,22 @@ class Field(AbstractField):
center = tuple([0] * self.spatial_dimensions)
return Field.Access(self, center)
+ def neighbor_vector(self, offset):
+ """Like neighbor, but returns the entire vector/tensor stored at offset."""
+ if self.spatial_dimensions == 2 and len(offset) == 3:
+ assert offset[2] == 0
+ offset = offset[:2]
+
+ if self.index_dimensions == 0:
+ return sp.Matrix([self.__getitem__(offset)])
+ elif self.index_dimensions == 1:
+ return sp.Matrix([self.__getitem__(offset)(i) for i in range(self.index_shape[0])])
+ elif self.index_dimensions == 2:
+ return sp.Matrix([[self.__getitem__(offset)(i, k) for k in range(self.index_shape[1])]
+ for i in range(self.index_shape[0])])
+ else:
+ raise NotImplementedError("neighbor_vector is not implemented for more than 2 index dimensions")
+
def __getitem__(self, offset):
if type(offset) is np.ndarray:
offset = tuple(offset)
diff --git a/pystencils_tests/test_fd_derivation.ipynb b/pystencils_tests/test_fd_derivation.ipynb
index d195585f091cb29622a3051c1468ca19bd5c9abe..f992c8ac5051dd83553cb7da3e6842c9be1aa8c3 100644
--- a/pystencils_tests/test_fd_derivation.ipynb
+++ b/pystencils_tests/test_fd_derivation.ipynb
@@ -60,17 +60,22 @@
"metadata": {},
"outputs": [
{
- "ename": "AttributeError",
- "evalue": "MutableDenseMatrix has no attribute tomatrix.",
- "output_type": "error",
- "traceback": [
- "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
- "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)",
- "\u001b[0;32m<ipython-input-4-ea41cd8e50a0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mstandard_2d_00\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_stencil\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
- "\u001b[0;32m~/pystencils/pystencils/pystencils/fd/derivation.py\u001b[0m in \u001b[0;36mas_matrix\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 185\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mdirection\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mweight\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mzip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstencil\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mweights\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 186\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mmax_offset\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mdirection\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmax_offset\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mdirection\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mweight\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 187\u001b[0;31m \u001b[0mresult\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtomatrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtomatrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 188\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"test\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 189\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mdim\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
- "\u001b[0;32m~/anaconda3/envs/pystencils/lib/python3.7/site-packages/sympy/matrices/matrices.py\u001b[0m in \u001b[0;36m__getattr__\u001b[0;34m(self, attr)\u001b[0m\n\u001b[1;32m 2139\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2140\u001b[0m raise AttributeError(\n\u001b[0;32m-> 2141\u001b[0;31m \"%s has no attribute %s.\" % (self.__class__.__name__, attr))\n\u001b[0m\u001b[1;32m 2142\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2143\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m__len__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
- "\u001b[0;31mAttributeError\u001b[0m: MutableDenseMatrix has no attribute tomatrix."
- ]
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGkAAABLCAMAAABEOFk1AAAAPFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAo1xBWAAAAE3RSTlMAMquZdlQiEEDpMESJZs3d77ts6LWnRAAAAAlwSFlzAAAOxAAADsQBlSsOGwAAAj9JREFUWAntmV17gyAMRlFQ1olf4///1wl0kpc+Qlba3QxvBnKas6TUaio6649evOtYg0CIzkp1HMO7RGJ24XvrTOPbJDHwgCajZtPF1csRk6MYmpajgnotZ8jkAAPTvLkkzHKZyn2BySEGpt1vv8HqgorJIQYm602TLe1CJocYNWkrXTKTVfmcmFyCUdNojVN0JROTSzA0+Zy6ILzOawy5l7gEo6Yk3UsVk0swahLhLRyYO6LIYTgw7bvLYy7vch6H4cCEH7XL6gkmhxiYxOquRlv5wsfkAEOTNsfVvSwSTA4wNF0XrH6lmWpq2KrXqocVaDsC6yG0lPvOuEz6l6XVG7bSLVi0aXcrZewUzzyMSDgw6d70xa/BGE36f2rzt6PxbBxhODAdUPkLN0baVjfuc88PJFyNafH59zbzTr3IFLJbc/V+pWnK3h2+0rRmH1qfNOmFHPc3R/rb3rhLklFi+rh9ngBZOs9lBiovolv564bPub8zDU40Zp4hSbiaXS46/2yinjOp3J5NyjiuxhjZ+89vsvQzJeEwJ9lvdpHzD1f4u4SGxrUJwqGpELpquZlqyteq16qHFfgHO4I2GTF5nDE5imH1oMmIsWHG5AADE7YqIDZMmBxiYMImIwSHCZNDDEzYZITgMGFyiFFT0iaD4HTC5BKMmpImIw0OYyaXYGhqPctSb7P1LGHTnRNoRp5nHwaA0b13dEZbz/KhXLkTWL0cWbvWTDUVDNX7q1+OtfthV6lcp6kmmaNFFATiG+4bKVe2WSLVAAAAAElFTkSuQmCC\n",
+ "text/latex": [
+ "$$\\left[\\begin{matrix}0 & 0 & 0\\\\1 & -2 & 1\\\\0 & 0 & 0\\end{matrix}\\right]$$"
+ ],
+ "text/plain": [
+ "⎡0 0 0⎤\n",
+ "⎢ ⎥\n",
+ "⎢1 -2 1⎥\n",
+ "⎢ ⎥\n",
+ "⎣0 0 0⎦"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
}
],
"source": [
@@ -88,9 +93,20 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 5,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Finite difference stencil of accuracy 2, isotropic error: True"
+ ]
+ },
+ "execution_count": 5,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"stencil = [(i, j) for i in (-1, 0, 1) for j in (-1, 0, 1)]\n",
"isotropic_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil)\n",
@@ -102,18 +118,48 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 6,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIIAAABNCAMAAABdY0yuAAAAYFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACACGXEAAAAH3RSTlMAMquZdlQiEEDpMES73YnNZu+F48P7apvf1/Fc0cls39hR/gAAAAlwSFlzAAAOxAAADsQBlSsOGwAAA8dJREFUaAXtmumW2yAMhcFrU6/pdJ0ufv+3LNjGBnQBkcn0NOfEPxIbpMsXGQMyEcI/ZO+XsK55bhwrWZYDq0nXiOcWsJLLehSHZHcLghA8N9tq2FoWQi5lpY763yPMut1i0QjN0fp6YoO6NdErnhuxqp8IKqzPKOi+BaNQ3/ZQ8tyIFUBoymEpqmjfR5U8N2AFEJD+e5Y9LkI7WYd8W4zcKGxDtvXJFz+cgi6HhTnZLV2EoPt7VhCEfTofy75sMxsex2qOu0BxD8FM56OauIoprufXTp1orn6hfR0Q9xDMpD8oBDKB2nL0vFLrjTYxnGyzpCceQuiEaBb1wT+u53oj6LQjuOIBBC1SL1mdYamrKtV9zrWCJR5BGBJRdX9rs4zmLroV9tWJYImHEcbSdk6et+ttS9y7A8EWDyLMeQRCbAjxp9IgOOIhhE6HNWvknXTrvCi44gShHmQhm08vdV33Wd3x60U2317iNwyKewh6Ov+sxvAv6zge1/NqO+XyPQodEPcQtGxXy6iS17S55LkBK4SQ1QUMAbPngP6FEFRulRpjjobPk47lBqwAQqPGzzk/t+a5ASuAoH9Y5hRlYsFz86wAgr4J21hnlFnfPDdgBRCuaqJu8qYozchzA1YAYVRRsMdwVgyUC8sNWAEEMY6lGp6zD54btUII2Y2/zeFxEe6dyny4fLRCaRIN9fol5+C5Ias/l5vHoRzAmC3sC+mcBEmOVVkwplgijhCSOQkCEKVaNZXpqYWKAwROTgIg9FgyJzsQEAcInJwEIOiidBSAOEDg5CQYob0mVztAnCKwchKIUE3J5AeJUwRWTgIRGDcCiVMEVk4SQGiWVG4LEh6AwMlJKMLaD5olnk0JAcQBQtfLvs58vaGQhkb28zU1NgFxgCD7H683rKB//vr9Wni7CjRYQBwgpHN0qryWmKQ1UL0XE6sngo7MMwrPKOxPCO4LZOPEWMe/eW7ECnRHsHESb3qr5bkBK4DAae+eNo+LcO9U5nzlcGQa5oQfcOMRXsAeFuZkF/8fbwTcOOEFgyQpxA2Ke1EIbJwQLVRAkxTPKiDuIZi1grdx4mnBS5CkELttreCJhxDcjROiBQpAkkKsdgRXPICgna2NE6IFCkCSQqzOFZMlHkGwNk6IFi1ASQqxOhEs8TBC5ks3lKSEEWzxIIKzcUK0QAFIUoiViYIjHkJwN06IFigASQqx2hFc8Q1hHTK3v1Zt07mcsndlOvV6I5UB+eLHX6ta/QenqlKPijDT+XWFIr8iWlCXYzwDouLrX6uqSvwFDr5NBGUlpz0AAAAASUVORK5CYII=\n",
+ "text/latex": [
+ "$$\\left[\\begin{matrix}\\frac{1}{12} & - \\frac{1}{6} & \\frac{1}{12}\\\\\\frac{5}{6} & - \\frac{5}{3} & \\frac{5}{6}\\\\\\frac{1}{12} & - \\frac{1}{6} & \\frac{1}{12}\\end{matrix}\\right]$$"
+ ],
+ "text/plain": [
+ "⎡1/12 -1/6 1/12⎤\n",
+ "⎢ ⎥\n",
+ "⎢5/6 -5/3 5/6 ⎥\n",
+ "⎢ ⎥\n",
+ "⎣1/12 -1/6 1/12⎦"
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"isotropic_2d_00_res.as_matrix()"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 7,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKQAAACaCAYAAAAnxeOcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAE+FJREFUeJztnXtwXNV5wH9nd7UPrWTJ8soWsh6uATPBIQ/KKzEloSngdcOAE9M0E6cUApmUuJPAJI1pnLSJA2RqbEKGkjKEAs0wDiUYOoAX8iyEDEljYGyS1ji2gm1JfuhaWkkrrR6rPf1jtbZkr/Z1X+de3d+MZ6x9nPt9+n66j71nvyOklHh4qILP7gA8PGYSsGpDIqHVARcBUUBU8FYJpIBdMh4bMSM2IxAJbQHwp0AtleVXiCxwBNgt47Gs3tjMQiS0ZcB5QE2Fb53iVH6zDtHCikO2SGg3AF8HgjqGGQe+IeOxp42JyjhEQvs08BUqL0wpeoCbZTz2jsHj6kIktBBwP3ClzqEOAjfJeKwn/4Dph2yR0JqBb6JPRoAQsFkktJj+qIxDJLR2YBPGywiwFNhswrh6+TT6ZQToBP5p5gNWnEP+mYHb8QOrDBrLKK4wefyLREKrNXkblfIhA8daJRKaP/+DFUI2GDxeo8Hj6cXo/E7HB9SbvI1KMTLnALnrCsAaIfWe4Js9nhNQLWfT4rH/Y5/hAR+f/3AHazvOZf9uveeZzsKtuevIy34hw9Esm5/s4ZKrh+0OxXLcmruOvOwXsiYITUum7A7DFtyau4687BfSw2MGnpAeSuEJ6aEUlt3LLsrGtUs5uDdMb1eQa9Yn+ejNQ3aHZBluzb3KvNQQ8tvP9JR+kUtxa+5V5uUdsj2UQo09pBu5Y3U7e1+PFHzu3Peluf8nhy2OyBHYI+Sa5hUVvX5n3z6TIjGPbS+WFq7Y78FpORtUU3uEfHLffjaubaO3K8SW5w5yznsn2PNqmIe/vphAjWThkgx3fv8INe65m1aQnX372PWzCD/ctohsVnDtLQNc+fGU3WFVhUE1teccstCtpZZlk2x5/jD3vXSYJR2TvPJsnS2xWcl4WrDjwSbu3tHNtsRhx8oIhtXUHiEL3Vpa3DZFuDY3fd3vlwgXXm9JCb1dZ5Psy00y3vNqhGA4y6Yb2vjaJ1rRev0lRlAXg2qqXtW7D9Tw5stRLr82xehwHccOdTIyVIcbvh2ZTkWR0sfIcCNSQv9xP8cOBfnWU92sXj/IY3cpNRu+JJmMn6S2iKOHOslOze3SzJqWQK2r7FTSx723tfDlB48QDEF2KsPkRIhk31kMapJowwDjo/XgrLqdJDXYhJQCgSSdqqOuIct5F6YJhuDiq0Z56oEmu0Msi6TWSGailfF0FCkFPn8G4Sv8ZbTTa1oCdfaQmUm46+az+NSXTrDs/EkAQpEx/IFJpBRksz5SySae2HKfEMLsWdrGk5kMMDkeBkBKH6nkQlZemqb7QBCZhbffCNHSMWlzlCURQnyefW+sYmy0LvfHJST1jf2IAnN2C9W0BPbtIU+/teTzw4E9YbZvXcT2rYuI35jkqk8OU7/wBIPaEqQUSCloWvJbwHnzB1ODDUh5qmoT4xGaWnxctjrF7fF2fALueOCojRGWywvAd5k5azy6YBAov6ZFsE/IQreW1tx45v3O2vphBrXFCAFSCk4cvcyK8Aynti4FUpAabEIISbShH38gw7oNSdZtSNodXgUsRspTR9ZI3RA+f+4Ev9yaFkGdQ/Zc+HySBU0akbpBWjq6ph+dEsJhl+HB8DiNzRoA4egwjbETBQ9zCiOEuAT4Df7ABC0dXQTDo9Qv7DdyG1YUNa17hPqFSZqWHCcQzHD95zZMP6qKlGMWbEP/71AnJ2UEjYs+8l8EghkWt3VTEzTivHc8/x8rCrrL0NGWr3yR3BfMQQ0pXzd5/P0yHhs0eRtFmSmjlLIZY3N+S8Zj1gkp47E/ADsMGu5pGY8dkFIeQhEpZTy2G0iYNHwG2GrS2GVRQEaAHwC9Bgw/Adw3a3sW9fbxAR+c/ldHdc2mfgW8NrP5khCig1x/GAC/lNKWxkzTnRcuBz5AqWZT27duonnpbv7ir58rMuQUcBR4ScZjfzQy1kqYQ8bccwmtAVhN9c2meoEXZTx2aNa4Tu8PqYqU5SKEkMBjUsqb7I6lGMVkNBO7z790o9Lh2y3YJSO4QEjwpDQSO2UElwgJnpRGYLeM4CIhwZNSDyrICC4TEjwpq0EVGcGFQoInZSWoJCO4VEjwpCwH1WQEFwsJnpTFUFFGcLmQ4ElZCFVlhHkgJHhSzkRlGWGeCAmelKC+jDCPhIT5LaUTZIR5JiTMTymdIiPMQyFhfknpJBlhngoJ80NKp8kI81hIcLeUTpQR5rmQ4E4pnSojeEIC7pLSyTKCJ+RJ3CCl02UET8hZOFlKN8gInpBn4EQp3SIjeEIWxElSuklG8IScEydI6TYZwROyKCpL6UYZwROyJCpK6VYZwROyLFSS0s0ygidk2aggpdtlBBs76E43aCrWpGhSxmNVrUpvFlLKQ0KITnK9hKaEELN6CYmEJoAgxZpN1YTAH/CLhBYusqmsjMcmZj7gBBlFQgsy905OAhMyHivaTMryZlMioV0L3Equa1Yp3gYelvFYsU5hlnN6gyt29q0A/gG4lFJ/5N37VxCpG2JRS6l+4v3A88C/sKb5/Sgq4/Qf1kYgDjSWePkE8Cpwt4zHCi69Z+lhRyS0VcAWypOR6ddtEQntA+ZFVTmzDt/h6BRSPg6swtgjThPwNxzc+z0UlXGazcAnKS0j5I4efw48Pn2EPAOrz4Ouo7LekEy//joTYtHFSSnf80HoOXCJKQs7jafDTE58Br9fSRlFQosA11Tx1qXARYWesFrI9irf12ZoFAYhpTzErZv/EYCeAysMlXI8Haavp4PauimeO3q+cQMbSjNQejWkwhSsqdVCVrs9tVYcm8nSswdp6cytDmGUlHkZfb4pWpcfAFRdA1GPPwVzsr/QPV0Bbr+mk7azc1eVmx7vPWMRR9UJ1GRo6ezi6MHl9BxYwdKz95W15Eeh3KMLak6T0VnorKf9QgK866I039huRBN1+6hWypm5n7lndCY66qnGB+P73ozwxavbeeirMdRuEV6cvJRQ/uE7n/v37mzheLfzZQRd9bRfyObWKR7Z1cV9Lx5mUPPzi6edvXB7JVLmc//2M30MHFvAmy9nHS+jznraL2QwLKmtyy3uveqjKbreqvaqTR3KlTIYlvj9IbQjHVx4ZZbjhwYsjNIcdNbTfiFTyVMxvPVahNblyi/RWxblSNl/LHLynPGd/0u6Ined9bT/ouaN/46wfWuMUCTL4rZJbvmGZndIhlHoQifPeDrMazvbee4RCNdOsLgt4IrcddbTfiGvuH6EK64fmfN5KaH3j21CNP+K3L33D1oXnAGcLiVAdspPX08H7//QFH95k+POGYUQFwA/5ewLfsrWRJBg6NREkFL1LIH9Qs7F5EQNqcFGRocb6PrdcnKxjkzPtlGHh3/dxJL20r/HWOshtN4OAMbTUXy+LIvbDpKZKP3ex+9pE2seUOnc+kIgSir5V/R1+/EHJqhrHKCuoaK1sQuhrpADx1qYGI8AkJ3Kn5dEgXfsCqkgzz4E1332zMfvvQ32vVn4PWe/B77ykI+jh5aXtY2f/+dvdURoHpLcESwzGSLZ1+IOIdc0r5jzuSf/cIxUsgmfT5C7ABuUUjZZFlsZiIR2G/CFM574zo9n/5z/0Hsmny0yiWln36nzzSd+/yEZj5WarmYZQoirgWcRQiBEiGBkhPrG3CcExeo5M6c5sF/InX372PWzCD/ctohsVnDtLQNc+fHUyefrGwcJRd4htzpqrV1h6mLmHZhs1k+4dpix0XrufQH+7c4x/AGJzy/Z+PARmpc64bZpF7CXWGuCls5P4A+cirlUPUtgv5DjacGOB5u4e0c3wTlOk5qWDEgpb7Q2MIM4/XZg9/4VCJ+kpbOL7NRyvvS9MO3n7uOFRxfwwqMN/O2mfrtDLoWUcj9woUhoy4B1s54sp55FsP9zyD2vRgiGs2y6oY2vfaIVrVfVmS2VU+zedKAmQ+vyLny+3OeUo8M+Ot81McdIzkFnPe0Xsv+4n2OHgnzrqW5Wrx/ksbtidodkCOVMlAjUZBjq7+buz0DiPxaz4v1jFkdpPDrrab+QdQ1ZzrswTTAEF181Svf+oN0h6aaSWTvnXzLKd3/axXW3wg/uWW7KzHMr0VlP+4VceWma7gNBZBbefiNES4ezb59VIuPEWG5+WqAmQ3PbEYJh42eeW43Oetp/UdPYnOWy1Sluj7fjE3DHA8p8vFExlc5n3LsrxL9vbsbng5pQltvvP0hmsrOi+ZSqobOe9gsJsG5DknUbknaHoYtqJte+5/IxvvPS7K+DZiYrn+SrGjrqaf8h2w0YOdO7mkm+LsJqITNVvk/d88rfvdZu+NcOTpcyqan6YbmeuAq6YLWQXVW+T8kZMUKIS9i+daMpXzvIS5lKwvqVx1XoulaAY0C1M3sKumB1kj+i8r3dJPC0CbHo4mSvnd//WqN1+cumbCRQk6Gx+Ztks6BIK8CZTPcfqqY2/wvsLvSEpRc1Mh7bLRLaLcCNwEpKNJsCfg88KuOxt6yIr1xmNX4aH2sWCa0V+DvgA+Tut899JTI8ADBOrnfPXGSBI8ALtHQ+BjzCHA2uFOAe4Di53j5LKN5sahj4JfCvczWdsrzZlNPR24VMCCGBx6SUN1X4vlkNrhST0jCUOgSojp0t8VToT2kFrkzKDFTozzgfpHRdQmaggox53C6lq5IxA5VkzONmKV2TiBmoKGMet0rpiiTMQGUZ87hRSscnYAZOkDGP26R0dPBm4CQZ87hJSscGbgZOlDGPW6R0ZNBm4GQZ87hBSscFbAZukDGP06V0VLBm4CYZ8zhZSscEagZulDGPU6V0RJBm4GYZ8zhRSuUDNIP5IGMep0mpdHBmMJ9kzOMkKZUNzAzmo4x5nCKlkkGZwXyWMY8TpFQuIDPwZDyF6lIqFYwZeDKeicpSKhOIGXgyzo2qUioRhBl4MpZGRSltD8AMPBnLRzUpXSekJ2PlqCSlq4T0ZKweVaR0jZCejPpRQUpXCOnJaBx2S+l4IT0ZjcdOKS1pNiUSWgi4DlgFVLTC/DQp4FfAs9Mt4HLjKiKjSGgR4Hpy3c+iRV/80hM3s6jlD1z0kV8WeZUEeoEXZDz2G8MCrZBiDa5EQusEPgacR/EudoXIAj3A8zIe2zVrmxYJ+QhwuQFDvSLjsVtBKRn9wOPAxWW9oXv/CiJ1QyxqKbcZ/FdlPPajauPTSyEpRUI7B9gOLNA5vAS+LOOx5/IPmL4rFgntvRgjI8AVIqG9WxUZp7mYcmWsjttMHLskcxy+16NfRsj10fzczAesODe4wNDRXv/5WtSREeDdJo+/VCQ0W1fAPV1KpDSypudMn/IA1gipf2WuZF+Mvp42xkYjvPzMN1FHRjAiPzW2UZRZUv7Pjz/GeDrEkXf+hIkxIxaWP3kOqv5VdnbKx8jQQibGImi97fj9owrJWD7pkVr6eltz/08t4MSRs8hmHbUIzSwp+3o6mcrUMDRg6N7bfiGHB3x8/sMdrO04l/27z9wTpIZy5ypS5ooXqT8shDDir9JaJsdDjI+e+oQhPVLPUH+gaO5q0oQQp9pJj4/WMZWZveJrqZoWwX4hw9Esm5/s4ZKrh894TkpIJZtOyiiEJJ1qQ5UVyCohumCIXH/xHKHICNEFk3Pmri4rAWblkhpsmPWKYjUtgf1C1gShaUnhBXjGRqJkpwIIIQnUTNDYfJTP3f0FKWW1a6PYhz8wRTA8CuSKWdc4UDR3RZFSPsEFq35OdMEAQmSRwMhg06wVx3TkpfaeRogs4doU9Qv7CUXya0k7qoCzqG/sZzwdRfiyhGtH7Q6naiLRMSJRjQWLTjA6XM/YSB1Sill7zSpRW8hwNE04mrY7DMMI1abxBzLU1icduajm6fh8krqGIeoahowaUm0h3YYQcNayapfXmxeoIeTGtUs5uDdMb1eQa9Yn+ejNhv3FKY9bc68yLzWE/PYzPXaHYBtuzb3KvOy/yvbwmIEae0g3csfqdva+Hin43LnvS3P/Tw5bHJEjsEfINc0rKnr9zr59JkViHtteLC1csd+D03I2qKb2CPnkvv1sXNtGb1eILc8d5Jz3TrDn1TAPf30xgRrJwiUZ7vz+EWqccjetSnb27WPXzyL8cNsislnBtbcMcOXHU3aHVRUG1dSec8hCt5Zalk2y5fnD3PfSYZZ0TPLKs9XMLHcW42nBjgebuHtHN9sShx0rIxhWU3v2kIVuLS1uO/Wz3y+xv4mC+ex5NUIwnGXTDW2EIlm+cN8xYq3OvBNlUE3Vq3r3gRrefDnK5dc6d29RLv3H/Rw7FORbT3Wzev0gj90VszskU6igpmoJmUr6uPe2Fr784BGCzpthVjF1DVnOuzBNMAQXXzVK9373nTRXWFN1hMxMwl03n8WnvnSCZedP2h2OJay8NE33gSAyC2+/EaKlw115V1FT+z6HPP3Wks8PB/aE2b51Edu3LiJ+Y5KrPumkeYKV09ic5bLVKW6Pt+MTcMcD5X4TUU0MqKl9Qha6tbTmRnfcx62EdRuSrNuQtDsMQzCgplYcsrOlX2LreHox/4vt6uVsdDwnf4dWCHlC8fH0YnY8GUC1I4eROaeBk5OVrRDyZWDcoLHGgWItSOzgF5g7i/0VGY+NlX6ZpfzEwLF+IeOxk78/04WU8dgQ8PfAcZ1DHQM2TI+nDDIe6wO+CGgmDL8L2GTCuHp5EngCmCj1wiJI4DXgn2c+aElvHwCR0ARwDtU3m9ov4zFrgq0CkdB85PIr3myqPLLAURmPHTNgLNMQCa0OWEZ1zaZ6p/+YZ49plZAeHuWgzgfjHh54Qnoohiekh1J4Qnooxf8DQobOZydRBWoAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ "<matplotlib.figure.Figure at 0x7f06c885ccf8>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
"source": [
"plt.figure(figsize=(2,2))\n",
"isotropic_2d_00_res.visualize()"
@@ -121,7 +167,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
@@ -131,18 +177,40 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 9,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "sympy.matrices.dense.MutableDenseMatrix"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"type(isotropic_2d_00_res.as_matrix())"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 10,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "sympy.matrices.dense.MutableDenseMatrix"
+ ]
+ },
+ "execution_count": 10,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"type(expected_result)"
]
@@ -156,9 +224,28 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 11,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHoAAABNCAMAAABe+IDDAAAAVFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACwjqETAAAAG3RSTlMAMquZdlQiEEDpMES73Ylmze/h0WzfvadcyfFyvev/AAAACXBIWXMAAA7EAAAOxAGVKw4bAAADhklEQVRoBe2a2XqlIBCE8YCajHuSWX3/9xyWg9LlQvfEzPLN8cKNpn5B0AJVKl2KNj3K7fe61XUuaE0/Ey+07tbI7F5fKXUbs2H3gB3xYvbLzUc0EnRn0cVsV8wlineB6DJrY5fyF9CNUtVsV8wlogfHu81w3TGVKaZUOfNvNhEv343uDPsq1bXoXvPJ16IHCflSdNMr1RTscp9XeCnpXMVYlmXLb2ZEHJpZpbv5xm83k++f3EKDOKC5KlfEPdBKhUerXfMqVBaO0X9Thfe9GXhldlG90Td+31JEHEs9Nqqa2Ghtr1LzvQUVB7SxL+2a36+1fZgNzIahFIgDegovbXaxbSC/1CAO6Lk0RmK2bB1N7Ec4iFN0NbvXgeQpbkb27UFxiq692RE4HkmFozhFq4AW9C7nzbjtA8QBPToqu9T+Plcz90pBHNDNi27bV3YDd2Z4mJ5Czza9sS3lZAFxQPdvuv/MRxeW9uVrGDUMllsGL3+AB3FASz29h4Qu0Tk7fmqMQRzRQk+/ouvZ9W94avjkZeUuLhkwANqFSTy9l/WlDuOfLtfkEvEdtMTTp2j3BsvmTQK2aJGnT9Fu1Ddlnm2p+AYt8/QrOtzrzCOBiCNa6OlXtMq3cEXFAS319B4djL2xQ6Dzfg3igJZ5eg9ejL3p0zvp0+gKxAFNYz/26N9B12OysO3Jfu25Uj89f7onokvfz0POLlniDklND2JAtJFvzzCXkgZ/7D7ea6GnjxfHy0ajAC309JHMywZRiBZ5+khWvKEARAHaifE9/YLmZ0vEt2iBp0/RvGxp1AYt8PQpmZeNRG3Qf7LCBZ4+LTYzWzJgoKUWevqI5mXDKIpWXVW0w8SfJ7h/BQhDgVw2EAf0t9fvP278qfX4FcANBfLZQBzQwiEumPp4Aw63ZPz8TjQ19YfImHAh2kkmpj4SDrcXoxNTf4iMCdeiM1YwQsP2UjQx9ZSzc3Qlmpr6HRg9dY4mE/U04/YITP02AM4Qcehci5+HPAeHYOoPouJpEAd0jPod2wf68RWA2c6oqc9menwFsK88Mtw7nXLbr8/E1O8HLGdBHNAwUb/kOt5JTf1xlE8BcYrGifqMlk0mpv48HMUpGifqz7VCKrvCUZyi/6+vAH6aw02jN3ZCv2T/XmNHDO4rQM76x9u2ii+/xNTuBxVj3FS6KnUv+M7Fs/4RvYr7X2KMUT8Bqn5A97SBuwoAAAAASUVORK5CYII=\n",
+ "text/latex": [
+ "$$\\left[\\begin{matrix}\\frac{1}{6} & \\frac{2}{3} & \\frac{1}{6}\\\\\\frac{2}{3} & - \\frac{10}{3} & \\frac{2}{3}\\\\\\frac{1}{6} & \\frac{2}{3} & \\frac{1}{6}\\end{matrix}\\right]$$"
+ ],
+ "text/plain": [
+ "⎡1/6 2/3 1/6⎤\n",
+ "⎢ ⎥\n",
+ "⎢2/3 -10/3 2/3⎥\n",
+ "⎢ ⎥\n",
+ "⎣1/6 2/3 1/6⎦"
+ ]
+ },
+ "execution_count": 11,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"isotropic_2d_11 = FiniteDifferenceStencilDerivation((1,1), stencil)\n",
"isotropic_2d_11_res = isotropic_2d_11.get_stencil(isotropic=True)\n",
@@ -168,13 +255,128 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"expected_result = sp.Matrix([[1, 4, 1], [4, -20, 4], [1, 4, 1]]) / 6\n",
"assert iso_laplacian == expected_result"
]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# stencils for staggered fields\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "half = sp.Rational(1, 2)\n",
+ "\n",
+ "fd_points_ex = (\n",
+ " (half, 0),\n",
+ " (-half, 0),\n",
+ " (half, 1),\n",
+ " (half, -1),\n",
+ " (-half, 1),\n",
+ " (-half, -1)\n",
+ ")\n",
+ "assert set(fd_points_ex) == set(FiniteDifferenceStaggeredStencilDerivation(\"E\", 2).stencil)\n",
+ "\n",
+ "fd_points_ey = (\n",
+ " (0, half),\n",
+ " (0, -half),\n",
+ " (-1,-half),\n",
+ " (-1, half),\n",
+ " (1, -half),\n",
+ " (1, half)\n",
+ ")\n",
+ "assert set(fd_points_ey) == set(FiniteDifferenceStaggeredStencilDerivation(\"N\",2).stencil)\n",
+ "\n",
+ "fd_points_c = (\n",
+ " (half, half),\n",
+ " (-half, -half),\n",
+ " (half, -half),\n",
+ " (-half, half)\n",
+ ")\n",
+ "assert set(fd_points_c) == set(FiniteDifferenceStaggeredStencilDerivation(\"NE\",2).stencil)\n",
+ "\n",
+ "assert len(FiniteDifferenceStaggeredStencilDerivation(\"E\",3).points) == 10\n",
+ "assert len(FiniteDifferenceStaggeredStencilDerivation(\"NE\",3).points) == 12\n",
+ "assert len(FiniteDifferenceStaggeredStencilDerivation(\"TNE\",3).points) == 8"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "c = ps.fields(\"c: [2D]\")\n",
+ "c3 = ps.fields(\"c3: [3D]\")\n",
+ "\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"E\", 2, (0,)).apply(c) == c[1, 0] - c[0, 0]\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"W\", 2, (0,)).apply(c) == c[0, 0] - c[-1, 0]\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"N\", 2, (1,)).apply(c) == c[0, 1] - c[0, 0]\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"S\", 2, (1,)).apply(c) == c[0, 0] - c[0, -1]\n",
+ "\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"E\", 3, (0,)).apply(c3) == c3[1, 0, 0] - c3[0, 0, 0]\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"W\", 3, (0,)).apply(c3) == c3[0, 0, 0] - c3[-1, 0, 0]\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"N\", 3, (1,)).apply(c3) == c3[0, 1, 0] - c3[0, 0, 0]\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"S\", 3, (1,)).apply(c3) == c3[0, 0, 0] - c3[0, -1, 0]\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"T\", 3, (2,)).apply(c3) == c3[0, 0, 1] - c3[0, 0, 0]\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"B\", 3, (2,)).apply(c3) == c3[0, 0, 0] - c3[0, 0, -1]\n",
+ "\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"NE\", 2, (0,)).apply(c) + \\\n",
+ " FiniteDifferenceStaggeredStencilDerivation(\"NE\", 2, (1,)).apply(c) == c[1, 1] - c[0, 0]\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"NE\", 3, (0,)).apply(c3) + \\\n",
+ " FiniteDifferenceStaggeredStencilDerivation(\"NE\", 3, (1,)).apply(c3) == c3[1, 1, 0] - c3[0, 0, 0]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAF3CAYAAABkPHbIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAE2dJREFUeJzt3X2MXeVh5/HfGYxtbPwCHrc2YIwDdhpapSnbJW22QF7YiCEJgsACWgkFpaVSQjc02apItbRRt6RpKrZJitIq2hKaRU2ysO1mlSqTVmmzZNWlf5SqadolNIaYGmgSxobBxuMZZubsH3f8Pi93Hs/c189HsvDcOWee548zfH2ec+65VV3XAYASA+2eAADda0W7JwCdoBoeWZHkp5JsTFIt0zDTSb6f5B/roUFLACypanhkIMkbkwxm6U4QJpPsqYcG/3nOcS1n0e+q4ZGfTfLpJBtaNOSzSe6qhwafbdF49LhqeOT1ST6bZOsyDfHNJB+qhwbHTv2G5Sz62swZSCsDkiTbk/xmC8ej992f5QtIklyd5Bdn+4aI0O9+Kq0NyFFXVMMj69owLj2mGh7ZmmRXC4Z662wvigj9bmObxh1Isr5NY9NbWnUMzzqOiNDvlusieqePTe9o1XE06zgiAgvZ862V+dUbLsqtl12a6zcvvGzw1BOr8oGrtufGi3bmA1dtz1NPrGrBLOFkjz6wMR+8entuuGBnPv4LW+bd9ku/c15uf/2lufmSy/Jbd23JxJGmwyQisJAVZ9f5ufcczN2//YMFt50YT+6788Jcc9MreWTPnrztltHcd+eFmRhvwUThBJu2TubWe/bnmve+Mu92j391Tb782fNz36P78rknnskP952dB399U7PDiAgs5JLLX8sNd41mx48vXIK//cs1mZpKbvvll7JydZ1b73k5dZ38zdfXtGCmcNzbbzmUt958KOvOm5p3u69/aUPedstoLnvjRDZsms7tH9mfx/6k6ZtNRASW0t4nV2XbrvFUJ/xqbds1nr1PWtKiMz23Z2Ve9xPH/4G0803jeeXAWXn5xab6ICKwlMZeHciaddMnvbZm3XTGDvldozMdOTyQteuPn62cu7Hx91cPiggU+drD63LTxTtz08U7c++NFy5q33PWnh6MsUMDOefc6Tn2gPZavWY6rx4869jXr442jt+165o6Zj07C0513R0Hc90dB4v2veQN4/nKg+elns6xJa19/7Qq737/S0s4Q1g6F102kWf+YVWSxjH/3b9bnfXnT2Xj5qYi4kwEFlJPJ+NjVV6baNz2OD5WzXkL5BVvP5yBgeSR392YiSNVHn2g8Qatn772cMvmC0ky+VrjWJ2eSqanGn+ffO307d5x22i+8eiGPP3tlRndP5AvffL8XPPe0WaHcSYCC3lh74rc9ebXHfv6pot3ZtPWyTz8988kSe698cJcfuVY3vdrB7JyVbL7oefz6Q9vyRfu35wLdkxk90PPZ6Xr6rTY5z+2KX/8meO36v7Vn67PzXfvz/V3juaDV+/I733ze9l6yWTe8q7DefY7B7L7lm2ZOFLlynceys9/dH+zw3iKL32tGh55Z5IH2jT8O+qhwefaNDY9ohoeuTzJ/2zBUP9SDw2+9dQXLWcBUExE6HeTfTo2vaNVx9Gs44gI/W5Pm8YdTfJim8amt+xLcqQF43x3thdFhL4287Gf32zD0F+shwbnfxwFNGHm0waX+5rIdJIvzPYNd2dB8qE0PrXtmiTnLeM4dRqfsf7VJH+0jOPQf34jjTOSdybZnKV7PPxkGmfrX6yHBv/PbBu4OwuAYpazACgmIlCgqqq1VVWNVVX1/9o9F2gnEYEytyaZSrK9qqo3tXsy0C4iAmV+JcnaJKvSuDAPfUlEYJGqqvrJJJfMfHlWkturqlrbvhnRT6qq2lpV1VC753GUiMDi/Yc0zkCOmkpye5vmQh+oGn6mqqovp3Er7yNVVZ3d7nklbvGFRZk54/hhklM/M/3Juq4vb8OU6GFVVa1K8u+S7E6yLY3jbjLJZ+q6/nA753aUNxvC4tyaxrt3T7W9qqo31XX9d62eEL2nqqqtSX5p5s9AknNP+PZUkk+1Y16zsZwFi/Mfc/Iv9FEusHNGTlmy+l4ax9r6nHy81Ukeq+v62XbMcTaWs6BJVVW9McnjOX0p66ixJJvrun61dbOi282xZDXXY0sOJXlPXdf/uzWzW5gzEWjeh3LyBfVTucBO06qq2lBV1cfSuMb2+0l+LI3bxud77tWLSR5rwfSa5kwEmjDPBfVTucBOU6qqui6Nh3E2+7DEQ0k+Utf1f12+WS2eMxFozlwX1E/lHew068+S3J/k8CL26binP4sINGeuC+qncoGdptSNZaB7k3wmC4dkIskf1HW9mOC0hOUsWMDMBfW/TnJOk7uMJfmRuq4PLd+s6BVVVVVJPpHk7sy9XHokyY910l1ZRzkTgYVdnMZ7qiZP+XPUbK//aCsnSPeaOSPZnbkD0nG39Z7ImQgUqqqqTvKtuq5dA6HYzONLJma+/IMk/z4nB6Xjbus9kXesA7TJKQHZkeTZJC/l5KWtjrut90QiAtAGpwakruu9M6/fO/Pa3TP//XjdwUtGIgLQYnMFJGlcIzkhJO9LB97WeyLXRKCQayKUmC8gp2xXJccuvHcsZyIALdJsQJLOj8dRbvEFaIHFBKSbiAjAMuvVgCQiArCsejkgiYgALJteD0giIgDLoh8CkogIwJLrl4AkIgKwpPopIImIACyZfgtIIiIAS6IfA5KICMAZ69eAJCICcEb6OSCJiAAU6/eAJCICUERAGkQEYJEE5DgRAVgEATmZiAA0SUBOJyIATRCQ2YkIwAIEZG4iAjAPAZmfiADMQUAWJiIAsxCQ5ogIwCkEpHkiAnACAVkcEQGYISCLJyIAEZBSIgL0PQEpJyJAXxOQMyMiQN8SkDMnIkBfEpClISJA3xGQpSMiQF8RkKUlIkDfEJClJyJAXxCQ5SEiQM8TkOUjIkBPE5DlJSJAzxKQ5SciQE8SkNYQEaDnCEjriAjQUwSktUQE6BkC0noiAvQEAWkPEQG6noC0j4gAXU1A2ktEgK4lIO0nIkBXEpDOICJA1xGQziEiQFcRkM4iIkDXEJDOIyJAVxCQziQiQMcTkM4lIkBHE5DOJiJAxxKQziciQEcSkO4gIkDHEZDuISJARxGQ7iIiQMcQkO4jIkBHEJDuJCJA2wlI9xIRoK0EpLuJCNA2AtL9RARoCwHpDSICtJyA9A4RAVpKQHqLiAAtIyC9R0SAlhCQ3iQiwLITkN4lIsCyEpDeJiLAshGQ3iciwLIQkP4gIsCSE5D+ISLAkhKQ/iIiwJIRkP4jIsCSEJD+JCLAGROQ/iUiwBkRkP4mIkAxAUFEgCICQiIiQAEB4SgRARZFQDiRiABNExBOJSJAUwSE2YgIsCABYS4iAsxLQJiPiABzEhAWIiLArASEZogIcBoBoVkiApxEQFgMEQGOERAWS0SAJAJCGREBBIRiIgJ9TkA4EyICfUxAOFMiAn1KQFgKIgJ9SEBYKiICfUZAWEoiAn1EQFhqIgJ9QkBYDiICfUBAWC4iAj1OQFhOIgI9TEBYbiICPUpAaAURgR4kILSKiECPERBaSUSghwgIrSYi0CMEhHYQEegBAkK7iAh0OQGhnUQEupiA0G4iAl1KQOgEIgJdSEDoFCICXUZA6CQiAl1EQOg0IgJdQkDoRCICXUBA6FQiAh1OQOhkIgIdTEDodCICHUpA6AYiAh1IQOgWIgIdRkDoJiICHURA6DYiAh1CQOhGIgIdQEDoViICbSYgdDMRgTYSELqdiECbCAi9QESgDQSEXiEi0GICQi8REWghAaHXiAi0iIDQi0QEWkBA6FUiAstMQOhlK9o9Aegk1fDIiiT/OsmWLPSPrH/z7mT12g3V8MjNc25z5PBZ+Znr/nt+sC95bs+OeuLI3qWcL1TDI+uSvDnJ+iTVMg0zleSFJE/UQ4NTJ41f1/UyjQndpRoeeX2SzyUZbGqH5/bsyoqV49ly8bOzfr+uk+ef3pUk2XLxM1mx8ttJ3l8PDR5YmhnT76rhkfck+XiSs1s05AtJ7qyHBo8d85az4LjfSLMBWcjpAZlM8oYk9yzJz6fvzZyBtDIgSXJBkv904gsiAkmq4ZENSX5ySX7Y7AE56uolGQMaS1itDMhRPzuz7JtEROCodUvyU+YPyNKNQ8+pqmpbVVUPVFX15qqqmrm2sX7ZJzW7s5KsOfqFiEDDmV+QXDggSzMOveotSX4xyV8kebqqqrurqtowz/btPJaOjS0iMJ9HH9iYD169PTdcsDMf/4Utc25X18kf3rcrv/Ku5J5rp3P/3YOZOCIYLNZ4krVJdiT5RJLvV1X1hUWcnZxsz7dW5ldvuCi3XnZprt+8a8Htn3piVT5w1fbceNHOfOCq7XnqiVUL7SIiMJ9NWydz6z37c817X5lzm7pO/vyPduVrDyf/+Yv78rknnskP952dB399UwtnSu9Zm2R1ktvS/NnJyVacXefn3nMwd//2DxbcdmI8ue/OC3PNTa/kkT178rZbRnPfnRdmYnze3UQE5vP2Ww7lrTcfyrrzpubc5vmnd+XxryZvu/nl7LpiLBs2Tef2j+zPY3/S/C87zG0gs5yd5PHhnVnoLRqXXP5abrhrNDt+fP4SJMnf/uWaTE0lt/3yS1m5us6t97ycuk7+5utr5ttNROBMTE40Tvdf2DuRS984duz1nW8azysHzsrLL/odYykdPzv5xv/4aL7/7I4cfGljpqfO/Djb++SqbNs1nuqEH7Vt13j2PjnvkpZ3rEOpqqpT14116iOHV2Zqcmue27M1STI5c01973cuy+Do8X3GXk11/bXe4ctcmj02BlJPr87UZDK6/0cyemBzLrr0u2c08tirA1mzbvqk19asm87YoXkDJSJQavOF/5yxV9emquqsPue8TE8dzvrzG8sGrxyokgxmcOtI1p9//H8Ma84dT/Jf2jNhOtzVSf5tknOb3L5OVSXVwFTO3fBSkuRrD6/LZ3c3bgDZdcXhfOLLzzc9+jlrTw/G2KGBnHPu9Bx7JBERKLdy9XhWrm5EY9vO1dm3ZzLrz2/8Mj/1xJqsP38qF1126iNODtV1/akWz5QuUFXVvyS5tolNx5IMZMOmb2fTlh/NqjVjOXrj1nV3HMx1dxwsmsAlbxjPVx48L/V0ji1p7funVXn3+1+abzfrtTCfydeS8bEq01PJ9FTj75Ovnb7dO24bzTce3ZCnv70yo/sH8qVPnp9r3jt6+oZQpE5yKI1nV300ybb80v2/ldVrjwdk1r2mG8fsaxONjcbHqjlvPb/i7YczMJA88rsbM3GkyqMPbEyS/PS1h+ebmDMRmM/nP7Ypf/yZ47fq/tWfrs/Nd+/P9XeO5oNX78jvffN72XrJZN7yrsN59jsHsvuWbZk4UuXKdx7Kz390fxtnTm9onHUkf57kd5I8Vs88NbcaHll47xf2rshdb37dsa9vunhnNm2dzMN//0yS5N4bL8zlV47lfb92ICtXJbsfej6f/vCWfOH+zblgx0R2P/R8Vs7/VhFP8YUk1fDItiRfb8FQh+qhwX/VgnHoMlVV3ZbkwTQeKfJqkleSfCrJH9Z1/eJp2zc+guA3WzrJ466shxp3jDgTAegMI0lWJvlKkk/mhLOOTiYi0DDbc66WwywXVCCp6/ovqqpaV9f1wm8MbGjnsXRsbBfWoeGHSVrxYVFPtWAMutQiApK071jaVw8NHrvYLiKQZOYjP//bMg8zneTzyzwGfaIeGnwqyf9tw9APnfiFC+twgmp45IYk70ozn7HevMkk30vyaD00+PgS/UxINTxyTpI7klyVZEOW/zPW/1c9NPhnJ81BRAAoZTkLgGIiAkAxEQGgmIgAUExEACgmIgAUExEAiokIAMVEBIBiIgJAMREBoJiIAFBMRAAoJiIAFBMRAIqJCADFRASAYiICQDERAaCYiABQTEQAKCYiABQTEQCKiQgAxUQEgGIiAkAxEQGgmIgAUExEACgmIgAUExEAiokIAMVEBIBiIgJAMREBoJiIAFBMRAAoJiIAFBMRAIqJCADFRASAYiICQDERAaCYiABQTEQAKCYiABQTEQCKiQgAxUQEgGIiAkAxEQGgmIgAUExEACgmIgAUExEAiokIAMVEBIBiIgJAMREBoJiIAFBMRAAoJiIAFBMRAIqJCADFRASAYiICQDERAaCYiABQTEQAKCYiABQTEQCKiQgAxUQEgGIiAkAxEQGgmIgAUExEACgmIgAUExEAiokIAMVEBIBiIgJAMREBoJiIAFBMRAAoJiIAFBMRAIqJCADFRASAYiICQDERAaCYiABQTEQAKCYiABQTEQCKiQgAxUQEgGIiAkAxEQGgmIgAUExEACgmIgAUExEAiokIAMVEBIBiIgJAMREBoJiIAFBMRAAoJiIAFBMRAIqJCADFRASAYiICQDERAaCYiABQTEQAKCYiABQTEQCKiQgAxUQEgGIiAkAxEQGgmIgAUExEACgmIgAUExEAiokIAMVEBIBiIgJAMREBoJiIAFBMRAAoJiIAFBMRAIqJCADFRASAYiICQDERAaCYiABQTEQAKCYiABQTEQCKiQgAxUQEgGIiAkAxEQGgmIgAUExEACgmIgAUExEAiokIAMVEBIBiIgJAMREBoJiIAFBMRAAoJiIAFBMRAIqJCADFRASAYiICQDERAaCYiABQTEQAKCYiABQTEQCKiQgAxUQEgGIiAkAxEQGgmIgAUExEACgmIgAUExEAiokIAMX+P2nvVhY2t2ylAAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ "<matplotlib.figure.Figure at 0x7f06c669c898>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "d = FiniteDifferenceStaggeredStencilDerivation(\"NE\", 2, (0, 1))\n",
+ "assert d.apply(c) == c[0,0] + c[1,1] - c[1,0] - c[0,1]\n",
+ "d.visualize()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "v3 = ps.fields(\"v(3): [3D]\")\n",
+ "assert FiniteDifferenceStaggeredStencilDerivation(\"E\", 3, (0,)).apply(v3) == \\\n",
+ " sp.Matrix([v3[1,0,0](i) - v3[0,0,0](i) for i in range(*v3.index_shape)])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
}
],
"metadata": {
@@ -193,7 +395,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.3"
+ "version": "3.6.9"
}
},
"nbformat": 4,