4 find_containing_tetrahedron,
9 import Data.Maybe ( fromJust )
10 import qualified Data.Vector as V (
18 import Prelude hiding ( LT )
19 import Test.Tasty ( TestTree, testGroup )
20 import Test.Tasty.QuickCheck (
21 Arbitrary( arbitrary ),
34 import Comparisons ( (~=), (~~=) )
35 import qualified Face ( Face(..), center )
36 import FunctionValues ( FunctionValues, eval, rotate )
37 import Misc ( all_equal, disjoint )
38 import Point ( Point( Point ), dot )
40 Tetrahedron(Tetrahedron, function_values, v0, v1, v2, v3),
45 data Cube = Cube { i :: !Int,
48 fv :: !FunctionValues,
49 tetrahedra_volume :: !Double }
53 instance Arbitrary Cube where
55 i' <- choose (coordmin, coordmax)
56 j' <- choose (coordmin, coordmax)
57 k' <- choose (coordmin, coordmax)
58 fv' <- arbitrary :: Gen FunctionValues
59 (Positive tet_vol) <- arbitrary :: Gen (Positive Double)
60 return (Cube i' j' k' fv' tet_vol)
62 -- The idea here is that, when cubed in the volume formula,
63 -- these numbers don't overflow 64 bits. This number is not
64 -- magic in any other sense than that it does not cause test
65 -- failures, while 2^23 does.
66 coordmax = 4194304 :: Int -- 2^22
70 instance Show Cube where
72 "Cube_" ++ subscript ++ "\n" ++
73 " Center: " ++ (show (center cube)) ++ "\n" ++
74 " xmin: " ++ (show (xmin cube)) ++ "\n" ++
75 " xmax: " ++ (show (xmax cube)) ++ "\n" ++
76 " ymin: " ++ (show (ymin cube)) ++ "\n" ++
77 " ymax: " ++ (show (ymax cube)) ++ "\n" ++
78 " zmin: " ++ (show (zmin cube)) ++ "\n" ++
79 " zmax: " ++ (show (zmax cube)) ++ "\n"
82 (show (i cube)) ++ "," ++ (show (j cube)) ++ "," ++ (show (k cube))
85 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
87 xmin :: Cube -> Double
88 xmin cube = (i' - 1/2)
90 i' = fromIntegral (i cube) :: Double
92 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
94 xmax :: Cube -> Double
95 xmax cube = (i' + 1/2)
97 i' = fromIntegral (i cube) :: Double
99 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
101 ymin :: Cube -> Double
102 ymin cube = (j' - 1/2)
104 j' = fromIntegral (j cube) :: Double
106 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
108 ymax :: Cube -> Double
109 ymax cube = (j' + 1/2)
111 j' = fromIntegral (j cube) :: Double
113 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
115 zmin :: Cube -> Double
116 zmin cube = (k' - 1/2)
118 k' = fromIntegral (k cube) :: Double
120 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
122 zmax :: Cube -> Double
123 zmax cube = (k' + 1/2)
125 k' = fromIntegral (k cube) :: Double
128 -- | The center of Cube_ijk coincides with v_ijk at
129 -- (i, j, k). See Sorokina and Zeilfelder, p. 76.
130 center :: Cube -> Point
134 x = fromIntegral (i cube) :: Double
135 y = fromIntegral (j cube) :: Double
136 z = fromIntegral (k cube) :: Double
141 -- | The top (in the direction of z) face of the cube.
142 top_face :: Cube -> Face.Face
143 top_face cube = Face.Face v0' v1' v2' v3'
145 delta = (1/2) :: Double
147 v0' = cc + ( Point delta (-delta) delta )
148 v1' = cc + ( Point delta delta delta )
149 v2' = cc + ( Point (-delta) delta delta )
150 v3' = cc + ( Point (-delta) (-delta) delta )
154 -- | The back (in the direction of x) face of the cube.
155 back_face :: Cube -> Face.Face
156 back_face cube = Face.Face v0' v1' v2' v3'
158 delta = (1/2) :: Double
160 v0' = cc + ( Point delta (-delta) (-delta) )
161 v1' = cc + ( Point delta delta (-delta) )
162 v2' = cc + ( Point delta delta delta )
163 v3' = cc + ( Point delta (-delta) delta )
166 -- The bottom face (in the direction of -z) of the cube.
167 down_face :: Cube -> Face.Face
168 down_face cube = Face.Face v0' v1' v2' v3'
170 delta = (1/2) :: Double
172 v0' = cc + ( Point (-delta) (-delta) (-delta) )
173 v1' = cc + ( Point (-delta) delta (-delta) )
174 v2' = cc + ( Point delta delta (-delta) )
175 v3' = cc + ( Point delta (-delta) (-delta) )
179 -- | The front (in the direction of -x) face of the cube.
180 front_face :: Cube -> Face.Face
181 front_face cube = Face.Face v0' v1' v2' v3'
183 delta = (1/2) :: Double
185 v0' = cc + ( Point (-delta) (-delta) delta )
186 v1' = cc + ( Point (-delta) delta delta )
187 v2' = cc + ( Point (-delta) delta (-delta) )
188 v3' = cc + ( Point (-delta) (-delta) (-delta) )
190 -- | The left (in the direction of -y) face of the cube.
191 left_face :: Cube -> Face.Face
192 left_face cube = Face.Face v0' v1' v2' v3'
194 delta = (1/2) :: Double
196 v0' = cc + ( Point delta (-delta) delta )
197 v1' = cc + ( Point (-delta) (-delta) delta )
198 v2' = cc + ( Point (-delta) (-delta) (-delta) )
199 v3' = cc + ( Point delta (-delta) (-delta) )
202 -- | The right (in the direction of y) face of the cube.
203 right_face :: Cube -> Face.Face
204 right_face cube = Face.Face v0' v1' v2' v3'
206 delta = (1/2) :: Double
208 v0' = cc + ( Point (-delta) delta delta)
209 v1' = cc + ( Point delta delta delta )
210 v2' = cc + ( Point delta delta (-delta) )
211 v3' = cc + ( Point (-delta) delta (-delta) )
214 tetrahedron :: Cube -> Int -> Tetrahedron
217 Tetrahedron (fv cube) v0' v1' v2' v3' vol
224 vol = tetrahedra_volume cube
227 Tetrahedron fv' v0' v1' v2' v3' vol
234 fv' = rotate ccwx (fv cube)
235 vol = tetrahedra_volume cube
238 Tetrahedron fv' v0' v1' v2' v3' vol
245 fv' = rotate ccwx $ rotate ccwx $ fv cube
246 vol = tetrahedra_volume cube
249 Tetrahedron fv' v0' v1' v2' v3' vol
256 fv' = rotate cwx (fv cube)
257 vol = tetrahedra_volume cube
260 Tetrahedron fv' v0' v1' v2' v3' vol
267 fv' = rotate cwy (fv cube)
268 vol = tetrahedra_volume cube
271 Tetrahedron fv' v0' v1' v2' v3' vol
278 fv' = rotate cwy $ rotate cwz $ fv cube
279 vol = tetrahedra_volume cube
282 Tetrahedron fv' v0' v1' v2' v3' vol
289 fv' = rotate cwy $ rotate cwz
292 vol = tetrahedra_volume cube
295 Tetrahedron fv' v0' v1' v2' v3' vol
302 fv' = rotate cwy $ rotate ccwz $ fv cube
303 vol = tetrahedra_volume cube
306 Tetrahedron fv' v0' v1' v2' v3' vol
313 fv' = rotate cwy $ rotate cwy $ fv cube
314 vol = tetrahedra_volume cube
317 Tetrahedron fv' v0' v1' v2' v3' vol
324 fv' = rotate cwy $ rotate cwy
327 vol = tetrahedra_volume cube
329 tetrahedron cube 10 =
330 Tetrahedron fv' v0' v1' v2' v3' vol
337 fv' = rotate cwy $ rotate cwy
342 vol = tetrahedra_volume cube
344 tetrahedron cube 11 =
345 Tetrahedron fv' v0' v1' v2' v3' vol
352 fv' = rotate cwy $ rotate cwy
355 vol = tetrahedra_volume cube
357 tetrahedron cube 12 =
358 Tetrahedron fv' v0' v1' v2' v3' vol
365 fv' = rotate ccwy $ fv cube
366 vol = tetrahedra_volume cube
368 tetrahedron cube 13 =
369 Tetrahedron fv' v0' v1' v2' v3' vol
376 fv' = rotate ccwy $ rotate ccwz $ fv cube
377 vol = tetrahedra_volume cube
379 tetrahedron cube 14 =
380 Tetrahedron fv' v0' v1' v2' v3' vol
387 fv' = rotate ccwy $ rotate ccwz
390 vol = tetrahedra_volume cube
392 tetrahedron cube 15 =
393 Tetrahedron fv' v0' v1' v2' v3' vol
400 fv' = rotate ccwy $ rotate cwz $ fv cube
401 vol = tetrahedra_volume cube
403 tetrahedron cube 16 =
404 Tetrahedron fv' v0' v1' v2' v3' vol
411 fv' = rotate ccwz $ fv cube
412 vol = tetrahedra_volume cube
414 tetrahedron cube 17 =
415 Tetrahedron fv' v0' v1' v2' v3' vol
422 fv' = rotate ccwz $ rotate cwy $ fv cube
423 vol = tetrahedra_volume cube
425 tetrahedron cube 18 =
426 Tetrahedron fv' v0' v1' v2' v3' vol
433 fv' = rotate ccwz $ rotate cwy
436 vol = tetrahedra_volume cube
438 tetrahedron cube 19 =
439 Tetrahedron fv' v0' v1' v2' v3' vol
446 fv' = rotate ccwz $ rotate ccwy
448 vol = tetrahedra_volume cube
450 tetrahedron cube 20 =
451 Tetrahedron fv' v0' v1' v2' v3' vol
458 fv' = rotate cwz $ fv cube
459 vol = tetrahedra_volume cube
461 tetrahedron cube 21 =
462 Tetrahedron fv' v0' v1' v2' v3' vol
469 fv' = rotate cwz $ rotate ccwy $ fv cube
470 vol = tetrahedra_volume cube
472 tetrahedron cube 22 =
473 Tetrahedron fv' v0' v1' v2' v3' vol
480 fv' = rotate cwz $ rotate ccwy
483 vol = tetrahedra_volume cube
485 tetrahedron cube 23 =
486 Tetrahedron fv' v0' v1' v2' v3' vol
493 fv' = rotate cwz $ rotate cwy
495 vol = tetrahedra_volume cube
498 -- Only used in tests, so we don't need the added speed
500 tetrahedra :: Cube -> [Tetrahedron]
501 tetrahedra cube = [ tetrahedron cube n | n <- [0..23] ]
503 front_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
504 front_left_top_tetrahedra cube =
505 V.singleton (tetrahedron cube 0) `V.snoc`
506 (tetrahedron cube 3) `V.snoc`
507 (tetrahedron cube 6) `V.snoc`
508 (tetrahedron cube 7) `V.snoc`
509 (tetrahedron cube 20) `V.snoc`
510 (tetrahedron cube 21)
512 front_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
513 front_left_down_tetrahedra cube =
514 V.singleton (tetrahedron cube 0) `V.snoc`
515 (tetrahedron cube 2) `V.snoc`
516 (tetrahedron cube 3) `V.snoc`
517 (tetrahedron cube 12) `V.snoc`
518 (tetrahedron cube 15) `V.snoc`
519 (tetrahedron cube 21)
521 front_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
522 front_right_top_tetrahedra cube =
523 V.singleton (tetrahedron cube 0) `V.snoc`
524 (tetrahedron cube 1) `V.snoc`
525 (tetrahedron cube 5) `V.snoc`
526 (tetrahedron cube 6) `V.snoc`
527 (tetrahedron cube 16) `V.snoc`
528 (tetrahedron cube 19)
530 front_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
531 front_right_down_tetrahedra cube =
532 V.singleton (tetrahedron cube 1) `V.snoc`
533 (tetrahedron cube 2) `V.snoc`
534 (tetrahedron cube 12) `V.snoc`
535 (tetrahedron cube 13) `V.snoc`
536 (tetrahedron cube 18) `V.snoc`
537 (tetrahedron cube 19)
539 back_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
540 back_left_top_tetrahedra cube =
541 V.singleton (tetrahedron cube 0) `V.snoc`
542 (tetrahedron cube 3) `V.snoc`
543 (tetrahedron cube 6) `V.snoc`
544 (tetrahedron cube 7) `V.snoc`
545 (tetrahedron cube 20) `V.snoc`
546 (tetrahedron cube 21)
548 back_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
549 back_left_down_tetrahedra cube =
550 V.singleton (tetrahedron cube 8) `V.snoc`
551 (tetrahedron cube 11) `V.snoc`
552 (tetrahedron cube 14) `V.snoc`
553 (tetrahedron cube 15) `V.snoc`
554 (tetrahedron cube 22) `V.snoc`
555 (tetrahedron cube 23)
557 back_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
558 back_right_top_tetrahedra cube =
559 V.singleton (tetrahedron cube 4) `V.snoc`
560 (tetrahedron cube 5) `V.snoc`
561 (tetrahedron cube 9) `V.snoc`
562 (tetrahedron cube 10) `V.snoc`
563 (tetrahedron cube 16) `V.snoc`
564 (tetrahedron cube 17)
566 back_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
567 back_right_down_tetrahedra cube =
568 V.singleton (tetrahedron cube 8) `V.snoc`
569 (tetrahedron cube 9) `V.snoc`
570 (tetrahedron cube 13) `V.snoc`
571 (tetrahedron cube 14) `V.snoc`
572 (tetrahedron cube 17) `V.snoc`
573 (tetrahedron cube 18)
575 in_top_half :: Cube -> Point -> Bool
576 in_top_half cube (Point _ _ z) =
577 distance_from_top <= distance_from_bottom
579 distance_from_top = abs $ (zmax cube) - z
580 distance_from_bottom = abs $ (zmin cube) - z
582 in_front_half :: Cube -> Point -> Bool
583 in_front_half cube (Point x _ _) =
584 distance_from_front <= distance_from_back
586 distance_from_front = abs $ (xmin cube) - x
587 distance_from_back = abs $ (xmax cube) - x
590 in_left_half :: Cube -> Point -> Bool
591 in_left_half cube (Point _ y _) =
592 distance_from_left <= distance_from_right
594 distance_from_left = abs $ (ymin cube) - y
595 distance_from_right = abs $ (ymax cube) - y
598 -- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
599 -- contain the given 'Point'. This should be faster than checking
600 -- every tetrahedron individually, since we determine which half
601 -- (hemisphere?) of the cube the point lies in three times: once in
602 -- each dimension. This allows us to eliminate non-candidates
605 -- This can throw an exception, but the use of 'head' might
606 -- save us some unnecessary computations.
608 {-# INLINE find_containing_tetrahedron #-}
609 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
610 find_containing_tetrahedron cube p =
611 candidates `V.unsafeIndex` (fromJust lucky_idx)
613 front_half = in_front_half cube p
614 top_half = in_top_half cube p
615 left_half = in_left_half cube p
617 candidates :: V.Vector Tetrahedron
622 front_left_top_tetrahedra cube
624 front_left_down_tetrahedra cube
627 front_right_top_tetrahedra cube
629 front_right_down_tetrahedra cube
631 | otherwise = -- back half
634 back_left_top_tetrahedra cube
636 back_left_down_tetrahedra cube
639 back_right_top_tetrahedra cube
641 back_right_down_tetrahedra cube
643 -- Use the dot product instead of Euclidean distance here to save
644 -- a sqrt(). So, "distances" below really means "distances
646 distances :: V.Vector Double
647 distances = V.map ((dot p) . barycenter) candidates
649 shortest_distance :: Double
650 shortest_distance = V.minimum distances
652 -- Compute the index of the tetrahedron with the center closest to
653 -- p. This is a bad algorithm, but don't change it! If you make it
654 -- smarter by finding the index of shortest_distance in distances
655 -- (this should give the same answer and avoids recomputing the
656 -- dot product), the program gets slower. Seriously!
657 lucky_idx :: Maybe Int
658 lucky_idx = V.findIndex
659 (\t -> (barycenter t) `dot` p == shortest_distance)
669 prop_opposite_octant_tetrahedra_disjoint1 :: Cube -> Bool
670 prop_opposite_octant_tetrahedra_disjoint1 cube =
671 disjoint (front_left_top_tetrahedra cube) (front_right_down_tetrahedra cube)
673 prop_opposite_octant_tetrahedra_disjoint2 :: Cube -> Bool
674 prop_opposite_octant_tetrahedra_disjoint2 cube =
675 disjoint (back_left_top_tetrahedra cube) (back_right_down_tetrahedra cube)
677 prop_opposite_octant_tetrahedra_disjoint3 :: Cube -> Bool
678 prop_opposite_octant_tetrahedra_disjoint3 cube =
679 disjoint (front_left_top_tetrahedra cube) (back_right_top_tetrahedra cube)
681 prop_opposite_octant_tetrahedra_disjoint4 :: Cube -> Bool
682 prop_opposite_octant_tetrahedra_disjoint4 cube =
683 disjoint (front_left_down_tetrahedra cube) (back_right_down_tetrahedra cube)
685 prop_opposite_octant_tetrahedra_disjoint5 :: Cube -> Bool
686 prop_opposite_octant_tetrahedra_disjoint5 cube =
687 disjoint (front_left_top_tetrahedra cube) (back_left_down_tetrahedra cube)
689 prop_opposite_octant_tetrahedra_disjoint6 :: Cube -> Bool
690 prop_opposite_octant_tetrahedra_disjoint6 cube =
691 disjoint (front_right_top_tetrahedra cube) (back_right_down_tetrahedra cube)
694 -- | Since the grid size is necessarily positive, all tetrahedra
695 -- (which comprise cubes of positive volume) must have positive
697 prop_all_volumes_positive :: Cube -> Bool
698 prop_all_volumes_positive cube =
702 volumes = map volume ts
705 -- | In fact, since all of the tetrahedra are identical, we should
706 -- already know their volumes. There's 24 tetrahedra to a cube, so
707 -- we'd expect the volume of each one to be 1/24.
708 prop_all_volumes_exact :: Cube -> Bool
709 prop_all_volumes_exact cube =
710 and [volume t ~~= 1/24 | t <- tetrahedra cube]
712 -- | All tetrahedron should have their v0 located at the center of the
714 prop_v0_all_equal :: Cube -> Bool
715 prop_v0_all_equal cube = (v0 t0) == (v0 t1)
717 t0 = head (tetrahedra cube) -- Doesn't matter which two we choose.
718 t1 = head $ tail (tetrahedra cube)
721 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Note that the
722 -- third and fourth indices of c-t3 have been switched. This is
723 -- because we store the triangles oriented such that their volume is
724 -- positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde point
725 -- in opposite directions, one of them has to have negative volume!
726 prop_c0120_identity1 :: Cube -> Bool
727 prop_c0120_identity1 cube =
728 c t0 0 1 2 0 ~= (c t0 0 0 2 1 + c t3 0 0 1 2) / 2
730 t0 = tetrahedron cube 0
731 t3 = tetrahedron cube 3
734 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
735 -- 'prop_c0120_identity1' with tetrahedrons 1 and 2.
736 prop_c0120_identity2 :: Cube -> Bool
737 prop_c0120_identity2 cube =
738 c t1 0 1 2 0 ~= (c t1 0 0 2 1 + c t0 0 0 1 2) / 2
740 t0 = tetrahedron cube 0
741 t1 = tetrahedron cube 1
743 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
744 -- 'prop_c0120_identity1' with tetrahedrons 1 and 2.
745 prop_c0120_identity3 :: Cube -> Bool
746 prop_c0120_identity3 cube =
747 c t2 0 1 2 0 ~= (c t2 0 0 2 1 + c t1 0 0 1 2) / 2
749 t1 = tetrahedron cube 1
750 t2 = tetrahedron cube 2
752 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
753 -- 'prop_c0120_identity1' with tetrahedrons 2 and 3.
754 prop_c0120_identity4 :: Cube -> Bool
755 prop_c0120_identity4 cube =
756 c t3 0 1 2 0 ~= (c t3 0 0 2 1 + c t2 0 0 1 2) / 2
758 t2 = tetrahedron cube 2
759 t3 = tetrahedron cube 3
762 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
763 -- 'prop_c0120_identity1' with tetrahedrons 4 and 5.
764 prop_c0120_identity5 :: Cube -> Bool
765 prop_c0120_identity5 cube =
766 c t5 0 1 2 0 ~= (c t5 0 0 2 1 + c t4 0 0 1 2) / 2
768 t4 = tetrahedron cube 4
769 t5 = tetrahedron cube 5
771 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
772 -- 'prop_c0120_identity1' with tetrahedrons 5 and 6.
773 prop_c0120_identity6 :: Cube -> Bool
774 prop_c0120_identity6 cube =
775 c t6 0 1 2 0 ~= (c t6 0 0 2 1 + c t5 0 0 1 2) / 2
777 t5 = tetrahedron cube 5
778 t6 = tetrahedron cube 6
781 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). Repeats
782 -- 'prop_c0120_identity1' with tetrahedrons 6 and 7.
783 prop_c0120_identity7 :: Cube -> Bool
784 prop_c0120_identity7 cube =
785 c t7 0 1 2 0 ~= (c t7 0 0 2 1 + c t6 0 0 1 2) / 2
787 t6 = tetrahedron cube 6
788 t7 = tetrahedron cube 7
791 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
792 -- 'prop_c0120_identity1'.
793 prop_c0210_identity1 :: Cube -> Bool
794 prop_c0210_identity1 cube =
795 c t0 0 2 1 0 ~= (c t0 0 1 1 1 + c t3 0 1 1 1) / 2
797 t0 = tetrahedron cube 0
798 t3 = tetrahedron cube 3
801 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
802 -- 'prop_c0120_identity1'.
803 prop_c0300_identity1 :: Cube -> Bool
804 prop_c0300_identity1 cube =
805 c t0 0 3 0 0 ~= (c t0 0 2 0 1 + c t3 0 2 1 0) / 2
807 t0 = tetrahedron cube 0
808 t3 = tetrahedron cube 3
811 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
812 -- 'prop_c0120_identity1'.
813 prop_c1110_identity :: Cube -> Bool
814 prop_c1110_identity cube =
815 c t0 1 1 1 0 ~= (c t0 1 0 1 1 + c t3 1 0 1 1) / 2
817 t0 = tetrahedron cube 0
818 t3 = tetrahedron cube 3
821 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
822 -- 'prop_c0120_identity1'.
823 prop_c1200_identity1 :: Cube -> Bool
824 prop_c1200_identity1 cube =
825 c t0 1 2 0 0 ~= (c t0 1 1 0 1 + c t3 1 1 1 0) / 2
827 t0 = tetrahedron cube 0
828 t3 = tetrahedron cube 3
831 -- | Given in Sorokina and Zeilfelder, p. 79, (2.6). See
832 -- 'prop_c0120_identity1'.
833 prop_c2100_identity1 :: Cube -> Bool
834 prop_c2100_identity1 cube =
835 c t0 2 1 0 0 ~= (c t0 2 0 0 1 + c t3 2 0 1 0) / 2
837 t0 = tetrahedron cube 0
838 t3 = tetrahedron cube 3
842 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). Note that the
843 -- third and fourth indices of c-t3 have been switched. This is
844 -- because we store the triangles oriented such that their volume is
845 -- positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde
846 -- point in opposite directions, one of them has to have negative
848 prop_c0102_identity1 :: Cube -> Bool
849 prop_c0102_identity1 cube =
850 c t0 0 1 0 2 ~= (c t0 0 0 1 2 + c t1 0 0 2 1) / 2
852 t0 = tetrahedron cube 0
853 t1 = tetrahedron cube 1
856 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
857 -- 'prop_c0102_identity1'.
858 prop_c0201_identity1 :: Cube -> Bool
859 prop_c0201_identity1 cube =
860 c t0 0 2 0 1 ~= (c t0 0 1 1 1 + c t1 0 1 1 1) / 2
862 t0 = tetrahedron cube 0
863 t1 = tetrahedron cube 1
866 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
867 -- 'prop_c0102_identity1'.
868 prop_c0300_identity2 :: Cube -> Bool
869 prop_c0300_identity2 cube =
870 c t0 0 3 0 0 ~= (c t0 0 2 1 0 + c t1 0 2 0 1) / 2
872 t0 = tetrahedron cube 0
873 t1 = tetrahedron cube 1
876 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
877 -- 'prop_c0102_identity1'.
878 prop_c1101_identity :: Cube -> Bool
879 prop_c1101_identity cube =
880 c t0 1 1 0 1 ~= (c t0 1 0 1 1 + c t1 1 0 1 1) / 2
882 t0 = tetrahedron cube 0
883 t1 = tetrahedron cube 1
886 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
887 -- 'prop_c0102_identity1'.
888 prop_c1200_identity2 :: Cube -> Bool
889 prop_c1200_identity2 cube =
890 c t0 1 2 0 0 ~= (c t0 1 1 1 0 + c t1 1 1 0 1) / 2
892 t0 = tetrahedron cube 0
893 t1 = tetrahedron cube 1
896 -- | Given in Sorokina and Zeilfelder, p. 79, (2.7). See
897 -- 'prop_c0102_identity1'.
898 prop_c2100_identity2 :: Cube -> Bool
899 prop_c2100_identity2 cube =
900 c t0 2 1 0 0 ~= (c t0 2 0 1 0 + c t1 2 0 0 1) / 2
902 t0 = tetrahedron cube 0
903 t1 = tetrahedron cube 1
906 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). The third and
907 -- fourth indices of c-t6 have been switched. This is because we
908 -- store the triangles oriented such that their volume is
909 -- positive. If T and T-tilde share \<v0,v1,v2\> and v3,v3-tilde
910 -- point in opposite directions, one of them has to have negative
912 prop_c3000_identity :: Cube -> Bool
913 prop_c3000_identity cube =
914 c t0 3 0 0 0 ~= c t0 2 1 0 0 + c t6 2 1 0 0
915 - ((c t0 2 0 1 0 + c t0 2 0 0 1)/ 2)
917 t0 = tetrahedron cube 0
918 t6 = tetrahedron cube 6
921 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
922 -- 'prop_c3000_identity'.
923 prop_c2010_identity :: Cube -> Bool
924 prop_c2010_identity cube =
925 c t0 2 0 1 0 ~= c t0 1 1 1 0 + c t6 1 1 0 1
926 - ((c t0 1 0 2 0 + c t0 1 0 1 1)/ 2)
928 t0 = tetrahedron cube 0
929 t6 = tetrahedron cube 6
932 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
933 -- 'prop_c3000_identity'.
934 prop_c2001_identity :: Cube -> Bool
935 prop_c2001_identity cube =
936 c t0 2 0 0 1 ~= c t0 1 1 0 1 + c t6 1 1 1 0
937 - ((c t0 1 0 0 2 + c t0 1 0 1 1)/ 2)
939 t0 = tetrahedron cube 0
940 t6 = tetrahedron cube 6
943 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
944 -- 'prop_c3000_identity'.
945 prop_c1020_identity :: Cube -> Bool
946 prop_c1020_identity cube =
947 c t0 1 0 2 0 ~= c t0 0 1 2 0 + c t6 0 1 0 2
948 - ((c t0 0 0 3 0 + c t0 0 0 2 1)/ 2)
950 t0 = tetrahedron cube 0
951 t6 = tetrahedron cube 6
954 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
955 -- 'prop_c3000_identity'.
956 prop_c1002_identity :: Cube -> Bool
957 prop_c1002_identity cube =
958 c t0 1 0 0 2 ~= c t0 0 1 0 2 + c t6 0 1 2 0
959 - ((c t0 0 0 0 3 + c t0 0 0 1 2)/ 2)
961 t0 = tetrahedron cube 0
962 t6 = tetrahedron cube 6
965 -- | Given in Sorokina and Zeilfelder, p. 79, (2.8). See
966 -- 'prop_c3000_identity'.
967 prop_c1011_identity :: Cube -> Bool
968 prop_c1011_identity cube =
969 c t0 1 0 1 1 ~= c t0 0 1 1 1 + c t6 0 1 1 1 -
970 ((c t0 0 0 1 2 + c t0 0 0 2 1)/ 2)
972 t0 = tetrahedron cube 0
973 t6 = tetrahedron cube 6
976 -- | The function values at the interior should be the same for all
978 prop_interior_values_all_identical :: Cube -> Bool
979 prop_interior_values_all_identical cube =
980 all_equal [ eval (function_values tet) I | tet <- tetrahedra cube ]
983 -- | We know what (c t6 2 1 0 0) should be from Sorokina and Zeilfelder, p. 87.
984 -- This test checks the rotation works as expected.
985 prop_c_tilde_2100_rotation_correct :: Cube -> Bool
986 prop_c_tilde_2100_rotation_correct cube =
989 t0 = tetrahedron cube 0
990 t6 = tetrahedron cube 6
992 -- What gets computed for c2100 of t6.
993 expr1 = eval (function_values t6) $
995 (1/12)*(T + R + L + D) +
996 (1/64)*(FT + FR + FL + FD) +
999 (1/96)*(RT + LD + LT + RD) +
1000 (1/192)*(BT + BR + BL + BD)
1002 -- What should be computed for c2100 of t6.
1003 expr2 = eval (function_values t0) $
1005 (1/12)*(F + R + L + B) +
1006 (1/64)*(FT + RT + LT + BT) +
1009 (1/96)*(FR + FL + BR + BL) +
1010 (1/192)*(FD + RD + LD + BD)
1013 -- | We know what (c t6 2 1 0 0) should be from Sorokina and
1014 -- Zeilfelder, p. 87. This test checks the actual value based on
1015 -- the FunctionValues of the cube.
1017 -- If 'prop_c_tilde_2100_rotation_correct' passes, then this test is
1019 prop_c_tilde_2100_correct :: Cube -> Bool
1020 prop_c_tilde_2100_correct cube =
1021 c t6 2 1 0 0 ~= expected
1023 t0 = tetrahedron cube 0
1024 t6 = tetrahedron cube 6
1025 fvs = function_values t0
1026 expected = eval fvs $
1028 (1/12)*(F + R + L + B) +
1029 (1/64)*(FT + RT + LT + BT) +
1032 (1/96)*(FR + FL + BR + BL) +
1033 (1/192)*(FD + RD + LD + BD)
1036 -- Tests to check that the correct edges are incidental.
1037 prop_t0_shares_edge_with_t1 :: Cube -> Bool
1038 prop_t0_shares_edge_with_t1 cube =
1039 (v1 t0) == (v1 t1) && (v3 t0) == (v2 t1)
1041 t0 = tetrahedron cube 0
1042 t1 = tetrahedron cube 1
1044 prop_t0_shares_edge_with_t3 :: Cube -> Bool
1045 prop_t0_shares_edge_with_t3 cube =
1046 (v1 t0) == (v1 t3) && (v2 t0) == (v3 t3)
1048 t0 = tetrahedron cube 0
1049 t3 = tetrahedron cube 3
1051 prop_t0_shares_edge_with_t6 :: Cube -> Bool
1052 prop_t0_shares_edge_with_t6 cube =
1053 (v2 t0) == (v3 t6) && (v3 t0) == (v2 t6)
1055 t0 = tetrahedron cube 0
1056 t6 = tetrahedron cube 6
1058 prop_t1_shares_edge_with_t2 :: Cube -> Bool
1059 prop_t1_shares_edge_with_t2 cube =
1060 (v1 t1) == (v1 t2) && (v3 t1) == (v2 t2)
1062 t1 = tetrahedron cube 1
1063 t2 = tetrahedron cube 2
1065 prop_t1_shares_edge_with_t19 :: Cube -> Bool
1066 prop_t1_shares_edge_with_t19 cube =
1067 (v2 t1) == (v3 t19) && (v3 t1) == (v2 t19)
1069 t1 = tetrahedron cube 1
1070 t19 = tetrahedron cube 19
1072 prop_t2_shares_edge_with_t3 :: Cube -> Bool
1073 prop_t2_shares_edge_with_t3 cube =
1074 (v1 t1) == (v1 t2) && (v3 t1) == (v2 t2)
1076 t1 = tetrahedron cube 1
1077 t2 = tetrahedron cube 2
1079 prop_t2_shares_edge_with_t12 :: Cube -> Bool
1080 prop_t2_shares_edge_with_t12 cube =
1081 (v2 t2) == (v3 t12) && (v3 t2) == (v2 t12)
1083 t2 = tetrahedron cube 2
1084 t12 = tetrahedron cube 12
1086 prop_t3_shares_edge_with_t21 :: Cube -> Bool
1087 prop_t3_shares_edge_with_t21 cube =
1088 (v2 t3) == (v3 t21) && (v3 t3) == (v2 t21)
1090 t3 = tetrahedron cube 3
1091 t21 = tetrahedron cube 21
1093 prop_t4_shares_edge_with_t5 :: Cube -> Bool
1094 prop_t4_shares_edge_with_t5 cube =
1095 (v1 t4) == (v1 t5) && (v3 t4) == (v2 t5)
1097 t4 = tetrahedron cube 4
1098 t5 = tetrahedron cube 5
1100 prop_t4_shares_edge_with_t7 :: Cube -> Bool
1101 prop_t4_shares_edge_with_t7 cube =
1102 (v1 t4) == (v1 t7) && (v2 t4) == (v3 t7)
1104 t4 = tetrahedron cube 4
1105 t7 = tetrahedron cube 7
1107 prop_t4_shares_edge_with_t10 :: Cube -> Bool
1108 prop_t4_shares_edge_with_t10 cube =
1109 (v2 t4) == (v3 t10) && (v3 t4) == (v2 t10)
1111 t4 = tetrahedron cube 4
1112 t10 = tetrahedron cube 10
1114 prop_t5_shares_edge_with_t6 :: Cube -> Bool
1115 prop_t5_shares_edge_with_t6 cube =
1116 (v1 t5) == (v1 t6) && (v3 t5) == (v2 t6)
1118 t5 = tetrahedron cube 5
1119 t6 = tetrahedron cube 6
1121 prop_t5_shares_edge_with_t16 :: Cube -> Bool
1122 prop_t5_shares_edge_with_t16 cube =
1123 (v2 t5) == (v3 t16) && (v3 t5) == (v2 t16)
1125 t5 = tetrahedron cube 5
1126 t16 = tetrahedron cube 16
1128 prop_t6_shares_edge_with_t7 :: Cube -> Bool
1129 prop_t6_shares_edge_with_t7 cube =
1130 (v1 t6) == (v1 t7) && (v3 t6) == (v2 t7)
1132 t6 = tetrahedron cube 6
1133 t7 = tetrahedron cube 7
1135 prop_t7_shares_edge_with_t20 :: Cube -> Bool
1136 prop_t7_shares_edge_with_t20 cube =
1137 (v2 t7) == (v3 t20) && (v2 t7) == (v3 t20)
1139 t7 = tetrahedron cube 7
1140 t20 = tetrahedron cube 20
1143 p79_26_properties :: TestTree
1145 testGroup "p. 79, Section (2.6) properties" [
1146 testProperty "c0120 identity1" prop_c0120_identity1,
1147 testProperty "c0120 identity2" prop_c0120_identity2,
1148 testProperty "c0120 identity3" prop_c0120_identity3,
1149 testProperty "c0120 identity4" prop_c0120_identity4,
1150 testProperty "c0120 identity5" prop_c0120_identity5,
1151 testProperty "c0120 identity6" prop_c0120_identity6,
1152 testProperty "c0120 identity7" prop_c0120_identity7,
1153 testProperty "c0210 identity1" prop_c0210_identity1,
1154 testProperty "c0300 identity1" prop_c0300_identity1,
1155 testProperty "c1110 identity" prop_c1110_identity,
1156 testProperty "c1200 identity1" prop_c1200_identity1,
1157 testProperty "c2100 identity1" prop_c2100_identity1]
1159 p79_27_properties :: TestTree
1161 testGroup "p. 79, Section (2.7) properties" [
1162 testProperty "c0102 identity1" prop_c0102_identity1,
1163 testProperty "c0201 identity1" prop_c0201_identity1,
1164 testProperty "c0300 identity2" prop_c0300_identity2,
1165 testProperty "c1101 identity" prop_c1101_identity,
1166 testProperty "c1200 identity2" prop_c1200_identity2,
1167 testProperty "c2100 identity2" prop_c2100_identity2 ]
1170 p79_28_properties :: TestTree
1172 testGroup "p. 79, Section (2.8) properties" [
1173 testProperty "c3000 identity" prop_c3000_identity,
1174 testProperty "c2010 identity" prop_c2010_identity,
1175 testProperty "c2001 identity" prop_c2001_identity,
1176 testProperty "c1020 identity" prop_c1020_identity,
1177 testProperty "c1002 identity" prop_c1002_identity,
1178 testProperty "c1011 identity" prop_c1011_identity ]
1181 edge_incidence_tests :: TestTree
1182 edge_incidence_tests =
1183 testGroup "Edge incidence tests" [
1184 testProperty "t0 shares edge with t6" prop_t0_shares_edge_with_t6,
1185 testProperty "t0 shares edge with t1" prop_t0_shares_edge_with_t1,
1186 testProperty "t0 shares edge with t3" prop_t0_shares_edge_with_t3,
1187 testProperty "t1 shares edge with t2" prop_t1_shares_edge_with_t2,
1188 testProperty "t1 shares edge with t19" prop_t1_shares_edge_with_t19,
1189 testProperty "t2 shares edge with t3" prop_t2_shares_edge_with_t3,
1190 testProperty "t2 shares edge with t12" prop_t2_shares_edge_with_t12,
1191 testProperty "t3 shares edge with t21" prop_t3_shares_edge_with_t21,
1192 testProperty "t4 shares edge with t5" prop_t4_shares_edge_with_t5,
1193 testProperty "t4 shares edge with t7" prop_t4_shares_edge_with_t7,
1194 testProperty "t4 shares edge with t10" prop_t4_shares_edge_with_t10,
1195 testProperty "t5 shares edge with t6" prop_t5_shares_edge_with_t6,
1196 testProperty "t5 shares edge with t16" prop_t5_shares_edge_with_t16,
1197 testProperty "t6 shares edge with t7" prop_t6_shares_edge_with_t7,
1198 testProperty "t7 shares edge with t20" prop_t7_shares_edge_with_t20 ]
1200 cube_properties :: TestTree
1202 testGroup "Cube properties" [
1206 edge_incidence_tests,
1207 testProperty "opposite octant tetrahedra are disjoint (1)"
1208 prop_opposite_octant_tetrahedra_disjoint1,
1209 testProperty "opposite octant tetrahedra are disjoint (2)"
1210 prop_opposite_octant_tetrahedra_disjoint2,
1211 testProperty "opposite octant tetrahedra are disjoint (3)"
1212 prop_opposite_octant_tetrahedra_disjoint3,
1213 testProperty "opposite octant tetrahedra are disjoint (4)"
1214 prop_opposite_octant_tetrahedra_disjoint4,
1215 testProperty "opposite octant tetrahedra are disjoint (5)"
1216 prop_opposite_octant_tetrahedra_disjoint5,
1217 testProperty "opposite octant tetrahedra are disjoint (6)"
1218 prop_opposite_octant_tetrahedra_disjoint6,
1219 testProperty "all volumes positive" prop_all_volumes_positive,
1220 testProperty "all volumes exact" prop_all_volumes_exact,
1221 testProperty "v0 all equal" prop_v0_all_equal,
1222 testProperty "interior values all identical"
1223 prop_interior_values_all_identical,
1224 testProperty "c-tilde_2100 rotation correct"
1225 prop_c_tilde_2100_rotation_correct,
1226 testProperty "c-tilde_2100 correct"
1227 prop_c_tilde_2100_correct ]