4 find_containing_tetrahedron,
9 import Data.Maybe ( fromJust )
10 import qualified Data.Vector as V (
18 import Prelude hiding ( LT )
19 import Test.Framework ( Test, testGroup )
20 import Test.Framework.Providers.QuickCheck2 ( testProperty )
21 import Test.QuickCheck ( Arbitrary(..), Gen, Positive(..), choose )
31 import Comparisons ( (~=), (~~=) )
32 import qualified Face ( Face(..), center )
33 import FunctionValues ( FunctionValues, eval, rotate )
34 import Misc ( all_equal, disjoint )
35 import Point ( Point(..), dot )
36 import Tetrahedron ( Tetrahedron(..), barycenter, c, volume )
38 data Cube = Cube { i :: !Int,
41 fv :: !FunctionValues,
42 tetrahedra_volume :: !Double }
46 instance Arbitrary Cube where
48 i' <- choose (coordmin, coordmax)
49 j' <- choose (coordmin, coordmax)
50 k' <- choose (coordmin, coordmax)
51 fv' <- arbitrary :: Gen FunctionValues
52 (Positive tet_vol) <- arbitrary :: Gen (Positive Double)
53 return (Cube i' j' k' fv' tet_vol)
55 -- The idea here is that, when cubed in the volume formula,
56 -- these numbers don't overflow 64 bits. This number is not
57 -- magic in any other sense than that it does not cause test
58 -- failures, while 2^23 does.
59 coordmax = 4194304 -- 2^22
63 instance Show Cube where
65 "Cube_" ++ subscript ++ "\n" ++
66 " Center: " ++ (show (center cube)) ++ "\n" ++
67 " xmin: " ++ (show (xmin cube)) ++ "\n" ++
68 " xmax: " ++ (show (xmax cube)) ++ "\n" ++
69 " ymin: " ++ (show (ymin cube)) ++ "\n" ++
70 " ymax: " ++ (show (ymax cube)) ++ "\n" ++
71 " zmin: " ++ (show (zmin cube)) ++ "\n" ++
72 " zmax: " ++ (show (zmax cube)) ++ "\n"
75 (show (i cube)) ++ "," ++ (show (j cube)) ++ "," ++ (show (k cube))
78 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
80 xmin :: Cube -> Double
81 xmin cube = (i' - 1/2)
83 i' = fromIntegral (i cube) :: Double
85 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
87 xmax :: Cube -> Double
88 xmax cube = (i' + 1/2)
90 i' = fromIntegral (i cube) :: Double
92 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
94 ymin :: Cube -> Double
95 ymin cube = (j' - 1/2)
97 j' = fromIntegral (j cube) :: Double
99 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
101 ymax :: Cube -> Double
102 ymax cube = (j' + 1/2)
104 j' = fromIntegral (j cube) :: Double
106 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
108 zmin :: Cube -> Double
109 zmin cube = (k' - 1/2)
111 k' = fromIntegral (k cube) :: Double
113 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
115 zmax :: Cube -> Double
116 zmax cube = (k' + 1/2)
118 k' = fromIntegral (k cube) :: Double
121 -- | The center of Cube_ijk coincides with v_ijk at
122 -- (i, j, k). See Sorokina and Zeilfelder, p. 76.
123 center :: Cube -> Point
127 x = fromIntegral (i cube) :: Double
128 y = fromIntegral (j cube) :: Double
129 z = fromIntegral (k cube) :: Double
134 -- | The top (in the direction of z) face of the cube.
135 top_face :: Cube -> Face.Face
136 top_face cube = Face.Face v0' v1' v2' v3'
140 v0' = cc + ( Point delta (-delta) delta )
141 v1' = cc + ( Point delta delta delta )
142 v2' = cc + ( Point (-delta) delta delta )
143 v3' = cc + ( Point (-delta) (-delta) delta )
147 -- | The back (in the direction of x) face of the cube.
148 back_face :: Cube -> Face.Face
149 back_face cube = Face.Face v0' v1' v2' v3'
153 v0' = cc + ( Point delta (-delta) (-delta) )
154 v1' = cc + ( Point delta delta (-delta) )
155 v2' = cc + ( Point delta delta delta )
156 v3' = cc + ( Point delta (-delta) delta )
159 -- The bottom face (in the direction of -z) of the cube.
160 down_face :: Cube -> Face.Face
161 down_face cube = Face.Face v0' v1' v2' v3'
165 v0' = cc + ( Point (-delta) (-delta) (-delta) )
166 v1' = cc + ( Point (-delta) delta (-delta) )
167 v2' = cc + ( Point delta delta (-delta) )
168 v3' = cc + ( Point delta (-delta) (-delta) )
172 -- | The front (in the direction of -x) face of the cube.
173 front_face :: Cube -> Face.Face
174 front_face cube = Face.Face v0' v1' v2' v3'
178 v0' = cc + ( Point (-delta) (-delta) delta )
179 v1' = cc + ( Point (-delta) delta delta )
180 v2' = cc + ( Point (-delta) delta (-delta) )
181 v3' = cc + ( Point (-delta) (-delta) (-delta) )
183 -- | The left (in the direction of -y) face of the cube.
184 left_face :: Cube -> Face.Face
185 left_face cube = Face.Face v0' v1' v2' v3'
189 v0' = cc + ( Point delta (-delta) delta )
190 v1' = cc + ( Point (-delta) (-delta) delta )
191 v2' = cc + ( Point (-delta) (-delta) (-delta) )
192 v3' = cc + ( Point delta (-delta) (-delta) )
195 -- | The right (in the direction of y) face of the cube.
196 right_face :: Cube -> Face.Face
197 right_face cube = Face.Face v0' v1' v2' v3'
201 v0' = cc + ( Point (-delta) delta delta)
202 v1' = cc + ( Point delta delta delta )
203 v2' = cc + ( Point delta delta (-delta) )
204 v3' = cc + ( Point (-delta) delta (-delta) )
207 tetrahedron :: Cube -> Int -> Tetrahedron
210 Tetrahedron (fv cube) v0' v1' v2' v3' vol
217 vol = tetrahedra_volume cube
220 Tetrahedron fv' v0' v1' v2' v3' vol
227 fv' = rotate ccwx (fv cube)
228 vol = tetrahedra_volume cube
231 Tetrahedron fv' v0' v1' v2' v3' vol
238 fv' = rotate ccwx $ rotate ccwx $ fv cube
239 vol = tetrahedra_volume cube
242 Tetrahedron fv' v0' v1' v2' v3' vol
249 fv' = rotate cwx (fv cube)
250 vol = tetrahedra_volume cube
253 Tetrahedron fv' v0' v1' v2' v3' vol
260 fv' = rotate cwy (fv cube)
261 vol = tetrahedra_volume cube
264 Tetrahedron fv' v0' v1' v2' v3' vol
271 fv' = rotate cwy $ rotate cwz $ fv cube
272 vol = tetrahedra_volume cube
275 Tetrahedron fv' v0' v1' v2' v3' vol
282 fv' = rotate cwy $ rotate cwz
285 vol = tetrahedra_volume cube
288 Tetrahedron fv' v0' v1' v2' v3' vol
295 fv' = rotate cwy $ rotate ccwz $ fv cube
296 vol = tetrahedra_volume cube
299 Tetrahedron fv' v0' v1' v2' v3' vol
306 fv' = rotate cwy $ rotate cwy $ fv cube
307 vol = tetrahedra_volume cube
310 Tetrahedron fv' v0' v1' v2' v3' vol
317 fv' = rotate cwy $ rotate cwy
320 vol = tetrahedra_volume cube
322 tetrahedron cube 10 =
323 Tetrahedron fv' v0' v1' v2' v3' vol
330 fv' = rotate cwy $ rotate cwy
335 vol = tetrahedra_volume cube
337 tetrahedron cube 11 =
338 Tetrahedron fv' v0' v1' v2' v3' vol
345 fv' = rotate cwy $ rotate cwy
348 vol = tetrahedra_volume cube
350 tetrahedron cube 12 =
351 Tetrahedron fv' v0' v1' v2' v3' vol
358 fv' = rotate ccwy $ fv cube
359 vol = tetrahedra_volume cube
361 tetrahedron cube 13 =
362 Tetrahedron fv' v0' v1' v2' v3' vol
369 fv' = rotate ccwy $ rotate ccwz $ fv cube
370 vol = tetrahedra_volume cube
372 tetrahedron cube 14 =
373 Tetrahedron fv' v0' v1' v2' v3' vol
380 fv' = rotate ccwy $ rotate ccwz
383 vol = tetrahedra_volume cube
385 tetrahedron cube 15 =
386 Tetrahedron fv' v0' v1' v2' v3' vol
393 fv' = rotate ccwy $ rotate cwz $ fv cube
394 vol = tetrahedra_volume cube
396 tetrahedron cube 16 =
397 Tetrahedron fv' v0' v1' v2' v3' vol
404 fv' = rotate ccwz $ fv cube
405 vol = tetrahedra_volume cube
407 tetrahedron cube 17 =
408 Tetrahedron fv' v0' v1' v2' v3' vol
415 fv' = rotate ccwz $ rotate cwy $ fv cube
416 vol = tetrahedra_volume cube
418 tetrahedron cube 18 =
419 Tetrahedron fv' v0' v1' v2' v3' vol
426 fv' = rotate ccwz $ rotate cwy
429 vol = tetrahedra_volume cube
431 tetrahedron cube 19 =
432 Tetrahedron fv' v0' v1' v2' v3' vol
439 fv' = rotate ccwz $ rotate ccwy
441 vol = tetrahedra_volume cube
443 tetrahedron cube 20 =
444 Tetrahedron fv' v0' v1' v2' v3' vol
451 fv' = rotate cwz $ fv cube
452 vol = tetrahedra_volume cube
454 tetrahedron cube 21 =
455 Tetrahedron fv' v0' v1' v2' v3' vol
462 fv' = rotate cwz $ rotate ccwy $ fv cube
463 vol = tetrahedra_volume cube
465 tetrahedron cube 22 =
466 Tetrahedron fv' v0' v1' v2' v3' vol
473 fv' = rotate cwz $ rotate ccwy
476 vol = tetrahedra_volume cube
478 tetrahedron cube 23 =
479 Tetrahedron fv' v0' v1' v2' v3' vol
486 fv' = rotate cwz $ rotate cwy
488 vol = tetrahedra_volume cube
491 -- Only used in tests, so we don't need the added speed
493 tetrahedra :: Cube -> [Tetrahedron]
494 tetrahedra cube = [ tetrahedron cube n | n <- [0..23] ]
496 front_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
497 front_left_top_tetrahedra cube =
498 V.singleton (tetrahedron cube 0) `V.snoc`
499 (tetrahedron cube 3) `V.snoc`
500 (tetrahedron cube 6) `V.snoc`
501 (tetrahedron cube 7) `V.snoc`
502 (tetrahedron cube 20) `V.snoc`
503 (tetrahedron cube 21)
505 front_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
506 front_left_down_tetrahedra cube =
507 V.singleton (tetrahedron cube 0) `V.snoc`
508 (tetrahedron cube 2) `V.snoc`
509 (tetrahedron cube 3) `V.snoc`
510 (tetrahedron cube 12) `V.snoc`
511 (tetrahedron cube 15) `V.snoc`
512 (tetrahedron cube 21)
514 front_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
515 front_right_top_tetrahedra cube =
516 V.singleton (tetrahedron cube 0) `V.snoc`
517 (tetrahedron cube 1) `V.snoc`
518 (tetrahedron cube 5) `V.snoc`
519 (tetrahedron cube 6) `V.snoc`
520 (tetrahedron cube 16) `V.snoc`
521 (tetrahedron cube 19)
523 front_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
524 front_right_down_tetrahedra cube =
525 V.singleton (tetrahedron cube 1) `V.snoc`
526 (tetrahedron cube 2) `V.snoc`
527 (tetrahedron cube 12) `V.snoc`
528 (tetrahedron cube 13) `V.snoc`
529 (tetrahedron cube 18) `V.snoc`
530 (tetrahedron cube 19)
532 back_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
533 back_left_top_tetrahedra cube =
534 V.singleton (tetrahedron cube 0) `V.snoc`
535 (tetrahedron cube 3) `V.snoc`
536 (tetrahedron cube 6) `V.snoc`
537 (tetrahedron cube 7) `V.snoc`
538 (tetrahedron cube 20) `V.snoc`
539 (tetrahedron cube 21)
541 back_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
542 back_left_down_tetrahedra cube =
543 V.singleton (tetrahedron cube 8) `V.snoc`
544 (tetrahedron cube 11) `V.snoc`
545 (tetrahedron cube 14) `V.snoc`
546 (tetrahedron cube 15) `V.snoc`
547 (tetrahedron cube 22) `V.snoc`
548 (tetrahedron cube 23)
550 back_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
551 back_right_top_tetrahedra cube =
552 V.singleton (tetrahedron cube 4) `V.snoc`
553 (tetrahedron cube 5) `V.snoc`
554 (tetrahedron cube 9) `V.snoc`
555 (tetrahedron cube 10) `V.snoc`
556 (tetrahedron cube 16) `V.snoc`
557 (tetrahedron cube 17)
559 back_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
560 back_right_down_tetrahedra cube =
561 V.singleton (tetrahedron cube 8) `V.snoc`
562 (tetrahedron cube 9) `V.snoc`
563 (tetrahedron cube 13) `V.snoc`
564 (tetrahedron cube 14) `V.snoc`
565 (tetrahedron cube 17) `V.snoc`
566 (tetrahedron cube 18)
568 in_top_half :: Cube -> Point -> Bool
569 in_top_half cube (Point _ _ z) =
570 distance_from_top <= distance_from_bottom
572 distance_from_top = abs $ (zmax cube) - z
573 distance_from_bottom = abs $ (zmin cube) - z
575 in_front_half :: Cube -> Point -> Bool
576 in_front_half cube (Point x _ _) =
577 distance_from_front <= distance_from_back
579 distance_from_front = abs $ (xmin cube) - x
580 distance_from_back = abs $ (xmax cube) - x
583 in_left_half :: Cube -> Point -> Bool
584 in_left_half cube (Point _ y _) =
585 distance_from_left <= distance_from_right
587 distance_from_left = abs $ (ymin cube) - y
588 distance_from_right = abs $ (ymax cube) - y
591 -- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
592 -- contain the given 'Point'. This should be faster than checking
593 -- every tetrahedron individually, since we determine which half
594 -- (hemisphere?) of the cube the point lies in three times: once in
595 -- each dimension. This allows us to eliminate non-candidates
598 -- This can throw an exception, but the use of 'head' might
599 -- save us some unnecessary computations.
601 {-# INLINE find_containing_tetrahedron #-}
602 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
603 find_containing_tetrahedron cube p =
604 candidates `V.unsafeIndex` (fromJust lucky_idx)
606 front_half = in_front_half cube p
607 top_half = in_top_half cube p
608 left_half = in_left_half cube p
610 candidates :: V.Vector Tetrahedron
615 front_left_top_tetrahedra cube
617 front_left_down_tetrahedra cube
620 front_right_top_tetrahedra cube
622 front_right_down_tetrahedra cube
624 | otherwise = -- back half
627 back_left_top_tetrahedra cube
629 back_left_down_tetrahedra cube
632 back_right_top_tetrahedra cube
634 back_right_down_tetrahedra cube
636 -- Use the dot product instead of Euclidean distance here to save
637 -- a sqrt(). So, "distances" below really means "distances
639 distances :: V.Vector Double
640 distances = V.map ((dot p) . barycenter) candidates
642 shortest_distance :: Double
643 shortest_distance = V.minimum distances
645 -- Compute the index of the tetrahedron with the center closest to
646 -- p. This is a bad algorithm, but don't change it! If you make it
647 -- smarter by finding the index of shortest_distance in distances
648 -- (this should give the same answer and avoids recomputing the
649 -- dot product), the program gets slower. Seriously!
650 lucky_idx :: Maybe Int
651 lucky_idx = V.findIndex
652 (\t -> (barycenter t) `dot` p == shortest_distance)
664 prop_opposite_octant_tetrahedra_disjoint1 :: Cube -> Bool
665 prop_opposite_octant_tetrahedra_disjoint1 cube =
666 disjoint (front_left_top_tetrahedra cube) (front_right_down_tetrahedra cube)
668 prop_opposite_octant_tetrahedra_disjoint2 :: Cube -> Bool
669 prop_opposite_octant_tetrahedra_disjoint2 cube =
670 disjoint (back_left_top_tetrahedra cube) (back_right_down_tetrahedra cube)
672 prop_opposite_octant_tetrahedra_disjoint3 :: Cube -> Bool
673 prop_opposite_octant_tetrahedra_disjoint3 cube =
674 disjoint (front_left_top_tetrahedra cube) (back_right_top_tetrahedra cube)
676 prop_opposite_octant_tetrahedra_disjoint4 :: Cube -> Bool
677 prop_opposite_octant_tetrahedra_disjoint4 cube =
678 disjoint (front_left_down_tetrahedra cube) (back_right_down_tetrahedra cube)
680 prop_opposite_octant_tetrahedra_disjoint5 :: Cube -> Bool
681 prop_opposite_octant_tetrahedra_disjoint5 cube =
682 disjoint (front_left_top_tetrahedra cube) (back_left_down_tetrahedra cube)
684 prop_opposite_octant_tetrahedra_disjoint6 :: Cube -> Bool
685 prop_opposite_octant_tetrahedra_disjoint6 cube =
686 disjoint (front_right_top_tetrahedra cube) (back_right_down_tetrahedra cube)
689 -- | Since the grid size is necessarily positive, all tetrahedra
690 -- (which comprise cubes of positive volume) must have positive
692 prop_all_volumes_positive :: Cube -> Bool
693 prop_all_volumes_positive cube =
697 volumes = map volume ts
700 -- | In fact, since all of the tetrahedra are identical, we should
701 -- already know their volumes. There's 24 tetrahedra to a cube, so
702 -- we'd expect the volume of each one to be 1/24.
703 prop_all_volumes_exact :: Cube -> Bool
704 prop_all_volumes_exact cube =
705 and [volume t ~~= 1/24 | t <- tetrahedra cube]
707 -- | All tetrahedron should have their v0 located at the center of the
709 prop_v0_all_equal :: Cube -> Bool
710 prop_v0_all_equal cube = (v0 t0) == (v0 t1)
712 t0 = head (tetrahedra cube) -- Doesn't matter which two we choose.
713 t1 = head $ tail (tetrahedra cube)
716 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Note that the
717 -- third and fourth indices of c-t3 have been switched. This is
718 -- because we store the triangles oriented such that their volume is
719 -- positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde point
720 -- in opposite directions, one of them has to have negative volume!
721 prop_c0120_identity1 :: Cube -> Bool
722 prop_c0120_identity1 cube =
723 c t0 0 1 2 0 ~= (c t0 0 0 2 1 + c t3 0 0 1 2) / 2
725 t0 = tetrahedron cube 0
726 t3 = tetrahedron cube 3
729 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
730 -- 'prop_c0120_identity1' with tetrahedrons 1 and 2.
731 prop_c0120_identity2 :: Cube -> Bool
732 prop_c0120_identity2 cube =
733 c t1 0 1 2 0 ~= (c t1 0 0 2 1 + c t0 0 0 1 2) / 2
735 t0 = tetrahedron cube 0
736 t1 = tetrahedron cube 1
738 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
739 -- 'prop_c0120_identity1' with tetrahedrons 1 and 2.
740 prop_c0120_identity3 :: Cube -> Bool
741 prop_c0120_identity3 cube =
742 c t2 0 1 2 0 ~= (c t2 0 0 2 1 + c t1 0 0 1 2) / 2
744 t1 = tetrahedron cube 1
745 t2 = tetrahedron cube 2
747 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
748 -- 'prop_c0120_identity1' with tetrahedrons 2 and 3.
749 prop_c0120_identity4 :: Cube -> Bool
750 prop_c0120_identity4 cube =
751 c t3 0 1 2 0 ~= (c t3 0 0 2 1 + c t2 0 0 1 2) / 2
753 t2 = tetrahedron cube 2
754 t3 = tetrahedron cube 3
757 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
758 -- 'prop_c0120_identity1' with tetrahedrons 4 and 5.
759 prop_c0120_identity5 :: Cube -> Bool
760 prop_c0120_identity5 cube =
761 c t5 0 1 2 0 ~= (c t5 0 0 2 1 + c t4 0 0 1 2) / 2
763 t4 = tetrahedron cube 4
764 t5 = tetrahedron cube 5
766 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
767 -- 'prop_c0120_identity1' with tetrahedrons 5 and 6.
768 prop_c0120_identity6 :: Cube -> Bool
769 prop_c0120_identity6 cube =
770 c t6 0 1 2 0 ~= (c t6 0 0 2 1 + c t5 0 0 1 2) / 2
772 t5 = tetrahedron cube 5
773 t6 = tetrahedron cube 6
776 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
777 -- 'prop_c0120_identity1' with tetrahedrons 6 and 7.
778 prop_c0120_identity7 :: Cube -> Bool
779 prop_c0120_identity7 cube =
780 c t7 0 1 2 0 ~= (c t7 0 0 2 1 + c t6 0 0 1 2) / 2
782 t6 = tetrahedron cube 6
783 t7 = tetrahedron cube 7
786 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
787 -- 'prop_c0120_identity1'.
788 prop_c0210_identity1 :: Cube -> Bool
789 prop_c0210_identity1 cube =
790 c t0 0 2 1 0 ~= (c t0 0 1 1 1 + c t3 0 1 1 1) / 2
792 t0 = tetrahedron cube 0
793 t3 = tetrahedron cube 3
796 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
797 -- 'prop_c0120_identity1'.
798 prop_c0300_identity1 :: Cube -> Bool
799 prop_c0300_identity1 cube =
800 c t0 0 3 0 0 ~= (c t0 0 2 0 1 + c t3 0 2 1 0) / 2
802 t0 = tetrahedron cube 0
803 t3 = tetrahedron cube 3
806 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
807 -- 'prop_c0120_identity1'.
808 prop_c1110_identity :: Cube -> Bool
809 prop_c1110_identity cube =
810 c t0 1 1 1 0 ~= (c t0 1 0 1 1 + c t3 1 0 1 1) / 2
812 t0 = tetrahedron cube 0
813 t3 = tetrahedron cube 3
816 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
817 -- 'prop_c0120_identity1'.
818 prop_c1200_identity1 :: Cube -> Bool
819 prop_c1200_identity1 cube =
820 c t0 1 2 0 0 ~= (c t0 1 1 0 1 + c t3 1 1 1 0) / 2
822 t0 = tetrahedron cube 0
823 t3 = tetrahedron cube 3
826 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
827 -- 'prop_c0120_identity1'.
828 prop_c2100_identity1 :: Cube -> Bool
829 prop_c2100_identity1 cube =
830 c t0 2 1 0 0 ~= (c t0 2 0 0 1 + c t3 2 0 1 0) / 2
832 t0 = tetrahedron cube 0
833 t3 = tetrahedron cube 3
837 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). Note that the
838 -- third and fourth indices of c-t3 have been switched. This is
839 -- because we store the triangles oriented such that their volume is
840 -- positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde
841 -- point in opposite directions, one of them has to have negative
843 prop_c0102_identity1 :: Cube -> Bool
844 prop_c0102_identity1 cube =
845 c t0 0 1 0 2 ~= (c t0 0 0 1 2 + c t1 0 0 2 1) / 2
847 t0 = tetrahedron cube 0
848 t1 = tetrahedron cube 1
851 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
852 -- 'prop_c0102_identity1'.
853 prop_c0201_identity1 :: Cube -> Bool
854 prop_c0201_identity1 cube =
855 c t0 0 2 0 1 ~= (c t0 0 1 1 1 + c t1 0 1 1 1) / 2
857 t0 = tetrahedron cube 0
858 t1 = tetrahedron cube 1
861 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
862 -- 'prop_c0102_identity1'.
863 prop_c0300_identity2 :: Cube -> Bool
864 prop_c0300_identity2 cube =
865 c t0 0 3 0 0 ~= (c t0 0 2 1 0 + c t1 0 2 0 1) / 2
867 t0 = tetrahedron cube 0
868 t1 = tetrahedron cube 1
871 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
872 -- 'prop_c0102_identity1'.
873 prop_c1101_identity :: Cube -> Bool
874 prop_c1101_identity cube =
875 c t0 1 1 0 1 ~= (c t0 1 0 1 1 + c t1 1 0 1 1) / 2
877 t0 = tetrahedron cube 0
878 t1 = tetrahedron cube 1
881 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
882 -- 'prop_c0102_identity1'.
883 prop_c1200_identity2 :: Cube -> Bool
884 prop_c1200_identity2 cube =
885 c t0 1 2 0 0 ~= (c t0 1 1 1 0 + c t1 1 1 0 1) / 2
887 t0 = tetrahedron cube 0
888 t1 = tetrahedron cube 1
891 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
892 -- 'prop_c0102_identity1'.
893 prop_c2100_identity2 :: Cube -> Bool
894 prop_c2100_identity2 cube =
895 c t0 2 1 0 0 ~= (c t0 2 0 1 0 + c t1 2 0 0 1) / 2
897 t0 = tetrahedron cube 0
898 t1 = tetrahedron cube 1
901 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). The third and
902 -- fourth indices of c-t6 have been switched. This is because we
903 -- store the triangles oriented such that their volume is
904 -- positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde
905 -- point in opposite directions, one of them has to have negative
907 prop_c3000_identity :: Cube -> Bool
908 prop_c3000_identity cube =
909 c t0 3 0 0 0 ~= c t0 2 1 0 0 + c t6 2 1 0 0
910 - ((c t0 2 0 1 0 + c t0 2 0 0 1)/ 2)
912 t0 = tetrahedron cube 0
913 t6 = tetrahedron cube 6
916 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
917 -- 'prop_c3000_identity'.
918 prop_c2010_identity :: Cube -> Bool
919 prop_c2010_identity cube =
920 c t0 2 0 1 0 ~= c t0 1 1 1 0 + c t6 1 1 0 1
921 - ((c t0 1 0 2 0 + c t0 1 0 1 1)/ 2)
923 t0 = tetrahedron cube 0
924 t6 = tetrahedron cube 6
927 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
928 -- 'prop_c3000_identity'.
929 prop_c2001_identity :: Cube -> Bool
930 prop_c2001_identity cube =
931 c t0 2 0 0 1 ~= c t0 1 1 0 1 + c t6 1 1 1 0
932 - ((c t0 1 0 0 2 + c t0 1 0 1 1)/ 2)
934 t0 = tetrahedron cube 0
935 t6 = tetrahedron cube 6
938 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
939 -- 'prop_c3000_identity'.
940 prop_c1020_identity :: Cube -> Bool
941 prop_c1020_identity cube =
942 c t0 1 0 2 0 ~= c t0 0 1 2 0 + c t6 0 1 0 2
943 - ((c t0 0 0 3 0 + c t0 0 0 2 1)/ 2)
945 t0 = tetrahedron cube 0
946 t6 = tetrahedron cube 6
949 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
950 -- 'prop_c3000_identity'.
951 prop_c1002_identity :: Cube -> Bool
952 prop_c1002_identity cube =
953 c t0 1 0 0 2 ~= c t0 0 1 0 2 + c t6 0 1 2 0
954 - ((c t0 0 0 0 3 + c t0 0 0 1 2)/ 2)
956 t0 = tetrahedron cube 0
957 t6 = tetrahedron cube 6
960 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
961 -- 'prop_c3000_identity'.
962 prop_c1011_identity :: Cube -> Bool
963 prop_c1011_identity cube =
964 c t0 1 0 1 1 ~= c t0 0 1 1 1 + c t6 0 1 1 1 -
965 ((c t0 0 0 1 2 + c t0 0 0 2 1)/ 2)
967 t0 = tetrahedron cube 0
968 t6 = tetrahedron cube 6
971 -- | The function values at the interior should be the same for all
973 prop_interior_values_all_identical :: Cube -> Bool
974 prop_interior_values_all_identical cube =
975 all_equal [ eval (function_values tet) I | tet <- tetrahedra cube ]
978 -- | We know what (c t6 2 1 0 0) should be from Sorokina and Zeilfelder, p. 87.
979 -- This test checks the rotation works as expected.
980 prop_c_tilde_2100_rotation_correct :: Cube -> Bool
981 prop_c_tilde_2100_rotation_correct cube =
984 t0 = tetrahedron cube 0
985 t6 = tetrahedron cube 6
987 -- What gets computed for c2100 of t6.
988 expr1 = eval (function_values t6) $
990 (1/12)*(T + R + L + D) +
991 (1/64)*(FT + FR + FL + FD) +
994 (1/96)*(RT + LD + LT + RD) +
995 (1/192)*(BT + BR + BL + BD)
997 -- What should be computed for c2100 of t6.
998 expr2 = eval (function_values t0) $
1000 (1/12)*(F + R + L + B) +
1001 (1/64)*(FT + RT + LT + BT) +
1004 (1/96)*(FR + FL + BR + BL) +
1005 (1/192)*(FD + RD + LD + BD)
1008 -- | We know what (c t6 2 1 0 0) should be from Sorokina and
1009 -- Zeilfelder, p. 87. This test checks the actual value based on
1010 -- the FunctionValues of the cube.
1012 -- If 'prop_c_tilde_2100_rotation_correct' passes, then this test is
1014 prop_c_tilde_2100_correct :: Cube -> Bool
1015 prop_c_tilde_2100_correct cube =
1016 c t6 2 1 0 0 ~= expected
1018 t0 = tetrahedron cube 0
1019 t6 = tetrahedron cube 6
1020 fvs = function_values t0
1021 expected = eval fvs $
1023 (1/12)*(F + R + L + B) +
1024 (1/64)*(FT + RT + LT + BT) +
1027 (1/96)*(FR + FL + BR + BL) +
1028 (1/192)*(FD + RD + LD + BD)
1031 -- Tests to check that the correct edges are incidental.
1032 prop_t0_shares_edge_with_t1 :: Cube -> Bool
1033 prop_t0_shares_edge_with_t1 cube =
1034 (v1 t0) == (v1 t1) && (v3 t0) == (v2 t1)
1036 t0 = tetrahedron cube 0
1037 t1 = tetrahedron cube 1
1039 prop_t0_shares_edge_with_t3 :: Cube -> Bool
1040 prop_t0_shares_edge_with_t3 cube =
1041 (v1 t0) == (v1 t3) && (v2 t0) == (v3 t3)
1043 t0 = tetrahedron cube 0
1044 t3 = tetrahedron cube 3
1046 prop_t0_shares_edge_with_t6 :: Cube -> Bool
1047 prop_t0_shares_edge_with_t6 cube =
1048 (v2 t0) == (v3 t6) && (v3 t0) == (v2 t6)
1050 t0 = tetrahedron cube 0
1051 t6 = tetrahedron cube 6
1053 prop_t1_shares_edge_with_t2 :: Cube -> Bool
1054 prop_t1_shares_edge_with_t2 cube =
1055 (v1 t1) == (v1 t2) && (v3 t1) == (v2 t2)
1057 t1 = tetrahedron cube 1
1058 t2 = tetrahedron cube 2
1060 prop_t1_shares_edge_with_t19 :: Cube -> Bool
1061 prop_t1_shares_edge_with_t19 cube =
1062 (v2 t1) == (v3 t19) && (v3 t1) == (v2 t19)
1064 t1 = tetrahedron cube 1
1065 t19 = tetrahedron cube 19
1067 prop_t2_shares_edge_with_t3 :: Cube -> Bool
1068 prop_t2_shares_edge_with_t3 cube =
1069 (v1 t1) == (v1 t2) && (v3 t1) == (v2 t2)
1071 t1 = tetrahedron cube 1
1072 t2 = tetrahedron cube 2
1074 prop_t2_shares_edge_with_t12 :: Cube -> Bool
1075 prop_t2_shares_edge_with_t12 cube =
1076 (v2 t2) == (v3 t12) && (v3 t2) == (v2 t12)
1078 t2 = tetrahedron cube 2
1079 t12 = tetrahedron cube 12
1081 prop_t3_shares_edge_with_t21 :: Cube -> Bool
1082 prop_t3_shares_edge_with_t21 cube =
1083 (v2 t3) == (v3 t21) && (v3 t3) == (v2 t21)
1085 t3 = tetrahedron cube 3
1086 t21 = tetrahedron cube 21
1088 prop_t4_shares_edge_with_t5 :: Cube -> Bool
1089 prop_t4_shares_edge_with_t5 cube =
1090 (v1 t4) == (v1 t5) && (v3 t4) == (v2 t5)
1092 t4 = tetrahedron cube 4
1093 t5 = tetrahedron cube 5
1095 prop_t4_shares_edge_with_t7 :: Cube -> Bool
1096 prop_t4_shares_edge_with_t7 cube =
1097 (v1 t4) == (v1 t7) && (v2 t4) == (v3 t7)
1099 t4 = tetrahedron cube 4
1100 t7 = tetrahedron cube 7
1102 prop_t4_shares_edge_with_t10 :: Cube -> Bool
1103 prop_t4_shares_edge_with_t10 cube =
1104 (v2 t4) == (v3 t10) && (v3 t4) == (v2 t10)
1106 t4 = tetrahedron cube 4
1107 t10 = tetrahedron cube 10
1109 prop_t5_shares_edge_with_t6 :: Cube -> Bool
1110 prop_t5_shares_edge_with_t6 cube =
1111 (v1 t5) == (v1 t6) && (v3 t5) == (v2 t6)
1113 t5 = tetrahedron cube 5
1114 t6 = tetrahedron cube 6
1116 prop_t5_shares_edge_with_t16 :: Cube -> Bool
1117 prop_t5_shares_edge_with_t16 cube =
1118 (v2 t5) == (v3 t16) && (v3 t5) == (v2 t16)
1120 t5 = tetrahedron cube 5
1121 t16 = tetrahedron cube 16
1123 prop_t6_shares_edge_with_t7 :: Cube -> Bool
1124 prop_t6_shares_edge_with_t7 cube =
1125 (v1 t6) == (v1 t7) && (v3 t6) == (v2 t7)
1127 t6 = tetrahedron cube 6
1128 t7 = tetrahedron cube 7
1130 prop_t7_shares_edge_with_t20 :: Cube -> Bool
1131 prop_t7_shares_edge_with_t20 cube =
1132 (v2 t7) == (v3 t20) && (v2 t7) == (v3 t20)
1134 t7 = tetrahedron cube 7
1135 t20 = tetrahedron cube 20
1138 p79_26_properties :: Test.Framework.Test
1140 testGroup "p. 79, Section (2.6) Properties" [
1141 testProperty "c0120 identity1" prop_c0120_identity1,
1142 testProperty "c0120 identity2" prop_c0120_identity2,
1143 testProperty "c0120 identity3" prop_c0120_identity3,
1144 testProperty "c0120 identity4" prop_c0120_identity4,
1145 testProperty "c0120 identity5" prop_c0120_identity5,
1146 testProperty "c0120 identity6" prop_c0120_identity6,
1147 testProperty "c0120 identity7" prop_c0120_identity7,
1148 testProperty "c0210 identity1" prop_c0210_identity1,
1149 testProperty "c0300 identity1" prop_c0300_identity1,
1150 testProperty "c1110 identity" prop_c1110_identity,
1151 testProperty "c1200 identity1" prop_c1200_identity1,
1152 testProperty "c2100 identity1" prop_c2100_identity1]
1154 p79_27_properties :: Test.Framework.Test
1156 testGroup "p. 79, Section (2.7) Properties" [
1157 testProperty "c0102 identity1" prop_c0102_identity1,
1158 testProperty "c0201 identity1" prop_c0201_identity1,
1159 testProperty "c0300 identity2" prop_c0300_identity2,
1160 testProperty "c1101 identity" prop_c1101_identity,
1161 testProperty "c1200 identity2" prop_c1200_identity2,
1162 testProperty "c2100 identity2" prop_c2100_identity2 ]
1165 p79_28_properties :: Test.Framework.Test
1167 testGroup "p. 79, Section (2.8) Properties" [
1168 testProperty "c3000 identity" prop_c3000_identity,
1169 testProperty "c2010 identity" prop_c2010_identity,
1170 testProperty "c2001 identity" prop_c2001_identity,
1171 testProperty "c1020 identity" prop_c1020_identity,
1172 testProperty "c1002 identity" prop_c1002_identity,
1173 testProperty "c1011 identity" prop_c1011_identity ]
1176 edge_incidence_tests :: Test.Framework.Test
1177 edge_incidence_tests =
1178 testGroup "Edge Incidence Tests" [
1179 testProperty "t0 shares edge with t6" prop_t0_shares_edge_with_t6,
1180 testProperty "t0 shares edge with t1" prop_t0_shares_edge_with_t1,
1181 testProperty "t0 shares edge with t3" prop_t0_shares_edge_with_t3,
1182 testProperty "t1 shares edge with t2" prop_t1_shares_edge_with_t2,
1183 testProperty "t1 shares edge with t19" prop_t1_shares_edge_with_t19,
1184 testProperty "t2 shares edge with t3" prop_t2_shares_edge_with_t3,
1185 testProperty "t2 shares edge with t12" prop_t2_shares_edge_with_t12,
1186 testProperty "t3 shares edge with t21" prop_t3_shares_edge_with_t21,
1187 testProperty "t4 shares edge with t5" prop_t4_shares_edge_with_t5,
1188 testProperty "t4 shares edge with t7" prop_t4_shares_edge_with_t7,
1189 testProperty "t4 shares edge with t10" prop_t4_shares_edge_with_t10,
1190 testProperty "t5 shares edge with t6" prop_t5_shares_edge_with_t6,
1191 testProperty "t5 shares edge with t16" prop_t5_shares_edge_with_t16,
1192 testProperty "t6 shares edge with t7" prop_t6_shares_edge_with_t7,
1193 testProperty "t7 shares edge with t20" prop_t7_shares_edge_with_t20 ]
1195 cube_properties :: Test.Framework.Test
1197 testGroup "Cube Properties" [
1201 edge_incidence_tests,
1202 testProperty "opposite octant tetrahedra are disjoint (1)"
1203 prop_opposite_octant_tetrahedra_disjoint1,
1204 testProperty "opposite octant tetrahedra are disjoint (2)"
1205 prop_opposite_octant_tetrahedra_disjoint2,
1206 testProperty "opposite octant tetrahedra are disjoint (3)"
1207 prop_opposite_octant_tetrahedra_disjoint3,
1208 testProperty "opposite octant tetrahedra are disjoint (4)"
1209 prop_opposite_octant_tetrahedra_disjoint4,
1210 testProperty "opposite octant tetrahedra are disjoint (5)"
1211 prop_opposite_octant_tetrahedra_disjoint5,
1212 testProperty "opposite octant tetrahedra are disjoint (6)"
1213 prop_opposite_octant_tetrahedra_disjoint6,
1214 testProperty "all volumes positive" prop_all_volumes_positive,
1215 testProperty "all volumes exact" prop_all_volumes_exact,
1216 testProperty "v0 all equal" prop_v0_all_equal,
1217 testProperty "interior values all identical"
1218 prop_interior_values_all_identical,
1219 testProperty "c-tilde_2100 rotation correct"
1220 prop_c_tilde_2100_rotation_correct,
1221 testProperty "c-tilde_2100 correct"
1222 prop_c_tilde_2100_correct ]