675 s6 = s8 + y[4] * y[0] * x[4] * z[3] - 2.0 * y[4] * z[0] * x[4] * y[7] +
676 2.0 * y[4] * x[0] * y[7] * z[4] - y[4] * z[0] * x[4] * y[3] -
677 y[4] * x[0] * y[7] * z[3] + y[4] * z[0] * x[3] * y[7] -
678 y[4] * y[0] * x[3] * z[4] + y[0] * x[4] * y[3] * z[7] -
679 y[0] * x[7] * y[3] * z[4] - y[0] * x[3] * y[4] * z[7] +
680 y[0] * x[7] * y[4] * z[3] + x[2] * y[7] * y[7] * z[3] -
681 z[2] * x[3] * y[7] * y[7] - 2.0 * z[2] * x[0] * y[3] * y[3] +
682 2.0 * y[0] * z[1] * x[0] * y[4] + s7;
683 s8 = -2.0 * y[0] * y[1] * x[0] * z[4] - y[0] * y[1] * x[0] * z[5] -
684 y[0] * y[0] * x[3] * z[7] - z[1] * x[0] * y[3] * y[3] -
685 y[0] * x[1] * y[5] * z[0] - 2.0 * z[0] * x[7] * y[3] * y[3] +
686 x[0] * y[3] * y[3] * z[4] + 2.0 * x[0] * y[3] * y[3] * z[7] -
687 z[0] * x[4] * y[3] * y[3] + 2.0 * x[2] * y[3] * y[3] * z[0] +
688 x[1] * y[3] * y[3] * z[0] + 2.0 * y[7] * z[6] * x[7] * y[3] +
689 2.0 * y[7] * y[6] * x[3] * z[7] - 2.0 * y[7] * y[6] * x[7] * z[3] -
690 2.0 * y[7] * x[6] * y[3] * z[7];
691 s7 = s8 + y[4] * x[4] * y[3] * z[7] - y[4] * x[4] * y[7] * z[3] +
692 y[4] * x[3] * y[7] * z[4] - y[4] * x[7] * y[3] * z[4] +
693 2.0 * y[4] * y[0] * x[4] * z[7] - 2.0 * y[4] * y[0] * x[7] * z[4] +
694 2.0 * x[6] * y[7] * y[7] * z[3] + y[4] * x[0] * y[3] * z[4] +
695 y[0] * y[1] * x[5] * z[0] + y[0] * z[1] * x[0] * y[5] -
696 x[2] * y[0] * y[0] * z[3] + x[4] * y[3] * y[3] * z[7] -
697 x[7] * y[3] * y[3] * z[4] - x[5] * y[4] * y[4] * z[1] +
698 y[3] * z[0] * x[3] * y[4];
699 s8 = y[3] * y[0] * x[4] * z[3] + 2.0 * y[3] * y[0] * x[7] * z[3] +
700 2.0 * y[3] * y[2] * x[0] * z[3] - 2.0 * y[3] * y[2] * x[3] * z[0] +
701 2.0 * y[3] * z[2] * x[3] * y[0] + y[3] * z[1] * x[3] * y[0] -
702 2.0 * y[3] * x[2] * y[0] * z[3] - y[3] * x[1] * y[0] * z[3] -
703 y[3] * y[1] * x[3] * z[0] - 2.0 * y[3] * x[0] * y[7] * z[3] -
704 y[3] * x[0] * y[4] * z[3] - 2.0 * y[3] * y[0] * x[3] * z[7] -
705 y[3] * y[0] * x[3] * z[4] + 2.0 * y[3] * z[0] * x[3] * y[7] +
706 y[3] * y[1] * x[0] * z[3] + z[5] * x[1] * y[4] * y[4];
707 s5 = s8 - 2.0 * y[0] * y[0] * x[3] * z[4] -
708 2.0 * y[0] * x[1] * y[4] * z[0] + y[3] * x[7] * y[4] * z[3] -
709 y[3] * x[4] * y[7] * z[3] + y[3] * x[3] * y[7] * z[4] -
710 y[3] * x[3] * y[4] * z[7] + y[3] * x[0] * y[7] * z[4] -
711 y[3] * z[0] * x[4] * y[7] - 2.0 * y[4] * y[5] * x[0] * z[4] + s6 +
712 y[7] * x[0] * y[3] * z[7] - y[7] * z[0] * x[7] * y[3] +
713 y[7] * y[0] * x[7] * z[3] - y[7] * y[0] * x[3] * z[7] +
714 2.0 * y[0] * y[1] * x[4] * z[0] + s7;
715 s8 = -2.0 * y[7] * x[7] * y[3] * z[4] -
716 2.0 * y[7] * x[3] * y[4] * z[7] + 2.0 * y[7] * x[4] * y[3] * z[7] +
717 y[7] * y[0] * x[4] * z[7] - y[7] * y[0] * x[7] * z[4] +
718 2.0 * y[7] * x[7] * y[4] * z[3] - y[7] * x[0] * y[4] * z[7] +
719 y[7] * z[0] * x[7] * y[4] + z[5] * x[4] * y[7] * y[7] +
720 2.0 * z[6] * x[4] * y[7] * y[7] - x[5] * y[7] * y[7] * z[4] -
721 2.0 * x[6] * y[7] * y[7] * z[4] + 2.0 * y[7] * x[6] * y[4] * z[7] -
722 2.0 * y[7] * z[6] * x[7] * y[4] + 2.0 * y[7] * y[6] * x[7] * z[4];
723 s7 = s8 - 2.0 * y[7] * y[6] * x[4] * z[7] - y[7] * z[5] * x[7] * y[4] -
724 y[7] * y[5] * x[4] * z[7] - x[0] * y[7] * y[7] * z[3] +
725 z[0] * x[3] * y[7] * y[7] + y[7] * x[5] * y[4] * z[7] +
726 y[7] * y[5] * x[7] * z[4] - y[4] * x[1] * y[5] * z[0] -
727 x[1] * y[0] * y[0] * z[2] - y[4] * y[5] * x[1] * z[4] -
728 2.0 * y[4] * z[5] * x[4] * y[0] - y[4] * y[1] * x[0] * z[4] +
729 y[4] * y[5] * x[4] * z[1] + y[0] * z[0] * x[3] * y[7] -
730 y[0] * z[1] * x[0] * y[2];
731 s8 = 2.0 * y[0] * x[1] * y[3] * z[0] + y[4] * y[1] * x[4] * z[0] +
732 2.0 * y[0] * y[1] * x[0] * z[3] + y[4] * x[1] * y[0] * z[5] -
733 y[4] * z[1] * x[5] * y[0] + y[4] * z[1] * x[0] * y[5] -
734 y[4] * z[1] * x[4] * y[0] + y[4] * x[1] * y[0] * z[4] -
735 y[4] * z[5] * x[4] * y[1] + x[5] * y[4] * y[4] * z[6] -
736 z[5] * x[6] * y[4] * y[4] + y[4] * x[5] * y[1] * z[4] -
737 y[0] * z[2] * x[0] * y[3] + y[0] * y[5] * x[4] * z[0] +
738 y[0] * x[1] * y[2] * z[0];
739 s6 = s8 - 2.0 * y[0] * z[0] * x[4] * y[3] -
740 2.0 * y[0] * x[0] * y[4] * z[3] - 2.0 * y[0] * z[1] * x[0] * y[3] -
741 y[0] * x[0] * y[7] * z[3] - 2.0 * y[0] * y[1] * x[3] * z[0] +
742 y[0] * x[2] * y[3] * z[0] - y[0] * y[1] * x[2] * z[0] +
743 y[0] * y[1] * x[0] * z[2] - y[0] * x[2] * y[1] * z[3] +
744 y[0] * x[0] * y[3] * z[7] + y[0] * x[2] * y[3] * z[1] -
745 y[0] * y[2] * x[3] * z[0] + y[0] * y[2] * x[0] * z[3] -
746 y[0] * y[5] * x[0] * z[4] - y[4] * y[5] * x[4] * z[6] + s7;
747 s8 = s6 + y[4] * z[6] * x[5] * y[7] - y[4] * x[6] * y[7] * z[5] +
748 y[4] * x[6] * y[5] * z[7] - y[4] * z[6] * x[7] * y[5] -
749 y[4] * x[5] * y[6] * z[4] + y[4] * z[5] * x[4] * y[6] +
750 y[4] * y[5] * x[6] * z[4] - 2.0 * y[1] * y[1] * x[0] * z[5] +
751 2.0 * y[1] * y[1] * x[5] * z[0] - 2.0 * y[2] * y[2] * x[6] * z[3] +
752 x[5] * y[1] * y[1] * z[4] - z[5] * x[4] * y[1] * y[1] -
753 x[6] * y[2] * y[2] * z[7] + z[6] * x[7] * y[2] * y[2];
754 s7 = s8 - x[1] * y[5] * y[5] * z[0] + z[1] * x[0] * y[5] * y[5] +
755 y[1] * y[5] * x[4] * z[1] - y[1] * y[5] * x[1] * z[4] -
756 2.0 * y[2] * z[2] * x[3] * y[6] + 2.0 * y[1] * z[1] * x[0] * y[5] -
757 2.0 * y[1] * z[1] * x[5] * y[0] + 2.0 * y[1] * x[1] * y[0] * z[5] -
758 y[2] * x[2] * y[3] * z[7] - y[2] * z[2] * x[3] * y[7] +
759 y[2] * x[2] * y[7] * z[3] + y[2] * z[2] * x[7] * y[3] -
760 2.0 * y[2] * x[2] * y[3] * z[6] + 2.0 * y[2] * x[2] * y[6] * z[3] +
761 2.0 * y[2] * z[2] * x[6] * y[3] - y[3] * y[2] * x[6] * z[3];
762 s8 = y[3] * y[2] * x[3] * z[6] + y[3] * x[2] * y[6] * z[3] -
763 y[3] * z[2] * x[3] * y[6] - y[2] * y[2] * x[7] * z[3] +
764 2.0 * y[2] * y[2] * x[3] * z[6] + y[2] * y[2] * x[3] * z[7] -
765 2.0 * y[1] * x[1] * y[5] * z[0] - x[2] * y[3] * y[3] * z[6] +
766 z[2] * x[6] * y[3] * y[3] + 2.0 * y[6] * x[2] * y[5] * z[6] +
767 2.0 * y[6] * x[6] * y[2] * z[5] - 2.0 * y[6] * x[5] * y[2] * z[6] +
768 2.0 * y[3] * x[2] * y[7] * z[3] - 2.0 * y[3] * z[2] * x[3] * y[7] -
769 y[0] * z[0] * x[7] * y[3] - y[0] * z[2] * x[1] * y[3];
770 s4 = s8 - y[2] * y[6] * x[7] * z[2] + y[0] * z[2] * x[3] * y[1] +
771 y[1] * z[5] * x[1] * y[4] - y[1] * x[5] * y[4] * z[1] +
772 2.0 * y[0] * z[0] * x[3] * y[4] + 2.0 * y[0] * x[0] * y[3] * z[4] +
773 2.0 * z[2] * x[7] * y[3] * y[3] - 2.0 * z[5] * x[7] * y[4] * y[4] +
774 x[6] * y[4] * y[4] * z[7] - z[6] * x[7] * y[4] * y[4] +
775 y[1] * y[1] * x[0] * z[3] + y[3] * x[6] * y[7] * z[2] -
776 y[3] * z[6] * x[2] * y[7] + 2.0 * y[3] * y[2] * x[3] * z[7] + s5 +
778 s8 = s4 + y[2] * x[6] * y[7] * z[2] - y[2] * y[6] * x[7] * z[3] +
779 y[2] * y[6] * x[2] * z[7] - y[2] * z[6] * x[2] * y[7] -
780 y[2] * x[6] * y[3] * z[7] + y[2] * y[6] * x[3] * z[7] +
781 y[2] * z[6] * x[7] * y[3] - 2.0 * y[3] * y[2] * x[7] * z[3] -
782 x[6] * y[3] * y[3] * z[7] + y[1] * y[1] * x[4] * z[0] -
783 y[1] * y[1] * x[3] * z[0] + x[2] * y[6] * y[6] * z[3] -
784 z[2] * x[3] * y[6] * y[6] - y[1] * y[1] * x[0] * z[4];
785 s7 = s8 + y[5] * x[1] * y[0] * z[5] + y[6] * x[2] * y[7] * z[3] -
786 y[6] * y[2] * x[6] * z[3] + y[6] * y[2] * x[3] * z[6] -
787 y[6] * x[2] * y[3] * z[6] + y[6] * z[2] * x[6] * y[3] -
788 y[5] * y[1] * x[0] * z[5] - y[5] * z[1] * x[5] * y[0] +
789 y[5] * y[1] * x[5] * z[0] - y[6] * z[2] * x[3] * y[7] -
790 y[7] * y[6] * x[7] * z[2] + 2.0 * y[6] * y[6] * x[2] * z[7] +
791 y[6] * y[6] * x[3] * z[7] + x[6] * y[7] * y[7] * z[2] -
792 z[6] * x[2] * y[7] * y[7];
793 s8 = -x[2] * y[1] * y[1] * z[3] + 2.0 * y[1] * y[1] * x[0] * z[2] -
794 2.0 * y[1] * y[1] * x[2] * z[0] + z[2] * x[3] * y[1] * y[1] -
795 z[1] * x[0] * y[2] * y[2] + x[1] * y[2] * y[2] * z[0] +
796 y[2] * y[2] * x[0] * z[3] - y[2] * y[2] * x[3] * z[0] -
797 2.0 * y[2] * y[2] * x[3] * z[1] + y[1] * x[1] * y[3] * z[0] -
798 2.0 * y[6] * y[6] * x[7] * z[2] + 2.0 * y[5] * y[5] * x[4] * z[1] -
799 2.0 * y[5] * y[5] * x[1] * z[4] - y[6] * y[6] * x[7] * z[3] -
800 2.0 * y[1] * x[1] * y[0] * z[2];
801 s6 = s8 + 2.0 * y[1] * z[1] * x[2] * y[0] -
802 2.0 * y[1] * z[1] * x[0] * y[2] + 2.0 * y[1] * x[1] * y[2] * z[0] +
803 y[1] * x[2] * y[3] * z[1] - y[1] * y[2] * x[3] * z[1] -
804 y[1] * z[2] * x[1] * y[3] + y[1] * y[2] * x[1] * z[3] -
805 y[2] * x[1] * y[0] * z[2] + y[2] * z[1] * x[2] * y[0] +
806 y[2] * x[2] * y[3] * z[0] - y[7] * x[6] * y[2] * z[7] +
807 y[7] * z[6] * x[7] * y[2] + y[7] * y[6] * x[2] * z[7] -
808 y[6] * x[6] * y[3] * z[7] + y[6] * x[6] * y[7] * z[3] + s7;
809 s8 = s6 - y[6] * z[6] * x[3] * y[7] + y[6] * z[6] * x[7] * y[3] +
810 2.0 * y[2] * y[2] * x[1] * z[3] + x[2] * y[3] * y[3] * z[1] -
811 z[2] * x[1] * y[3] * y[3] + y[1] * x[1] * y[0] * z[4] +
812 y[1] * z[1] * x[3] * y[0] - y[1] * x[1] * y[0] * z[3] +
813 2.0 * y[5] * x[5] * y[1] * z[4] - 2.0 * y[5] * x[5] * y[4] * z[1] +
814 2.0 * y[5] * z[5] * x[1] * y[4] - 2.0 * y[5] * z[5] * x[4] * y[1] -
815 2.0 * y[6] * x[6] * y[2] * z[7] + 2.0 * y[6] * x[6] * y[7] * z[2];
816 s7 = s8 + 2.0 * y[6] * z[6] * x[7] * y[2] -
817 2.0 * y[6] * z[6] * x[2] * y[7] - y[1] * z[1] * x[4] * y[0] +
818 y[1] * z[1] * x[0] * y[4] - y[1] * z[1] * x[0] * y[3] +
819 2.0 * y[6] * y[6] * x[7] * z[5] + 2.0 * y[5] * y[5] * x[6] * z[4] -
820 2.0 * y[5] * y[5] * x[4] * z[6] + x[6] * y[5] * y[5] * z[7] -
821 y[3] * x[2] * y[1] * z[3] - y[3] * y[2] * x[3] * z[1] +
822 y[3] * z[2] * x[3] * y[1] + y[3] * y[2] * x[1] * z[3] -
823 y[2] * x[2] * y[0] * z[3] + y[2] * z[2] * x[3] * y[0];
824 s8 = s7 + 2.0 * y[2] * x[2] * y[3] * z[1] -
825 2.0 * y[2] * x[2] * y[1] * z[3] + y[2] * y[1] * x[0] * z[2] -
826 y[2] * y[1] * x[2] * z[0] + 2.0 * y[2] * z[2] * x[3] * y[1] -
827 2.0 * y[2] * z[2] * x[1] * y[3] - y[2] * z[2] * x[0] * y[3] +
828 y[5] * z[6] * x[5] * y[7] - y[5] * x[6] * y[7] * z[5] -
829 y[5] * y[6] * x[4] * z[7] - y[5] * y[6] * x[5] * z[7] -
830 2.0 * y[5] * x[5] * y[6] * z[4] + 2.0 * y[5] * x[5] * y[4] * z[6] -
831 2.0 * y[5] * z[5] * x[6] * y[4] + 2.0 * y[5] * z[5] * x[4] * y[6];
832 s5 = s8 - y[1] * y[5] * x[0] * z[4] - z[6] * x[7] * y[5] * y[5] +
833 y[6] * y[6] * x[7] * z[4] - y[6] * y[6] * x[4] * z[7] -
834 2.0 * y[6] * y[6] * x[5] * z[7] - x[5] * y[6] * y[6] * z[4] +
835 z[5] * x[4] * y[6] * y[6] + z[6] * x[5] * y[7] * y[7] -
836 x[6] * y[7] * y[7] * z[5] + y[1] * y[5] * x[4] * z[0] +
837 y[7] * y[6] * x[7] * z[5] + y[6] * y[5] * x[7] * z[4] +
838 y[5] * y[6] * x[7] * z[5] + y[6] * y[5] * x[6] * z[4] -
839 y[6] * y[5] * x[4] * z[6] + 2.0 * y[6] * z[6] * x[5] * y[7];
840 s8 = s5 - 2.0 * y[6] * x[6] * y[7] * z[5] +
841 2.0 * y[6] * x[6] * y[5] * z[7] - 2.0 * y[6] * z[6] * x[7] * y[5] -
842 y[6] * x[5] * y[7] * z[4] - y[6] * x[6] * y[7] * z[4] +
843 y[6] * x[6] * y[4] * z[7] - y[6] * z[6] * x[7] * y[4] +
844 y[6] * z[5] * x[4] * y[7] + y[6] * z[6] * x[4] * y[7] +
845 y[6] * x[5] * y[4] * z[6] - y[6] * z[5] * x[6] * y[4] +
846 y[7] * x[6] * y[5] * z[7] - y[7] * z[6] * x[7] * y[5] -
847 2.0 * y[6] * x[6] * y[5] * z[2];
848 s7 = s8 - y[7] * y[6] * x[5] * z[7] + 2.0 * y[4] * y[5] * x[4] * z[0] +
849 2.0 * x[3] * y[7] * y[7] * z[4] - 2.0 * x[4] * y[7] * y[7] * z[3] -
850 z[0] * x[4] * y[7] * y[7] + x[0] * y[7] * y[7] * z[4] -
851 y[0] * z[5] * x[4] * y[1] + y[0] * x[5] * y[1] * z[4] -
852 y[0] * x[5] * y[4] * z[0] + y[0] * z[5] * x[0] * y[4] -
853 y[5] * y[5] * x[0] * z[4] + y[5] * y[5] * x[4] * z[0] +
854 2.0 * y[1] * y[1] * x[2] * z[5] - 2.0 * y[1] * y[1] * x[5] * z[2] +
855 z[1] * x[5] * y[2] * y[2];
856 s8 = s7 - x[1] * y[2] * y[2] * z[5] - y[5] * z[5] * x[4] * y[0] +
857 y[5] * z[5] * x[0] * y[4] - y[5] * x[5] * y[4] * z[0] -
858 y[2] * x[1] * y[6] * z[5] - y[2] * y[1] * x[5] * z[6] +
859 y[2] * z[1] * x[5] * y[6] + y[2] * y[1] * x[6] * z[5] -
860 y[1] * z[1] * x[6] * y[5] - y[1] * x[1] * y[6] * z[5] +
861 y[1] * x[1] * y[5] * z[6] + y[1] * z[1] * x[5] * y[6] +
862 y[5] * x[5] * y[0] * z[4] + y[2] * y[1] * x[2] * z[5] -
863 y[2] * z[1] * x[2] * y[5];
864 s6 = s8 + y[2] * x[1] * y[5] * z[2] - y[2] * y[1] * x[5] * z[2] -
865 y[1] * y[1] * x[5] * z[6] + y[1] * y[1] * x[6] * z[5] -
866 z[1] * x[2] * y[5] * y[5] + x[1] * y[5] * y[5] * z[2] +
867 2.0 * y[1] * z[1] * x[5] * y[2] - 2.0 * y[1] * x[1] * y[2] * z[5] -
868 2.0 * y[1] * z[1] * x[2] * y[5] + 2.0 * y[1] * x[1] * y[5] * z[2] -
869 y[1] * y[1] * x[6] * z[2] + y[1] * y[1] * x[2] * z[6] -
870 2.0 * y[5] * x[1] * y[6] * z[5] - 2.0 * y[5] * y[1] * x[5] * z[6] +
871 2.0 * y[5] * z[1] * x[5] * y[6] + 2.0 * y[5] * y[1] * x[6] * z[5];
872 s8 = s6 - y[6] * z[1] * x[6] * y[5] - y[6] * y[1] * x[5] * z[6] +
873 y[6] * x[1] * y[5] * z[6] + y[6] * y[1] * x[6] * z[5] -
874 2.0 * z[1] * x[6] * y[5] * y[5] + 2.0 * x[1] * y[5] * y[5] * z[6] -
875 x[1] * y[6] * y[6] * z[5] + z[1] * x[5] * y[6] * y[6] +
876 y[5] * z[1] * x[5] * y[2] - y[5] * x[1] * y[2] * z[5] +
877 y[5] * y[1] * x[2] * z[5] - y[5] * y[1] * x[5] * z[2] -
878 y[6] * z[1] * x[2] * y[5] + y[6] * x[1] * y[5] * z[2];
879 s7 = s8 - y[1] * z[1] * x[2] * y[6] - y[1] * x[1] * y[2] * z[6] +
880 y[1] * x[1] * y[6] * z[2] + y[1] * z[1] * x[6] * y[2] +
881 y[5] * x[5] * y[6] * z[2] - y[5] * x[2] * y[6] * z[5] +
882 y[5] * x[6] * y[2] * z[5] - y[5] * x[5] * y[2] * z[6] -
883 x[6] * y[5] * y[5] * z[2] + x[2] * y[5] * y[5] * z[6] -
884 y[5] * y[5] * x[4] * z[7] + y[5] * y[5] * x[7] * z[4] -
885 y[1] * x[6] * y[5] * z[2] + y[1] * x[2] * y[5] * z[6] -
886 y[2] * x[6] * y[5] * z[2] - 2.0 * y[2] * y[1] * x[6] * z[2];
887 s8 = s7 - 2.0 * y[2] * z[1] * x[2] * y[6] +
888 2.0 * y[2] * x[1] * y[6] * z[2] + 2.0 * y[2] * y[1] * x[2] * z[6] -
889 2.0 * x[1] * y[2] * y[2] * z[6] + 2.0 * z[1] * x[6] * y[2] * y[2] +
890 x[6] * y[2] * y[2] * z[5] - x[5] * y[2] * y[2] * z[6] +
891 2.0 * x[5] * y[6] * y[6] * z[2] - 2.0 * x[2] * y[6] * y[6] * z[5] -
892 z[1] * x[2] * y[6] * y[6] - y[6] * y[1] * x[6] * z[2] -
893 y[6] * x[1] * y[2] * z[6] + y[6] * z[1] * x[6] * y[2] +
894 y[6] * y[1] * x[2] * z[6] + x[1] * y[6] * y[6] * z[2];
895 s3 = s8 + y[2] * x[5] * y[6] * z[2] + y[2] * x[2] * y[5] * z[6] -
896 y[2] * x[2] * y[6] * z[5] + y[5] * z[5] * x[4] * y[7] +
897 y[5] * x[5] * y[4] * z[7] - y[5] * z[5] * x[7] * y[4] -
898 y[5] * x[5] * y[7] * z[4] + 2.0 * y[4] * x[5] * y[0] * z[4] -
899 y[3] * z[6] * x[3] * y[7] + y[3] * y[6] * x[3] * z[7] +
900 y[3] * x[6] * y[7] * z[3] - y[3] * y[6] * x[7] * z[3] -
901 y[2] * y[1] * x[3] * z[0] - y[2] * z[1] * x[0] * y[3] +
902 y[2] * y[1] * x[0] * z[3] + y[2] * x[1] * y[3] * z[0];
903 s8 = y[1] * x[0] * z[3] + x[1] * y[3] * z[0] - y[0] * x[3] * z[7] -
904 x[1] * y[5] * z[0] - y[0] * x[3] * z[4] - x[1] * y[0] * z[2] +
905 z[1] * x[2] * y[0] - y[1] * x[0] * z[5] - z[1] * x[0] * y[2] -
906 y[1] * x[0] * z[4] + z[1] * x[5] * y[2] + z[0] * x[7] * y[4] +
907 z[0] * x[3] * y[7] + z[1] * x[0] * y[4] - x[1] * y[2] * z[5] +
908 x[2] * y[3] * z[0] + y[1] * x[2] * z[5] - x[2] * y[3] * z[7];
909 s7 = s8 - z[1] * x[2] * y[5] - y[1] * x[3] * z[0] - x[0] * y[7] * z[3] -
910 z[1] * x[0] * y[3] + y[5] * x[4] * z[0] - x[0] * y[4] * z[3] +
911 y[5] * x[7] * z[4] - z[0] * x[4] * y[3] + x[1] * y[0] * z[4] -
912 z[2] * x[3] * y[7] - y[6] * x[7] * z[2] + x[1] * y[5] * z[2] +
913 y[6] * x[7] * z[5] + x[0] * y[7] * z[4] + x[1] * y[2] * z[0] -
914 z[1] * x[4] * y[0] - z[0] * x[4] * y[7] - z[2] * x[0] * y[3];
915 s8 = x[5] * y[0] * z[4] + z[1] * x[0] * y[5] - x[2] * y[0] * z[3] -
916 z[1] * x[5] * y[0] + y[1] * x[5] * z[0] - x[1] * y[0] * z[3] -
917 x[1] * y[4] * z[0] - y[1] * x[5] * z[2] + x[2] * y[7] * z[3] +
918 y[0] * x[4] * z[3] - x[0] * y[4] * z[7] + x[1] * y[0] * z[5] -
919 y[1] * x[6] * z[2] - y[2] * x[6] * z[3] + y[0] * x[7] * z[3] -
920 y[2] * x[7] * z[3] + z[2] * x[7] * y[3] + y[2] * x[0] * z[3];
921 s6 = s8 + y[2] * x[3] * z[7] - y[2] * x[3] * z[0] - x[6] * y[5] * z[2] -
922 y[5] * x[0] * z[4] + z[2] * x[3] * y[0] + x[2] * y[3] * z[1] +
923 x[0] * y[3] * z[7] - x[2] * y[1] * z[3] + y[1] * x[4] * z[0] +
924 y[1] * x[0] * z[2] - z[1] * x[2] * y[6] + y[2] * x[3] * z[6] -
925 y[1] * x[2] * z[0] + z[1] * x[3] * y[0] - x[1] * y[2] * z[6] -
926 x[2] * y[3] * z[6] + x[0] * y[3] * z[4] + z[0] * x[3] * y[4] + s7;
927 s8 = x[5] * y[4] * z[7] + s6 + y[5] * x[6] * z[4] - y[5] * x[4] * z[6] +
928 z[6] * x[5] * y[7] - x[6] * y[2] * z[7] - x[6] * y[7] * z[5] +
929 x[5] * y[6] * z[2] + x[6] * y[5] * z[7] + x[6] * y[7] * z[2] +
930 y[6] * x[7] * z[4] - y[6] * x[4] * z[7] - y[6] * x[7] * z[3] +
931 z[6] * x[7] * y[2] + x[2] * y[5] * z[6] - x[2] * y[6] * z[5] +
932 y[6] * x[2] * z[7] + x[6] * y[2] * z[5];
933 s7 = s8 - x[5] * y[2] * z[6] - z[6] * x[7] * y[5] - z[5] * x[7] * y[4] +
934 z[5] * x[0] * y[4] - y[5] * x[4] * z[7] + y[0] * x[4] * z[7] -
935 z[6] * x[2] * y[7] - x[5] * y[4] * z[0] - x[5] * y[7] * z[4] -
936 y[0] * x[7] * z[4] + y[5] * x[4] * z[1] - x[6] * y[7] * z[4] +
937 x[7] * y[4] * z[3] - x[4] * y[7] * z[3] + x[3] * y[7] * z[4] -
938 x[7] * y[3] * z[4] - x[6] * y[3] * z[7] + x[6] * y[4] * z[7];
939 s8 = -x[3] * y[4] * z[7] + x[4] * y[3] * z[7] - z[6] * x[7] * y[4] -
940 z[1] * x[6] * y[5] + x[6] * y[7] * z[3] - x[1] * y[6] * z[5] -
941 y[1] * x[5] * z[6] + z[5] * x[4] * y[7] - z[5] * x[4] * y[0] +
942 x[1] * y[5] * z[6] - y[6] * x[5] * z[7] - y[2] * x[3] * z[1] +
943 z[1] * x[5] * y[6] - y[5] * x[1] * z[4] + z[6] * x[4] * y[7] +
944 x[5] * y[1] * z[4] - x[5] * y[6] * z[4] + y[6] * x[3] * z[7] -
946 s5 = s8 + x[5] * y[4] * z[6] + z[5] * x[1] * y[4] + y[1] * x[6] * z[5] -
947 z[6] * x[3] * y[7] + z[6] * x[7] * y[3] - z[5] * x[6] * y[4] -
948 z[5] * x[4] * y[1] + z[5] * x[4] * y[6] + x[1] * y[6] * z[2] +
949 x[2] * y[6] * z[3] + z[2] * x[6] * y[3] + z[1] * x[6] * y[2] +
950 z[2] * x[3] * y[1] - z[2] * x[1] * y[3] - z[2] * x[3] * y[6] +
951 y[2] * x[1] * z[3] + y[1] * x[2] * z[6] - z[0] * x[7] * y[3] + s7;
954 const double unknown1 = s1 * s2;
956 s8 = -z[2] * x[1] * y[2] * z[5] + z[2] * y[1] * x[2] * z[5] -
957 z[2] * z[1] * x[2] * y[5] + z[2] * z[1] * x[5] * y[2] +
958 2.0 * y[5] * x[7] * z[4] * z[4] - y[1] * x[2] * z[0] * z[0] +
959 x[0] * y[3] * z[7] * z[7] - 2.0 * z[5] * z[5] * x[4] * y[1] +
960 2.0 * z[5] * z[5] * x[1] * y[4] + z[5] * z[5] * x[0] * y[4] -
961 2.0 * z[2] * z[2] * x[1] * y[3] + 2.0 * z[2] * z[2] * x[3] * y[1] -
962 x[0] * y[4] * z[7] * z[7] - y[0] * x[3] * z[7] * z[7] +
963 x[1] * y[0] * z[5] * z[5];
964 s7 = s8 - y[1] * x[0] * z[5] * z[5] + z[1] * y[1] * x[2] * z[6] +
965 y[1] * x[0] * z[2] * z[2] + z[2] * z[2] * x[3] * y[0] -
966 z[2] * z[2] * x[0] * y[3] - x[1] * y[0] * z[2] * z[2] +
967 2.0 * z[5] * z[5] * x[4] * y[6] - 2.0 * z[5] * z[5] * x[6] * y[4] -
968 z[5] * z[5] * x[7] * y[4] - x[6] * y[7] * z[5] * z[5] +
969 2.0 * z[2] * y[1] * x[2] * z[6] - 2.0 * z[2] * x[1] * y[2] * z[6] +
970 2.0 * z[2] * z[1] * x[6] * y[2] - y[6] * x[5] * z[7] * z[7] +
971 2.0 * x[6] * y[4] * z[7] * z[7];
972 s8 = -2.0 * y[6] * x[4] * z[7] * z[7] + x[6] * y[5] * z[7] * z[7] -
973 2.0 * z[2] * z[1] * x[2] * y[6] + z[4] * y[6] * x[7] * z[5] +
974 x[5] * y[4] * z[6] * z[6] + z[6] * z[6] * x[4] * y[7] -
975 z[6] * z[6] * x[7] * y[4] - 2.0 * z[6] * z[6] * x[7] * y[5] +
976 2.0 * z[6] * z[6] * x[5] * y[7] - y[5] * x[4] * z[6] * z[6] +
977 2.0 * z[0] * z[0] * x[3] * y[4] - x[6] * y[5] * z[2] * z[2] +
978 z[1] * z[1] * x[5] * y[6] - z[1] * z[1] * x[6] * y[5] -
979 z[5] * z[5] * x[4] * y[0];
980 s6 = s8 + 2.0 * x[1] * y[3] * z[0] * z[0] +
981 2.0 * x[1] * y[6] * z[2] * z[2] - 2.0 * y[1] * x[6] * z[2] * z[2] -
982 y[1] * x[5] * z[2] * z[2] - z[1] * z[1] * x[2] * y[6] -
983 2.0 * z[1] * z[1] * x[2] * y[5] + 2.0 * z[1] * z[1] * x[5] * y[2] +
984 z[1] * y[1] * x[6] * z[5] + y[1] * x[2] * z[5] * z[5] +
985 z[2] * z[1] * x[2] * y[0] + z[1] * x[1] * y[5] * z[6] -
986 z[1] * x[1] * y[6] * z[5] - z[1] * y[1] * x[5] * z[6] -
987 z[1] * x[2] * y[6] * z[5] + z[1] * x[6] * y[2] * z[5] + s7;
988 s8 = -x[1] * y[2] * z[5] * z[5] + z[1] * x[5] * y[6] * z[2] -
989 2.0 * z[2] * z[2] * x[3] * y[6] + 2.0 * z[2] * z[2] * x[6] * y[3] +
990 z[2] * z[2] * x[7] * y[3] - z[2] * z[2] * x[3] * y[7] -
991 z[1] * x[6] * y[5] * z[2] + 2.0 * z[1] * x[1] * y[5] * z[2] -
992 2.0 * x[3] * y[4] * z[7] * z[7] + 2.0 * x[4] * y[3] * z[7] * z[7] +
993 x[5] * y[6] * z[2] * z[2] + y[1] * x[2] * z[6] * z[6] +
994 y[0] * x[4] * z[7] * z[7] + z[2] * x[2] * y[3] * z[0] -
995 x[1] * y[2] * z[6] * z[6];
996 s7 = s8 - z[7] * z[2] * x[3] * y[7] + x[2] * y[6] * z[3] * z[3] -
997 y[2] * x[6] * z[3] * z[3] - z[6] * x[2] * y[3] * z[7] -
998 z[2] * z[1] * x[0] * y[2] + z[6] * z[2] * x[6] * y[3] -
999 z[6] * z[2] * x[3] * y[6] + z[6] * x[2] * y[6] * z[3] +
1000 z[2] * x[1] * y[2] * z[0] + z[6] * y[2] * x[3] * z[7] -
1001 z[4] * z[5] * x[6] * y[4] + z[4] * z[5] * x[4] * y[6] -
1002 z[4] * y[6] * x[5] * z[7] + z[4] * z[6] * x[4] * y[7] +
1003 z[4] * x[5] * y[4] * z[6];
1004 s8 = -z[6] * y[2] * x[6] * z[3] - z[4] * y[5] * x[4] * z[6] -
1005 z[2] * y[1] * x[5] * z[6] + z[2] * x[1] * y[5] * z[6] +
1006 z[4] * x[6] * y[4] * z[7] + 2.0 * z[4] * z[5] * x[4] * y[7] -
1007 z[4] * z[6] * x[7] * y[4] + x[6] * y[7] * z[3] * z[3] -
1008 2.0 * z[4] * z[5] * x[7] * y[4] - 2.0 * z[4] * y[5] * x[4] * z[7] -
1009 z[4] * y[6] * x[4] * z[7] + z[4] * x[6] * y[5] * z[7] -
1010 z[4] * x[6] * y[7] * z[5] + 2.0 * z[4] * x[5] * y[4] * z[7] +
1011 z[2] * x[2] * y[5] * z[6] - z[2] * x[2] * y[6] * z[5];
1012 s5 = s8 + z[2] * x[6] * y[2] * z[5] - z[2] * x[5] * y[2] * z[6] -
1013 z[2] * x[2] * y[3] * z[7] - x[2] * y[3] * z[7] * z[7] +
1014 2.0 * z[2] * x[2] * y[3] * z[1] - z[2] * y[2] * x[3] * z[0] +
1015 z[2] * y[2] * x[0] * z[3] - z[2] * x[2] * y[0] * z[3] -
1016 z[7] * y[2] * x[7] * z[3] + z[7] * z[2] * x[7] * y[3] +
1017 z[7] * x[2] * y[7] * z[3] + z[6] * y[1] * x[2] * z[5] -
1018 z[6] * x[1] * y[2] * z[5] + z[5] * x[1] * y[5] * z[2] + s6 + s7;
1019 s8 = z[5] * z[1] * x[5] * y[2] - z[5] * z[1] * x[2] * y[5] -
1020 y[6] * x[7] * z[2] * z[2] + 2.0 * z[2] * x[2] * y[6] * z[3] -
1021 2.0 * z[2] * x[2] * y[3] * z[6] + 2.0 * z[2] * y[2] * x[3] * z[6] +
1022 y[2] * x[3] * z[6] * z[6] + y[6] * x[7] * z[5] * z[5] +
1023 z[2] * y[2] * x[3] * z[7] - z[2] * y[2] * x[7] * z[3] -
1024 2.0 * z[2] * y[2] * x[6] * z[3] + z[2] * x[2] * y[7] * z[3] +
1025 x[6] * y[2] * z[5] * z[5] - 2.0 * z[2] * x[2] * y[1] * z[3] -
1026 x[2] * y[6] * z[5] * z[5];
1027 s7 = s8 - y[1] * x[5] * z[6] * z[6] + z[6] * x[1] * y[6] * z[2] -
1028 z[3] * z[2] * x[3] * y[6] + z[6] * z[1] * x[6] * y[2] -
1029 z[6] * z[1] * x[2] * y[6] - z[6] * y[1] * x[6] * z[2] -
1030 2.0 * x[5] * y[2] * z[6] * z[6] + z[4] * z[1] * x[0] * y[4] -
1031 z[3] * x[2] * y[3] * z[6] - z[5] * y[1] * x[5] * z[2] +
1032 z[3] * y[2] * x[3] * z[6] + 2.0 * x[2] * y[5] * z[6] * z[6] -
1033 z[5] * x[1] * y[5] * z[0] + y[2] * x[3] * z[7] * z[7] -
1034 x[2] * y[3] * z[6] * z[6];
1035 s8 = z[5] * y[5] * x[4] * z[0] + z[3] * z[2] * x[6] * y[3] +
1036 x[1] * y[5] * z[6] * z[6] + z[5] * y[5] * x[7] * z[4] -
1037 z[1] * x[1] * y[2] * z[6] + z[1] * x[1] * y[6] * z[2] +
1038 2.0 * z[6] * y[6] * x[7] * z[5] - z[7] * y[6] * x[7] * z[2] -
1039 z[3] * y[6] * x[7] * z[2] + x[6] * y[7] * z[2] * z[2] -
1040 2.0 * z[6] * y[6] * x[7] * z[2] - 2.0 * x[6] * y[3] * z[7] * z[7] -
1041 x[6] * y[2] * z[7] * z[7] - z[5] * x[6] * y[5] * z[2] +
1042 y[6] * x[2] * z[7] * z[7];
1043 s6 = s8 + 2.0 * y[6] * x[3] * z[7] * z[7] + z[6] * z[6] * x[7] * y[3] -
1044 y[6] * x[7] * z[3] * z[3] + z[5] * x[5] * y[0] * z[4] +
1045 2.0 * z[6] * z[6] * x[7] * y[2] - 2.0 * z[6] * z[6] * x[2] * y[7] -
1046 z[6] * z[6] * x[3] * y[7] + z[7] * y[6] * x[7] * z[5] +
1047 z[7] * y[5] * x[7] * z[4] - 2.0 * z[7] * x[7] * y[3] * z[4] +
1048 2.0 * z[7] * x[3] * y[7] * z[4] - 2.0 * z[7] * x[4] * y[7] * z[3] +
1049 2.0 * z[7] * x[7] * y[4] * z[3] - z[7] * y[0] * x[7] * z[4] -
1050 2.0 * z[7] * z[6] * x[3] * y[7] + s7;
1051 s8 = s6 + 2.0 * z[7] * z[6] * x[7] * y[3] +
1052 2.0 * z[7] * x[6] * y[7] * z[3] + z[7] * x[6] * y[7] * z[2] -
1053 2.0 * z[7] * y[6] * x[7] * z[3] + z[7] * z[6] * x[7] * y[2] -
1054 z[7] * z[6] * x[2] * y[7] + z[5] * y[1] * x[5] * z[0] -
1055 z[5] * z[1] * x[5] * y[0] + 2.0 * y[1] * x[6] * z[5] * z[5] -
1056 2.0 * x[1] * y[6] * z[5] * z[5] + z[5] * z[1] * x[0] * y[5] +
1057 z[6] * y[6] * x[3] * z[7] + 2.0 * z[6] * x[6] * y[7] * z[2] -
1058 z[6] * y[6] * x[7] * z[3];
1059 s7 = s8 + 2.0 * z[6] * y[6] * x[2] * z[7] - z[6] * x[6] * y[3] * z[7] +
1060 z[6] * x[6] * y[7] * z[3] - 2.0 * z[6] * x[6] * y[2] * z[7] -
1061 2.0 * z[1] * y[1] * x[5] * z[2] - z[1] * y[1] * x[6] * z[2] -
1062 z[7] * z[0] * x[7] * y[3] - 2.0 * z[6] * x[6] * y[5] * z[2] -
1063 z[2] * z[6] * x[3] * y[7] + z[2] * x[6] * y[7] * z[3] -
1064 z[2] * z[6] * x[2] * y[7] + y[5] * x[6] * z[4] * z[4] +
1065 z[2] * y[6] * x[2] * z[7] + y[6] * x[7] * z[4] * z[4] +
1066 z[2] * z[6] * x[7] * y[2] - 2.0 * x[5] * y[7] * z[4] * z[4];
1067 s8 = -x[6] * y[7] * z[4] * z[4] - z[5] * y[5] * x[0] * z[4] -
1068 z[2] * x[6] * y[2] * z[7] - x[5] * y[6] * z[4] * z[4] -
1069 2.0 * z[5] * y[1] * x[5] * z[6] + 2.0 * z[5] * z[1] * x[5] * y[6] +
1070 2.0 * z[5] * x[1] * y[5] * z[6] - 2.0 * z[5] * z[1] * x[6] * y[5] -
1071 z[5] * x[5] * y[2] * z[6] + z[5] * x[5] * y[6] * z[2] +
1072 z[5] * x[2] * y[5] * z[6] + z[5] * z[5] * x[4] * y[7] -
1073 y[5] * x[4] * z[7] * z[7] + x[5] * y[4] * z[7] * z[7] +
1074 z[6] * z[1] * x[5] * y[6] + z[6] * y[1] * x[6] * z[5];
1075 s4 = s8 - z[6] * z[1] * x[6] * y[5] - z[6] * x[1] * y[6] * z[5] +
1076 z[2] * z[6] * x[7] * y[3] + 2.0 * z[6] * x[6] * y[2] * z[5] +
1077 2.0 * z[6] * x[5] * y[6] * z[2] - 2.0 * z[6] * x[2] * y[6] * z[5] +
1078 z[7] * z[0] * x[3] * y[7] + z[7] * z[0] * x[7] * y[4] +
1079 z[3] * z[6] * x[7] * y[3] - z[3] * z[6] * x[3] * y[7] -
1080 z[3] * x[6] * y[3] * z[7] + z[3] * y[6] * x[2] * z[7] -
1081 z[3] * x[6] * y[2] * z[7] + z[5] * x[5] * y[4] * z[7] + s5 + s7;
1082 s8 = s4 + z[3] * y[6] * x[3] * z[7] - z[7] * x[0] * y[7] * z[3] +
1083 z[6] * x[5] * y[4] * z[7] + z[7] * y[0] * x[7] * z[3] +
1084 z[5] * z[6] * x[4] * y[7] - 2.0 * z[5] * x[5] * y[6] * z[4] +
1085 2.0 * z[5] * x[5] * y[4] * z[6] - z[5] * x[5] * y[7] * z[4] -
1086 z[5] * y[6] * x[5] * z[7] - z[5] * z[6] * x[7] * y[4] -
1087 z[7] * z[0] * x[4] * y[7] - z[5] * z[6] * x[7] * y[5] -
1088 z[5] * y[5] * x[4] * z[7] + z[7] * x[0] * y[7] * z[4];
1089 s7 = s8 - 2.0 * z[5] * y[5] * x[4] * z[6] + z[5] * z[6] * x[5] * y[7] +
1090 z[5] * x[6] * y[5] * z[7] + 2.0 * z[5] * y[5] * x[6] * z[4] +
1091 z[6] * z[5] * x[4] * y[6] - z[6] * x[5] * y[6] * z[4] -
1092 z[6] * z[5] * x[6] * y[4] - z[6] * x[6] * y[7] * z[4] -
1093 2.0 * z[6] * y[6] * x[5] * z[7] + z[6] * x[6] * y[4] * z[7] -
1094 z[6] * y[5] * x[4] * z[7] - z[6] * y[6] * x[4] * z[7] +
1095 z[6] * y[6] * x[7] * z[4] + z[6] * y[5] * x[6] * z[4] +
1096 2.0 * z[6] * x[6] * y[5] * z[7];
1097 s8 = -2.0 * z[6] * x[6] * y[7] * z[5] - z[2] * y[1] * x[2] * z[0] +
1098 2.0 * z[7] * z[6] * x[4] * y[7] - 2.0 * z[7] * x[6] * y[7] * z[4] -
1099 2.0 * z[7] * z[6] * x[7] * y[4] + z[7] * z[5] * x[4] * y[7] -
1100 z[7] * z[5] * x[7] * y[4] - z[7] * x[5] * y[7] * z[4] +
1101 2.0 * z[7] * y[6] * x[7] * z[4] - z[7] * z[6] * x[7] * y[5] +
1102 z[7] * z[6] * x[5] * y[7] - z[7] * x[6] * y[7] * z[5] +
1103 z[1] * z[1] * x[6] * y[2] + s7 + x[1] * y[5] * z[2] * z[2];
1104 s6 = s8 + 2.0 * z[2] * y[2] * x[1] * z[3] -
1105 2.0 * z[2] * y[2] * x[3] * z[1] - 2.0 * x[1] * y[4] * z[0] * z[0] +
1106 2.0 * y[1] * x[4] * z[0] * z[0] + 2.0 * x[2] * y[7] * z[3] * z[3] -
1107 2.0 * y[2] * x[7] * z[3] * z[3] - x[1] * y[5] * z[0] * z[0] +
1108 z[0] * z[0] * x[7] * y[4] + z[0] * z[0] * x[3] * y[7] +
1109 x[2] * y[3] * z[0] * z[0] - 2.0 * y[1] * x[3] * z[0] * z[0] +
1110 y[5] * x[4] * z[0] * z[0] - 2.0 * z[0] * z[0] * x[4] * y[3] +
1111 x[1] * y[2] * z[0] * z[0] - z[0] * z[0] * x[4] * y[7] +
1112 y[1] * x[5] * z[0] * z[0];
1113 s8 = s6 - y[2] * x[3] * z[0] * z[0] + y[1] * x[0] * z[3] * z[3] -
1114 2.0 * x[0] * y[7] * z[3] * z[3] - x[0] * y[4] * z[3] * z[3] -
1115 2.0 * x[2] * y[0] * z[3] * z[3] - x[1] * y[0] * z[3] * z[3] +
1116 y[0] * x[4] * z[3] * z[3] - 2.0 * z[0] * y[1] * x[0] * z[4] +
1117 2.0 * z[0] * z[1] * x[0] * y[4] + 2.0 * z[0] * x[1] * y[0] * z[4] -
1118 2.0 * z[0] * z[1] * x[4] * y[0] - 2.0 * z[3] * x[2] * y[3] * z[7] -
1119 2.0 * z[3] * z[2] * x[3] * y[7] + 2.0 * z[3] * z[2] * x[7] * y[3];
1120 s7 = s8 + 2.0 * z[3] * y[2] * x[3] * z[7] +
1121 2.0 * z[5] * y[5] * x[4] * z[1] + 2.0 * z[0] * y[1] * x[0] * z[3] -
1122 z[0] * y[0] * x[3] * z[7] - 2.0 * z[0] * y[0] * x[3] * z[4] -
1123 z[0] * x[1] * y[0] * z[2] + z[0] * z[1] * x[2] * y[0] -
1124 z[0] * y[1] * x[0] * z[5] - z[0] * z[1] * x[0] * y[2] -
1125 z[0] * x[0] * y[7] * z[3] - 2.0 * z[0] * z[1] * x[0] * y[3] -
1126 z[5] * x[5] * y[4] * z[0] - 2.0 * z[0] * x[0] * y[4] * z[3] +
1127 z[0] * x[0] * y[7] * z[4] - z[0] * z[2] * x[0] * y[3];
1128 s8 = s7 + z[0] * x[5] * y[0] * z[4] + z[0] * z[1] * x[0] * y[5] -
1129 z[0] * x[2] * y[0] * z[3] - z[0] * z[1] * x[5] * y[0] -
1130 2.0 * z[0] * x[1] * y[0] * z[3] + 2.0 * z[0] * y[0] * x[4] * z[3] -
1131 z[0] * x[0] * y[4] * z[7] + z[0] * x[1] * y[0] * z[5] +
1132 z[0] * y[0] * x[7] * z[3] + z[0] * y[2] * x[0] * z[3] -
1133 z[0] * y[5] * x[0] * z[4] + z[0] * z[2] * x[3] * y[0] +
1134 z[0] * x[2] * y[3] * z[1] + z[0] * x[0] * y[3] * z[7] -
1135 z[0] * x[2] * y[1] * z[3];
1136 s5 = s8 + z[0] * y[1] * x[0] * z[2] + z[3] * x[1] * y[3] * z[0] -
1137 2.0 * z[3] * y[0] * x[3] * z[7] - z[3] * y[0] * x[3] * z[4] -
1138 z[3] * x[1] * y[0] * z[2] + z[3] * z[0] * x[7] * y[4] +
1139 2.0 * z[3] * z[0] * x[3] * y[7] + 2.0 * z[3] * x[2] * y[3] * z[0] -
1140 z[3] * y[1] * x[3] * z[0] - z[3] * z[1] * x[0] * y[3] -
1141 z[3] * z[0] * x[4] * y[3] + z[3] * x[1] * y[2] * z[0] -
1142 z[3] * z[0] * x[4] * y[7] - 2.0 * z[3] * z[2] * x[0] * y[3] -
1143 z[3] * x[0] * y[4] * z[7] - 2.0 * z[3] * y[2] * x[3] * z[0];
1144 s8 = s5 + 2.0 * z[3] * z[2] * x[3] * y[0] + z[3] * x[2] * y[3] * z[1] +
1145 2.0 * z[3] * x[0] * y[3] * z[7] + z[3] * y[1] * x[0] * z[2] -
1146 z[4] * y[0] * x[3] * z[7] - z[4] * x[1] * y[5] * z[0] -
1147 z[4] * y[1] * x[0] * z[5] + 2.0 * z[4] * z[0] * x[7] * y[4] +
1148 z[4] * z[0] * x[3] * y[7] + 2.0 * z[4] * y[5] * x[4] * z[0] +
1149 2.0 * y[0] * x[7] * z[3] * z[3] + 2.0 * y[2] * x[0] * z[3] * z[3] -
1150 x[2] * y[1] * z[3] * z[3] - y[0] * x[3] * z[4] * z[4];
1151 s7 = s8 - y[1] * x[0] * z[4] * z[4] + x[1] * y[0] * z[4] * z[4] +
1152 2.0 * x[0] * y[7] * z[4] * z[4] + 2.0 * x[5] * y[0] * z[4] * z[4] -
1153 2.0 * y[5] * x[0] * z[4] * z[4] + 2.0 * z[1] * z[1] * x[2] * y[0] -
1154 2.0 * z[1] * z[1] * x[0] * y[2] + z[1] * z[1] * x[0] * y[4] -
1155 z[1] * z[1] * x[0] * y[3] - z[1] * z[1] * x[4] * y[0] +
1156 2.0 * z[1] * z[1] * x[0] * y[5] - 2.0 * z[1] * z[1] * x[5] * y[0] +
1157 x[2] * y[3] * z[1] * z[1] - x[5] * y[4] * z[0] * z[0] -
1158 z[0] * z[0] * x[7] * y[3];
1159 s8 = s7 + x[7] * y[4] * z[3] * z[3] - x[4] * y[7] * z[3] * z[3] +
1160 y[2] * x[1] * z[3] * z[3] + x[0] * y[3] * z[4] * z[4] -
1161 2.0 * y[0] * x[7] * z[4] * z[4] + x[3] * y[7] * z[4] * z[4] -
1162 x[7] * y[3] * z[4] * z[4] - y[5] * x[1] * z[4] * z[4] +
1163 x[5] * y[1] * z[4] * z[4] + z[1] * z[1] * x[3] * y[0] +
1164 y[5] * x[4] * z[1] * z[1] - y[2] * x[3] * z[1] * z[1] -
1165 x[5] * y[4] * z[1] * z[1] - z[4] * x[0] * y[4] * z[3] -
1166 z[4] * z[0] * x[4] * y[3];
1167 s6 = s8 - z[4] * z[1] * x[4] * y[0] - 2.0 * z[4] * z[0] * x[4] * y[7] +
1168 z[4] * y[1] * x[5] * z[0] - 2.0 * z[5] * x[5] * y[4] * z[1] -
1169 z[4] * x[1] * y[4] * z[0] + z[4] * y[0] * x[4] * z[3] -
1170 2.0 * z[4] * x[0] * y[4] * z[7] + z[4] * x[1] * y[0] * z[5] -
1171 2.0 * z[1] * x[1] * y[2] * z[5] + z[4] * x[0] * y[3] * z[7] +
1172 2.0 * z[5] * x[5] * y[1] * z[4] + z[4] * y[1] * x[4] * z[0] +
1173 z[1] * y[1] * x[0] * z[3] + z[1] * x[1] * y[3] * z[0] -
1174 2.0 * z[1] * x[1] * y[5] * z[0] - 2.0 * z[1] * x[1] * y[0] * z[2];
1175 s8 = s6 - 2.0 * z[1] * y[1] * x[0] * z[5] - z[1] * y[1] * x[0] * z[4] +
1176 2.0 * z[1] * y[1] * x[2] * z[5] - z[1] * y[1] * x[3] * z[0] -
1177 2.0 * z[5] * y[5] * x[1] * z[4] + z[1] * y[5] * x[4] * z[0] +
1178 z[1] * x[1] * y[0] * z[4] + 2.0 * z[1] * x[1] * y[2] * z[0] -
1179 z[1] * z[2] * x[0] * y[3] + 2.0 * z[1] * y[1] * x[5] * z[0] -
1180 z[1] * x[1] * y[0] * z[3] - z[1] * x[1] * y[4] * z[0] +
1181 2.0 * z[1] * x[1] * y[0] * z[5] - z[1] * y[2] * x[3] * z[0];
1182 s7 = s8 + z[1] * z[2] * x[3] * y[0] - z[1] * x[2] * y[1] * z[3] +
1183 z[1] * y[1] * x[4] * z[0] + 2.0 * z[1] * y[1] * x[0] * z[2] +
1184 2.0 * z[0] * z[1] * x[3] * y[0] + 2.0 * z[0] * x[0] * y[3] * z[4] +
1185 z[0] * z[5] * x[0] * y[4] + z[0] * y[0] * x[4] * z[7] -
1186 z[0] * y[0] * x[7] * z[4] - z[0] * x[7] * y[3] * z[4] -
1187 z[0] * z[5] * x[4] * y[0] - z[0] * x[5] * y[4] * z[1] +
1188 z[3] * z[1] * x[3] * y[0] + z[3] * x[0] * y[3] * z[4] +
1189 z[3] * z[0] * x[3] * y[4] + z[3] * y[0] * x[4] * z[7];
1190 s8 = s7 + z[3] * x[3] * y[7] * z[4] - z[3] * x[7] * y[3] * z[4] -
1191 z[3] * x[3] * y[4] * z[7] + z[3] * x[4] * y[3] * z[7] -
1192 z[3] * y[2] * x[3] * z[1] + z[3] * z[2] * x[3] * y[1] -
1193 z[3] * z[2] * x[1] * y[3] - 2.0 * z[3] * z[0] * x[7] * y[3] +
1194 z[4] * z[0] * x[3] * y[4] + 2.0 * z[4] * z[5] * x[0] * y[4] +
1195 2.0 * z[4] * y[0] * x[4] * z[7] - 2.0 * z[4] * x[5] * y[4] * z[0] +
1196 z[4] * y[5] * x[4] * z[1] + z[4] * x[7] * y[4] * z[3] -
1197 z[4] * x[4] * y[7] * z[3];
1198 s3 = s8 - z[4] * x[3] * y[4] * z[7] + z[4] * x[4] * y[3] * z[7] -
1199 2.0 * z[4] * z[5] * x[4] * y[0] - z[4] * x[5] * y[4] * z[1] +
1200 z[4] * z[5] * x[1] * y[4] - z[4] * z[5] * x[4] * y[1] -
1201 2.0 * z[1] * y[1] * x[2] * z[0] + z[1] * z[5] * x[0] * y[4] -
1202 z[1] * z[5] * x[4] * y[0] - z[1] * y[5] * x[1] * z[4] +
1203 z[1] * x[5] * y[1] * z[4] + z[1] * z[5] * x[1] * y[4] -
1204 z[1] * z[5] * x[4] * y[1] + z[1] * z[2] * x[3] * y[1] -
1205 z[1] * z[2] * x[1] * y[3] + z[1] * y[2] * x[1] * z[3];
1206 s8 = y[1] * x[0] * z[3] + x[1] * y[3] * z[0] - y[0] * x[3] * z[7] -
1207 x[1] * y[5] * z[0] - y[0] * x[3] * z[4] - x[1] * y[0] * z[2] +
1208 z[1] * x[2] * y[0] - y[1] * x[0] * z[5] - z[1] * x[0] * y[2] -
1209 y[1] * x[0] * z[4] + z[1] * x[5] * y[2] + z[0] * x[7] * y[4] +
1210 z[0] * x[3] * y[7] + z[1] * x[0] * y[4] - x[1] * y[2] * z[5] +
1211 x[2] * y[3] * z[0] + y[1] * x[2] * z[5] - x[2] * y[3] * z[7];
1212 s7 = s8 - z[1] * x[2] * y[5] - y[1] * x[3] * z[0] - x[0] * y[7] * z[3] -
1213 z[1] * x[0] * y[3] + y[5] * x[4] * z[0] - x[0] * y[4] * z[3] +
1214 y[5] * x[7] * z[4] - z[0] * x[4] * y[3] + x[1] * y[0] * z[4] -
1215 z[2] * x[3] * y[7] - y[6] * x[7] * z[2] + x[1] * y[5] * z[2] +
1216 y[6] * x[7] * z[5] + x[0] * y[7] * z[4] + x[1] * y[2] * z[0] -
1217 z[1] * x[4] * y[0] - z[0] * x[4] * y[7] - z[2] * x[0] * y[3];
1218 s8 = x[5] * y[0] * z[4] + z[1] * x[0] * y[5] - x[2] * y[0] * z[3] -
1219 z[1] * x[5] * y[0] + y[1] * x[5] * z[0] - x[1] * y[0] * z[3] -
1220 x[1] * y[4] * z[0] - y[1] * x[5] * z[2] + x[2] * y[7] * z[3] +
1221 y[0] * x[4] * z[3] - x[0] * y[4] * z[7] + x[1] * y[0] * z[5] -
1222 y[1] * x[6] * z[2] - y[2] * x[6] * z[3] + y[0] * x[7] * z[3] -
1223 y[2] * x[7] * z[3] + z[2] * x[7] * y[3] + y[2] * x[0] * z[3];
1224 s6 = s8 + y[2] * x[3] * z[7] - y[2] * x[3] * z[0] - x[6] * y[5] * z[2] -
1225 y[5] * x[0] * z[4] + z[2] * x[3] * y[0] + x[2] * y[3] * z[1] +
1226 x[0] * y[3] * z[7] - x[2] * y[1] * z[3] + y[1] * x[4] * z[0] +
1227 y[1] * x[0] * z[2] - z[1] * x[2] * y[6] + y[2] * x[3] * z[6] -
1228 y[1] * x[2] * z[0] + z[1] * x[3] * y[0] - x[1] * y[2] * z[6] -
1229 x[2] * y[3] * z[6] + x[0] * y[3] * z[4] + z[0] * x[3] * y[4] + s7;
1230 s8 = x[5] * y[4] * z[7] + s6 + y[5] * x[6] * z[4] - y[5] * x[4] * z[6] +
1231 z[6] * x[5] * y[7] - x[6] * y[2] * z[7] - x[6] * y[7] * z[5] +
1232 x[5] * y[6] * z[2] + x[6] * y[5] * z[7] + x[6] * y[7] * z[2] +
1233 y[6] * x[7] * z[4] - y[6] * x[4] * z[7] - y[6] * x[7] * z[3] +
1234 z[6] * x[7] * y[2] + x[2] * y[5] * z[6] - x[2] * y[6] * z[5] +
1235 y[6] * x[2] * z[7] + x[6] * y[2] * z[5];
1236 s7 = s8 - x[5] * y[2] * z[6] - z[6] * x[7] * y[5] - z[5] * x[7] * y[4] +
1237 z[5] * x[0] * y[4] - y[5] * x[4] * z[7] + y[0] * x[4] * z[7] -
1238 z[6] * x[2] * y[7] - x[5] * y[4] * z[0] - x[5] * y[7] * z[4] -
1239 y[0] * x[7] * z[4] + y[5] * x[4] * z[1] - x[6] * y[7] * z[4] +
1240 x[7] * y[4] * z[3] - x[4] * y[7] * z[3] + x[3] * y[7] * z[4] -
1241 x[7] * y[3] * z[4] - x[6] * y[3] * z[7] + x[6] * y[4] * z[7];
1242 s8 = -x[3] * y[4] * z[7] + x[4] * y[3] * z[7] - z[6] * x[7] * y[4] -
1243 z[1] * x[6] * y[5] + x[6] * y[7] * z[3] - x[1] * y[6] * z[5] -
1244 y[1] * x[5] * z[6] + z[5] * x[4] * y[7] - z[5] * x[4] * y[0] +
1245 x[1] * y[5] * z[6] - y[6] * x[5] * z[7] - y[2] * x[3] * z[1] +
1246 z[1] * x[5] * y[6] - y[5] * x[1] * z[4] + z[6] * x[4] * y[7] +
1247 x[5] * y[1] * z[4] - x[5] * y[6] * z[4] + y[6] * x[3] * z[7] -
1249 s5 = s8 + x[5] * y[4] * z[6] + z[5] * x[1] * y[4] + y[1] * x[6] * z[5] -
1250 z[6] * x[3] * y[7] + z[6] * x[7] * y[3] - z[5] * x[6] * y[4] -
1251 z[5] * x[4] * y[1] + z[5] * x[4] * y[6] + x[1] * y[6] * z[2] +
1252 x[2] * y[6] * z[3] + z[2] * x[6] * y[3] + z[1] * x[6] * y[2] +
1253 z[2] * x[3] * y[1] - z[2] * x[1] * y[3] - z[2] * x[3] * y[6] +
1254 y[2] * x[1] * z[3] + y[1] * x[2] * z[6] - z[0] * x[7] * y[3] + s7;
1257 const double unknown2 = s1 * s2;
1259 return {unknown0, unknown1, unknown2};
1275 template <
int structdim,
int dim,
int spacedim>
1287 template <
int dim,
int spacedim>
1293 return (accessor.
vertex(1) - accessor.
vertex(0)).norm();
1361 for (
unsigned int d = 0;
d < 3; ++
d)
1363 if (a_max < 1e-14 * accessor.
diameter())
1364 return cross_product_3d(b0, b1).
norm();
1369 constexpr unsigned int n_qp = 4;
1370 const double c1 = 2.0 / 7.0 *
std::sqrt(6.0 / 5.0);
1371 const double w0 = (18.0 -
std::sqrt(30)) / 72.0;
1372 const double w1 = (18.0 +
std::sqrt(30)) / 72.0;
1374 const std::array<double, n_qp> q{{
1380 const std::array<double, n_qp>
w{{w0, w1, w1, w0}};
1383 for (
unsigned int i = 0; i < n_qp; ++i)
1384 for (
unsigned int j = 0; j < n_qp; ++j)
1385 area += cross_product_3d(q[i] * a + b0, q[j] * a + b1).
norm() *
1395 return 0.5 * cross_product_3d(v01, v02).
norm();
1404 template <
int structdim,
int dim,
int spacedim>
1411 return std::numeric_limits<double>::quiet_NaN();
1415 template <
int structdim,
int dim,
int spacedim>
1418 const bool use_interpolation)
1420 if (use_interpolation)
1423 const auto points_and_weights =
1427 points_and_weights.first.end()),
1429 points_and_weights.second.end()));
1434 if constexpr (structdim == 1)
1436 else if constexpr (structdim == 2)
1438 else if constexpr (structdim == 3)
1453template <
int structdim,
int dim,
int spacedim>
1456template <
int structdim,
int dim,
int spacedim>
1459template <
int structdim,
int dim,
int spacedim>
1467template <
int structdim,
int dim,
int spacedim>
1470 const std::initializer_list<int> &new_indices)
const
1473 this->objects().get_bounding_object_indices(this->present_index);
1478 for (
const auto &new_index : new_indices)
1480 bounding_object_index_ref[i] = new_index;
1487template <
int structdim,
int dim,
int spacedim>
1490 const std::initializer_list<unsigned int> &new_indices)
const
1493 this->objects().get_bounding_object_indices(this->present_index);
1498 for (
const auto &new_index : new_indices)
1500 bounding_object_index_ref[i] = new_index;
1507template <
int structdim,
int dim,
int spacedim>
1513 return ::barycenter(*
this);
1518template <
int structdim,
int dim,
int spacedim>
1524 return ::measure(*
this);
1529template <
int structdim,
int dim,
int spacedim>
1534 std::make_pair(this->vertex(0), this->vertex(0));
1536 const unsigned int n_vertices = this->n_vertices();
1537 for (
unsigned int v = 1; v < n_vertices; ++v)
1540 for (
unsigned int k = 0; k < spacedim; ++k)
1542 boundary_points.first[k] =
std::min(boundary_points.first[k], x[k]);
1543 boundary_points.second[k] =
std::max(boundary_points.second[k], x[k]);
1552template <
int structdim,
int dim,
int spacedim>
1555 const unsigned int )
const
1558 return std::numeric_limits<double>::signaling_NaN();
1569 return this->diameter();
1579 return this->diameter();
1587 const unsigned int lines[2][2] = {
1593 return std::max(this->line(lines[axis][0])->diameter(),
1594 this->line(lines[axis][1])->diameter());
1601 const unsigned int lines[2][2] = {
1607 return std::max(this->line(lines[axis][0])->diameter(),
1608 this->line(lines[axis][1])->diameter());
1616 const unsigned int lines[3][4] = {
1623 double lengths[4] = {this->line(lines[axis][0])->diameter(),
1624 this->line(lines[axis][1])->diameter(),
1625 this->line(lines[axis][2])->diameter(),
1626 this->line(lines[axis][3])->diameter()};
1628 return std::max({lengths[0], lengths[1], lengths[2], lengths[3]});
1638 set_manifold_id(manifold_ind);
1640 if (this->has_children())
1641 for (
unsigned int c = 0; c < this->n_children(); ++c)
1642 this->child(c)->set_all_manifold_ids(manifold_ind);
1647 for (
const unsigned int i : this->face_indices())
1648 this->quad(i)->set_manifold_id(manifold_ind);
1649 for (
const unsigned int i : this->line_indices())
1650 this->line(i)->set_manifold_id(manifold_ind);
1654template <
int structdim,
int dim,
int spacedim>
1661 std::array<double, GeometryInfo<structdim>::vertices_per_cell> w;
1665 p[i] = this->vertex(i);
1669 return this->get_manifold().get_new_point(
make_array_view(p.begin(), p.end()),
1676template <
int structdim,
int dim,
int spacedim>
1684 vertices[v] = this->vertex(v);
1687 GridTools::affine_cell_approximation<structdim, spacedim>(vertices);
1695template <
int structdim,
int dim,
int spacedim>
1698 const bool respect_manifold,
1699 const bool use_interpolation)
const
1701 if (respect_manifold ==
false)
1707 return p / this->n_vertices();
1710 return get_new_point_on_object(*
this, use_interpolation);
1647 for (
const unsigned int i : this->face_indices()) {
…}
1717template <
int spacedim>
1728template <
int spacedim>
1738template <
int spacedim>
1748template <
int spacedim>
1754 if (this->has_children())
1755 for (
unsigned int c = 0; c < this->n_children(); ++c)
1756 this->child(c)->recursively_set_user_flag();
1761template <
int spacedim>
1767 if (this->has_children())
1768 for (
unsigned int c = 0; c < this->n_children(); ++c)
1769 this->child(c)->recursively_clear_user_flag();
1774template <
int spacedim>
1784template <
int spacedim>
1794template <
int spacedim>
1804template <
int spacedim>
1815template <
int spacedim>
1819 set_user_pointer(p);
1821 if (this->has_children())
1822 for (
unsigned int c = 0; c < this->n_children(); ++c)
1823 this->child(c)->recursively_set_user_pointer(p);
1828template <
int spacedim>
1832 clear_user_pointer();
1834 if (this->has_children())
1835 for (
unsigned int c = 0; c < this->n_children(); ++c)
1836 this->child(c)->recursively_clear_user_pointer();
1841template <
int spacedim>
1851template <
int spacedim>
1861template <
int spacedim>
1872template <
int spacedim>
1878 if (this->has_children())
1879 for (
unsigned int c = 0; c < this->n_children(); ++c)
1880 this->child(c)->recursively_set_user_index(p);
1885template <
int spacedim>
1891 if (this->has_children())
1892 for (
unsigned int c = 0; c < this->n_children(); ++c)
1893 this->child(c)->recursively_clear_user_index();
1906 return (this->vertex(0)[0] <= p[0]) && (p[0] <= this->vertex(1)[0]);
1939 static const int direction[4] = {-1, 1, 1, -1};
1940 for (
unsigned int f = 0; f < 4; ++f)
1964 if ((-face[1] * to_p[0] + face[0] * to_p[1]) < 0)
1995 const unsigned int dim = 3;
1996 const unsigned int spacedim = 3;
2000 for (
unsigned int v = 1; v < this->n_vertices(); ++v)
2001 for (
unsigned int d = 0; d < dim; ++d)
2003 maxp[d] =
std::max(maxp[d], this->vertex(v)[d]);
2004 minp[d] =
std::min(minp[d], this->vertex(v)[d]);
2009 for (
unsigned int d = 0; d < dim; ++d)
2010 if ((p[d] < minp[d]) || (p[d] > maxp[d]))
2025 .
template get_default_linear_mapping<dim, spacedim>()
2026 .transform_real_to_unit_cell(cell_iterator, p)));
2039template <
int dim,
int spacedim>
2045 ExcMessage(
"The current iterator points to an inactive cell. "
2046 "You cannot convert it to an iterator to an active cell."));
2048 ExcMessage(
"The triangulation associated with the iterator does not "
2049 "match that of the DoFHandler."));
2060template <
int dim,
int spacedim>
2066 ExcMessage(
"The triangulation associated with the iterator does not "
2067 "match that of the DoFHandler."));
2082template <
int dim,
int spacedim>
2083template <
int dim_,
int spacedim_>
2092 this->reference_cell()
2093 .template get_default_linear_mapping<dim_, spacedim_>()
2094 .transform_real_to_unit_cell(cell_iterator, p);
2105 return point_inside_codim<1, 2>(p);
2113 return point_inside_codim<1, 3>(p);
2123 return point_inside_codim<2, 3>(p);
2128template <
int dim,
int spacedim>
2132 for (
const auto face : this->face_indices())
2133 if (at_boundary(face))
2141template <
int dim,
int spacedim>
2146 return this->tria->levels[this->present_level]
2147 ->cells.boundary_or_material_id[this->present_index]
2153template <
int dim,
int spacedim>
2160 this->tria->levels[this->present_level]
2161 ->cells.boundary_or_material_id[this->present_index]
2162 .material_id = mat_id;
2167template <
int dim,
int spacedim>
2172 set_material_id(mat_id);
2174 if (this->has_children())
2175 for (
unsigned int c = 0; c < this->n_children(); ++c)
2176 this->child(c)->recursively_set_material_id(mat_id);
2181template <
int dim,
int spacedim>
2187 Assert(this->is_active(),
2188 ExcMessage(
"set_subdomain_id() can only be called on active cells!"));
2189 this->tria->levels[this->present_level]->subdomain_ids[this->present_index] =
2195template <
int dim,
int spacedim>
2201 this->tria->levels[this->present_level]
2202 ->level_subdomain_ids[this->present_index] = new_level_subdomain_id;
2206template <
int dim,
int spacedim>
2211 if constexpr (dim == spacedim)
2213 else if constexpr (dim == spacedim - 1)
2214 return this->tria->levels[this->present_level]
2215 ->direction_flags[this->present_index];
2219 ExcMessage(
"This function cannot be called if dim<spacedim-1."));
2226template <
int dim,
int spacedim>
2229 const bool new_direction_flag)
const
2235 if constexpr (dim == spacedim)
2236 Assert(new_direction_flag ==
true,
2237 ExcMessage(
"If dim==spacedim, direction flags are always true and "
2238 "can not be set to anything else."));
2239 else if constexpr (dim == spacedim - 1)
2240 this->tria->levels[this->present_level]
2241 ->direction_flags[this->present_index] = new_direction_flag;
2243 Assert(new_direction_flag ==
true,
2244 ExcMessage(
"If dim<spacedim-1, then this function can be called "
2245 "only if the argument is 'true'."));
2250template <
int dim,
int spacedim>
2256 this->tria->levels[this->present_level]->parents[this->present_index / 2] =
2262template <
int dim,
int spacedim>
2271 return this->tria->levels[this->present_level]
2272 ->parents[this->present_index / 2];
2277template <
int dim,
int spacedim>
2280 const unsigned int active_cell_index)
const
2282 this->tria->levels[this->present_level]
2283 ->active_cell_indices[this->present_index] = active_cell_index;
2288template <
int dim,
int spacedim>
2293 this->tria->levels[this->present_level]
2294 ->global_active_cell_indices[this->present_index] = index;
2299template <
int dim,
int spacedim>
2304 this->tria->levels[this->present_level]
2305 ->global_level_cell_indices[this->present_index] = index;
2310template <
int dim,
int spacedim>
2317 this->present_level - 1,
2324template <
int dim,
int spacedim>
2329 if (this->has_children())
2330 for (
unsigned int c = 0; c < this->n_children(); ++c)
2331 this->child(c)->recursively_set_subdomain_id(new_subdomain_id);
2333 set_subdomain_id(new_subdomain_id);
2338template <
int dim,
int spacedim>
2341 const unsigned int i,
2347 this->tria->levels[this->present_level]
2348 ->neighbors[this->present_index * ReferenceCells::max_n_faces<dim>() + i];
2351 neighbor.first = pointer->present_level;
2352 neighbor.second = pointer->present_index;
2356 neighbor.first = -1;
2357 neighbor.second = -1;
2363template <
int dim,
int spacedim>
2367 std::array<unsigned char, 30> id;
2370 const unsigned int n_child_indices = ptr.
level();
2372 while (ptr.
level() > 0)
2375 const unsigned int n_children = parent->n_children();
2378 unsigned char v =
static_cast<unsigned char>(-1);
2379 for (
unsigned int c = 0; c < n_children; ++c)
2381 if (parent->child_index(c) == ptr.
index())
2389 id[ptr.
level() - 1] = v;
2395 const unsigned int coarse_index = ptr.
index();
2397 return {this->tria->coarse_cell_index_to_coarse_cell_id(coarse_index),
2404template <
int dim,
int spacedim>
2407 const unsigned int neighbor)
const
2420 if ((dim == 1) && (spacedim == dim))
2424 this->neighbor(neighbor);
2439 const unsigned int this_face_index = face_index(neighbor);
2441 const unsigned int neighbor_guess =
2444 if (neighbor_guess < neighbor_cell->n_faces() &&
2445 neighbor_cell->face_index(neighbor_guess) == this_face_index)
2446 return neighbor_guess;
2453 for (
const unsigned int face_no : neighbor_cell->face_indices())
2454 if (neighbor_cell->face_index(face_no) == this_face_index)
2470template <
int dim,
int spacedim>
2473 const unsigned int face_no)
const
2475 const unsigned int n2 = neighbor_of_neighbor_internal(face_no);
2484template <
int dim,
int spacedim>
2487 const unsigned int face_no)
const
2489 return neighbor_of_neighbor_internal(face_no) ==
2495template <
int dim,
int spacedim>
2496std::pair<unsigned int, unsigned int>
2498 const unsigned int neighbor)
const
2503 Assert(neighbor_is_coarser(neighbor),
2510 const int this_face_index = face_index(neighbor);
2512 this->neighbor(neighbor);
2527 const unsigned int face_no_guess =
2531 neighbor_cell->face(face_no_guess);
2533 if (face_guess->has_children())
2534 for (
unsigned int subface_no = 0;
2535 subface_no < face_guess->n_children();
2537 if (face_guess->child_index(subface_no) == this_face_index)
2538 return std::make_pair(face_no_guess, subface_no);
2544 for (
const unsigned int face_no : neighbor_cell->face_indices())
2546 if (face_no != face_no_guess)
2549 face = neighbor_cell->face(face_no);
2550 if (face->has_children())
2551 for (
unsigned int subface_no = 0;
2552 subface_no < face->n_children();
2554 if (face->child_index(subface_no) == this_face_index)
2555 return std::make_pair(face_no, subface_no);
2569 const int this_face_index = face_index(neighbor);
2571 this->neighbor(neighbor);
2579 const unsigned int face_no_guess =
2583 neighbor_cell->face(face_no_guess);
2585 if (face_guess->has_children())
2586 for (
unsigned int subface_no = 0;
2587 subface_no < face_guess->n_children();
2590 if (face_guess->child_index(subface_no) == this_face_index)
2594 return std::make_pair(face_no_guess,
2595 translate_subface_no(face_guess,
2598 if (face_guess->child(subface_no)->has_children())
2599 for (
unsigned int subsub_no = 0;
2600 subsub_no < face_guess->child(subface_no)->n_children();
2602 if (face_guess->child(subface_no)->child_index(subsub_no) ==
2607 return std::make_pair(face_no_guess,
2608 translate_subface_no(face_guess,
2615 for (
const unsigned int face_no : neighbor_cell->face_indices())
2617 if (face_no == face_no_guess)
2621 neighbor_cell->face(face_no);
2623 if (!face->has_children())
2626 for (
unsigned int subface_no = 0; subface_no < face->n_children();
2629 if (face->child_index(subface_no) == this_face_index)
2633 return std::make_pair(face_no,
2634 translate_subface_no(face,
2637 if (face->child(subface_no)->has_children())
2638 for (
unsigned int subsub_no = 0;
2639 subsub_no < face->child(subface_no)->n_children();
2641 if (face->child(subface_no)->child_index(subsub_no) ==
2646 return std::make_pair(face_no,
2647 translate_subface_no(face,
2671template <
int dim,
int spacedim>
2674 const unsigned int i_face)
const
2689 cell_iterator current_cell(*
this);
2690 if (this->tria->periodic_face_map.find(
2691 std::make_pair(current_cell, i_face)) !=
2692 this->tria->periodic_face_map.end())
2699template <
int dim,
int spacedim>
2713 cell_iterator current_cell(*
this);
2716 this->tria->periodic_face_map.find(std::make_pair(current_cell, i_face));
2719 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2721 return my_face_pair->second.first.first;
2726template <
int dim,
int spacedim>
2729 const unsigned int i_face)
const
2731 if (!(this->face(i_face)->at_boundary()))
2732 return this->neighbor(i_face);
2733 else if (this->has_periodic_neighbor(i_face))
2734 return this->periodic_neighbor(i_face);
2738 return this->neighbor(i_face);
2743template <
int dim,
int spacedim>
2746 const unsigned int i_face,
2747 const unsigned int i_subface)
const
2763 cell_iterator my_it(*
this);
2766 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2771 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2773 cell_iterator parent_nb_it = my_face_pair->second.first.first;
2774 unsigned int nb_face_num = my_face_pair->second.first.second;
2776 parent_nb_it->face(nb_face_num);
2783 const auto [orientation, rotation, flip] =
2786 unsigned int sub_neighbor_num =
2793 nb_parent_face_it->refinement_case());
2794 return parent_nb_it->child(sub_neighbor_num);
2799template <
int dim,
int spacedim>
2800std::pair<unsigned int, unsigned int>
2802 const unsigned int i_face)
const
2818 const int my_face_index = this->face_index(i_face);
2819 cell_iterator my_it(*
this);
2822 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2827 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2829 cell_iterator nb_it = my_face_pair->second.first.first;
2830 unsigned int face_num_of_nb = my_face_pair->second.first.second;
2833 this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2839 Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2841 cell_iterator p_nb_of_p_nb = nb_face_pair->second.first.first;
2843 p_nb_of_p_nb->face(nb_face_pair->second.first.second);
2844 for (
unsigned int i_subface = 0; i_subface < parent_face_it->n_children();
2846 if (parent_face_it->child_index(i_subface) == my_face_index)
2847 return std::make_pair(face_num_of_nb, i_subface);
2860template <
int dim,
int spacedim>
2863 const unsigned int i_face)
const
2865 return periodic_neighbor(i_face)->index();
2870template <
int dim,
int spacedim>
2873 const unsigned int i_face)
const
2875 return periodic_neighbor(i_face)->level();
2880template <
int dim,
int spacedim>
2883 const unsigned int i_face)
const
2885 return periodic_neighbor_face_no(i_face);
2890template <
int dim,
int spacedim>
2893 const unsigned int i_face)
const
2905 cell_iterator my_it(*
this);
2908 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2914 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2916 return my_face_pair->second.first.second;
2921template <
int dim,
int spacedim>
2924 const unsigned int i_face)
const
2945 cell_iterator my_it(*
this);
2948 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2953 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2956 cell_iterator nb_it = my_face_pair->second.first.first;
2957 unsigned int face_num_of_nb = my_face_pair->second.first.second;
2960 this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2966 Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2968 const unsigned int my_level = this->
level();
2969 const unsigned int neighbor_level = nb_face_pair->second.first.first->level();
2971 return my_level > neighbor_level;
2976template <
int dim,
int spacedim>
2983 return (neighbor_index(i) == -1);
2988template <
int dim,
int spacedim>
2993 return at_boundary();
2996 for (
unsigned int l = 0; l < this->n_lines(); ++l)
2997 if (this->line(l)->at_boundary())
3006template <
int dim,
int spacedim>
3009 const unsigned int face,
3010 const unsigned int subface)
const
3012 Assert(!this->has_children(),
3013 ExcMessage(
"The present cell must not have children!"));
3014 Assert(!this->at_boundary(face),
3015 ExcMessage(
"The present cell must have a valid neighbor!"));
3016 Assert(this->neighbor(face)->has_children() ==
true,
3017 ExcMessage(
"The neighbor must have children!"));
3025 const unsigned int neighbor_neighbor =
3026 this->neighbor_of_neighbor(face);
3027 const auto neighbor = this->neighbor(face);
3029 Assert(neighbor->refinement_case() ==
3035 const auto relative_orientation =
3036 neighbor->combined_face_orientation(neighbor_neighbor) ==
3037 this->combined_face_orientation(face) ?
3040 const unsigned int neighbor_child_index =
3041 neighbor->reference_cell().child_cell_on_face(
3042 neighbor_neighbor, subface, relative_orientation);
3043 auto child = neighbor->child(neighbor_child_index);
3049 const unsigned int neighbor_neighbor =
3050 this->neighbor_of_neighbor(face);
3051 const unsigned int neighbor_child_index =
3053 this->neighbor(face)->refinement_case(),
3058 this->neighbor(face)->child(neighbor_child_index);
3064 while (sub_neighbor->has_children())
3067 sub_neighbor->refinement_case(),
3068 neighbor_neighbor) ==
3073 sub_neighbor->refinement_case(), neighbor_neighbor, 0));
3076 return sub_neighbor;
3144 mother_face = this->face(face);
3145 const unsigned int total_children =
3146 mother_face->n_active_descendants();
3151 unsigned int neighbor_neighbor;
3154 this->neighbor(face);
3158 mother_face->refinement_case();
3159 if (mother_face_ref_case ==
3168 neighbor_neighbor = this->neighbor_of_neighbor(face);
3172 const unsigned int neighbor_child_index =
3174 neighbor->refinement_case(),
3177 neighbor->face_orientation(neighbor_neighbor),
3178 neighbor->face_flip(neighbor_neighbor),
3179 neighbor->face_rotation(neighbor_neighbor));
3180 neighbor_child = neighbor->child(neighbor_child_index);
3184 Assert((this->face(face)->child(subface) ==
3185 neighbor_child->face(neighbor_neighbor)),
3195 unsigned int first_child_to_find;
3196 unsigned int neighbor_child_index;
3197 if (total_children == 2)
3198 first_child_to_find = subface;
3201 first_child_to_find = subface / 2;
3202 if (total_children == 3 && subface == 1 &&
3203 !mother_face->child(0)->has_children())
3204 first_child_to_find = 1;
3206 if (neighbor_is_coarser(face))
3208 std::pair<unsigned int, unsigned int> indices =
3209 neighbor_of_coarser_neighbor(face);
3210 neighbor_neighbor = indices.first;
3255 unsigned int iso_subface;
3256 if (neighbor->face(neighbor_neighbor)
3258 iso_subface = 2 * first_child_to_find + indices.second;
3261 Assert(neighbor->face(neighbor_neighbor)
3262 ->refinement_case() ==
3266 first_child_to_find + 2 * indices.second;
3268 neighbor_child_index =
3270 neighbor->refinement_case(),
3273 neighbor->face_orientation(neighbor_neighbor),
3274 neighbor->face_flip(neighbor_neighbor),
3275 neighbor->face_rotation(neighbor_neighbor));
3279 neighbor_neighbor = neighbor_of_neighbor(face);
3280 neighbor_child_index =
3282 neighbor->refinement_case(),
3284 first_child_to_find,
3285 neighbor->face_orientation(neighbor_neighbor),
3286 neighbor->face_flip(neighbor_neighbor),
3287 neighbor->face_rotation(neighbor_neighbor),
3288 mother_face_ref_case);
3291 neighbor_child = neighbor->child(neighbor_child_index);
3297 neighbor_child->has_children() &&
3299 neighbor_child->refinement_case(), neighbor_neighbor) ==
3301 neighbor_child = neighbor_child->child(
3303 neighbor_child->refinement_case(),
3313 if (total_children == 3)
3315 if (mother_face->child(0)->has_children())
3318 neighbor_child = neighbor_child->child(
3320 neighbor_child->refinement_case(),
3323 neighbor_child->face_orientation(
3325 neighbor_child->face_flip(neighbor_neighbor),
3326 neighbor_child->face_rotation(
3328 mother_face->child(0)->refinement_case()));
3332 Assert(mother_face->child(1)->has_children(),
3335 neighbor_child = neighbor_child->child(
3337 neighbor_child->refinement_case(),
3340 neighbor_child->face_orientation(
3342 neighbor_child->face_flip(neighbor_neighbor),
3343 neighbor_child->face_rotation(
3345 mother_face->child(1)->refinement_case()));
3348 else if (total_children == 4)
3350 neighbor_child = neighbor_child->child(
3352 neighbor_child->refinement_case(),
3355 neighbor_child->face_orientation(neighbor_neighbor),
3356 neighbor_child->face_flip(neighbor_neighbor),
3357 neighbor_child->face_rotation(neighbor_neighbor),
3358 mother_face->child(subface / 2)->refinement_case()));
3366 while (neighbor_child->has_children())
3369 neighbor_child->refinement_case(), neighbor_neighbor, 0));
3377 switch (this->subface_case(face))
3382 requested = mother_face->child(subface);
3387 mother_face->child(subface / 2)->child(subface % 2);
3396 requested = mother_face->child(0)->child(subface);
3399 requested = mother_face->child(1);
3410 requested = mother_face->child(0);
3415 mother_face->child(1)->child(subface - 1);
3425 Assert(requested == neighbor_child->face(neighbor_neighbor),
3429 return neighbor_child;
3446template <
int structdim,
int dim,
int spacedim>
3455template <
int structdim,
int dim,
int spacedim>
3464template <
int structdim,
int dim,
int spacedim>
3473#include "grid/tria_accessor.inst"