from pypol import *
-x, y = symbols('x y')
+x, y, z = symbols('x y z')
diam = Ge(y, x - 1) & Le(y, x + 1) & Ge(y, -x - 1) & Le(y, -x + 1)
print('diamond:', diam)
-print('projected on x:', diam.drop_dims('y'))
+print()
+rhom1 = Le(0, x) & Le(x, 3) & Le(0, y) & Le(y, 3) & Le(0, z) & Le(z, 3) \
+& Le(z - 2, x) & Ge(z + 2, x) & Ge(z - 1, -x) & Le(z - 5, -x) \
+& Le(z - 2, y) & Ge(z + 2, y) & Ge(z - 1, -y) & Le(z - 5, -y) \
+& Le(y - 2, x) & Ge(y + 2, x) & Ge(y - 1, -x) & Le(y - 5, -x)
+rhom1.plot()
+rhom2 = rhom1 & Le(x + y + z, 7) & Ge(-2, -x - y - z ) \
+& Le(x + y - z, 4) & Ge(x + y - z, -1) \
+& Le(x - y + z, 4) & Ge(x - y + z, -1) \
+& Le(-x + y + z, 4) & Ge(-x + y + z, -1)
+rhom2.plot()