~ubuntu-branches/ubuntu/raring/ffc/raring
» Revision
8
: test/regression/references/QuadratureElement.h
« back to all changes in this revision
Viewing changes to test/regression/references/QuadratureElement.h
- browse files at revision 8
- compare with another revision
- download diff
- download tarball
- view history from revision 8
- Committer: Bazaar Package Importer
- Author(s): Johannes Ring
- Date: 2010-07-01 19:54:32 UTC
- mfrom: (1.1.5 upstream)
- Revision ID: james.westby@ubuntu.com-20100701195432-xz3gw5nprdj79jcb
Tags: 0.9.3-1
* New upstream release.
* debian/control:
- Minor fix in Vcs fields.
- Bump Standards-Version to 3.9.0 (no changes needed).
- Update version for python-ufc, python-fiat, and python-ufl in
Depends field.
* Switch to dpkg-source 3.0 (quilt) format.
* Update debian/copyright and debian/copyright_hints.
* debian/control:
- Minor fix in Vcs fields.
- Bump Standards-Version to 3.9.0 (no changes needed).
- Update version for python-ufc, python-fiat, and python-ufl in
Depends field.
* Switch to dpkg-source 3.0 (quilt) format.
* Update debian/copyright and debian/copyright_hints.
- files added:
- bench/MassH1_2D_1.ufl
- bench/results
- debian/source
- files removed:
- bench/NavierStokes.form
- files modified:
- ChangeLog
added
removed
test/regression/references/QuadratureElement.h
1
// This code conforms with the UFC specification version 1.4
2
// and was automatically generated by FFC version 0.9.2.
1
// This code conforms with the UFC specification version 1.4.1
2
// and was automatically generated by FFC version 0.9.3.
3
3
//
4
4
// This code was generated with the following parameters:
5
5
//
546
546
547
547
// Reset values.
548
548
*values = 0.00000000;
549
550
// Map degree of freedom to element degree of freedom
551
const unsigned int dof = i;
552
553
// Array of basisvalues.
554
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
555
556
// Declare helper variables.
557
unsigned int rr = 0;
558
unsigned int ss = 0;
559
unsigned int tt = 0;
560
double tmp5 = 0.00000000;
561
double tmp6 = 0.00000000;
562
double tmp7 = 0.00000000;
563
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
564
double tmp1 = (1.00000000 - Y)/2.00000000;
565
double tmp2 = tmp1*tmp1;
566
567
// Compute basisvalues.
568
basisvalues[0] = 1.00000000;
569
basisvalues[1] = tmp0;
570
for (unsigned int r = 1; r < 2; r++)
571
{
572
rr = (r + 1)*((r + 1) + 1)/2;
573
ss = r*(r + 1)/2;
574
tt = (r - 1)*((r - 1) + 1)/2;
575
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
576
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
577
}// end loop over 'r'
578
for (unsigned int r = 0; r < 2; r++)
579
{
580
rr = (r + 1)*(r + 1 + 1)/2 + 1;
581
ss = r*(r + 1)/2;
582
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
583
}// end loop over 'r'
584
for (unsigned int r = 0; r < 1; r++)
585
{
586
for (unsigned int s = 1; s < 2 - r; s++)
587
{
588
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
589
ss = (r + s)*(r + s + 1)/2 + s;
590
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
591
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
592
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
593
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
594
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
595
}// end loop over 's'
596
}// end loop over 'r'
597
for (unsigned int r = 0; r < 3; r++)
598
{
599
for (unsigned int s = 0; s < 3 - r; s++)
600
{
601
rr = (r + s)*(r + s + 1)/2 + s;
602
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
603
}// end loop over 's'
604
}// end loop over 'r'
605
606
// Table(s) of coefficients.
607
static const double coefficients0[6][6] = \
608
{{0.00000000, -0.17320508, -0.10000000, 0.12171612, 0.09428090, 0.05443311},
609
{0.00000000, 0.17320508, -0.10000000, 0.12171612, -0.09428090, 0.05443311},
610
{0.00000000, 0.00000000, 0.20000000, 0.00000000, 0.00000000, 0.16329932},
611
{0.47140452, 0.23094011, 0.13333333, 0.00000000, 0.18856181, -0.16329932},
612
{0.47140452, -0.23094011, 0.13333333, 0.00000000, -0.18856181, -0.16329932},
613
{0.47140452, 0.00000000, -0.26666667, -0.24343225, 0.00000000, 0.05443311}};
614
615
// Compute value(s).
616
for (unsigned int r = 0; r < 6; r++)
617
{
618
*values += coefficients0[dof][r]*basisvalues[r];
619
}// end loop over 'r'
549
switch (i)
550
{
551
case 0:
552
{
553
554
// Array of basisvalues.
555
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
556
557
// Declare helper variables.
558
unsigned int rr = 0;
559
unsigned int ss = 0;
560
unsigned int tt = 0;
561
double tmp5 = 0.00000000;
562
double tmp6 = 0.00000000;
563
double tmp7 = 0.00000000;
564
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
565
double tmp1 = (1.00000000 - Y)/2.00000000;
566
double tmp2 = tmp1*tmp1;
567
568
// Compute basisvalues.
569
basisvalues[0] = 1.00000000;
570
basisvalues[1] = tmp0;
571
for (unsigned int r = 1; r < 2; r++)
572
{
573
rr = (r + 1)*((r + 1) + 1)/2;
574
ss = r*(r + 1)/2;
575
tt = (r - 1)*((r - 1) + 1)/2;
576
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
577
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
578
}// end loop over 'r'
579
for (unsigned int r = 0; r < 2; r++)
580
{
581
rr = (r + 1)*(r + 1 + 1)/2 + 1;
582
ss = r*(r + 1)/2;
583
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
584
}// end loop over 'r'
585
for (unsigned int r = 0; r < 1; r++)
586
{
587
for (unsigned int s = 1; s < 2 - r; s++)
588
{
589
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
590
ss = (r + s)*(r + s + 1)/2 + s;
591
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
592
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
593
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
594
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
595
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
596
}// end loop over 's'
597
}// end loop over 'r'
598
for (unsigned int r = 0; r < 3; r++)
599
{
600
for (unsigned int s = 0; s < 3 - r; s++)
601
{
602
rr = (r + s)*(r + s + 1)/2 + s;
603
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
604
}// end loop over 's'
605
}// end loop over 'r'
606
607
// Table(s) of coefficients.
608
static const double coefficients0[6] = \
609
{0.00000000, -0.17320508, -0.10000000, 0.12171612, 0.09428090, 0.05443311};
610
611
// Compute value(s).
612
for (unsigned int r = 0; r < 6; r++)
613
{
614
*values += coefficients0[r]*basisvalues[r];
615
}// end loop over 'r'
616
break;
617
}
618
case 1:
619
{
620
621
// Array of basisvalues.
622
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
623
624
// Declare helper variables.
625
unsigned int rr = 0;
626
unsigned int ss = 0;
627
unsigned int tt = 0;
628
double tmp5 = 0.00000000;
629
double tmp6 = 0.00000000;
630
double tmp7 = 0.00000000;
631
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
632
double tmp1 = (1.00000000 - Y)/2.00000000;
633
double tmp2 = tmp1*tmp1;
634
635
// Compute basisvalues.
636
basisvalues[0] = 1.00000000;
637
basisvalues[1] = tmp0;
638
for (unsigned int r = 1; r < 2; r++)
639
{
640
rr = (r + 1)*((r + 1) + 1)/2;
641
ss = r*(r + 1)/2;
642
tt = (r - 1)*((r - 1) + 1)/2;
643
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
644
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
645
}// end loop over 'r'
646
for (unsigned int r = 0; r < 2; r++)
647
{
648
rr = (r + 1)*(r + 1 + 1)/2 + 1;
649
ss = r*(r + 1)/2;
650
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
651
}// end loop over 'r'
652
for (unsigned int r = 0; r < 1; r++)
653
{
654
for (unsigned int s = 1; s < 2 - r; s++)
655
{
656
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
657
ss = (r + s)*(r + s + 1)/2 + s;
658
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
659
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
660
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
661
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
662
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
663
}// end loop over 's'
664
}// end loop over 'r'
665
for (unsigned int r = 0; r < 3; r++)
666
{
667
for (unsigned int s = 0; s < 3 - r; s++)
668
{
669
rr = (r + s)*(r + s + 1)/2 + s;
670
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
671
}// end loop over 's'
672
}// end loop over 'r'
673
674
// Table(s) of coefficients.
675
static const double coefficients0[6] = \
676
{0.00000000, 0.17320508, -0.10000000, 0.12171612, -0.09428090, 0.05443311};
677
678
// Compute value(s).
679
for (unsigned int r = 0; r < 6; r++)
680
{
681
*values += coefficients0[r]*basisvalues[r];
682
}// end loop over 'r'
683
break;
684
}
685
case 2:
686
{
687
688
// Array of basisvalues.
689
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
690
691
// Declare helper variables.
692
unsigned int rr = 0;
693
unsigned int ss = 0;
694
unsigned int tt = 0;
695
double tmp5 = 0.00000000;
696
double tmp6 = 0.00000000;
697
double tmp7 = 0.00000000;
698
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
699
double tmp1 = (1.00000000 - Y)/2.00000000;
700
double tmp2 = tmp1*tmp1;
701
702
// Compute basisvalues.
703
basisvalues[0] = 1.00000000;
704
basisvalues[1] = tmp0;
705
for (unsigned int r = 1; r < 2; r++)
706
{
707
rr = (r + 1)*((r + 1) + 1)/2;
708
ss = r*(r + 1)/2;
709
tt = (r - 1)*((r - 1) + 1)/2;
710
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
711
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
712
}// end loop over 'r'
713
for (unsigned int r = 0; r < 2; r++)
714
{
715
rr = (r + 1)*(r + 1 + 1)/2 + 1;
716
ss = r*(r + 1)/2;
717
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
718
}// end loop over 'r'
719
for (unsigned int r = 0; r < 1; r++)
720
{
721
for (unsigned int s = 1; s < 2 - r; s++)
722
{
723
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
724
ss = (r + s)*(r + s + 1)/2 + s;
725
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
726
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
727
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
728
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
729
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
730
}// end loop over 's'
731
}// end loop over 'r'
732
for (unsigned int r = 0; r < 3; r++)
733
{
734
for (unsigned int s = 0; s < 3 - r; s++)
735
{
736
rr = (r + s)*(r + s + 1)/2 + s;
737
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
738
}// end loop over 's'
739
}// end loop over 'r'
740
741
// Table(s) of coefficients.
742
static const double coefficients0[6] = \
743
{0.00000000, 0.00000000, 0.20000000, 0.00000000, 0.00000000, 0.16329932};
744
745
// Compute value(s).
746
for (unsigned int r = 0; r < 6; r++)
747
{
748
*values += coefficients0[r]*basisvalues[r];
749
}// end loop over 'r'
750
break;
751
}
752
case 3:
753
{
754
755
// Array of basisvalues.
756
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
757
758
// Declare helper variables.
759
unsigned int rr = 0;
760
unsigned int ss = 0;
761
unsigned int tt = 0;
762
double tmp5 = 0.00000000;
763
double tmp6 = 0.00000000;
764
double tmp7 = 0.00000000;
765
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
766
double tmp1 = (1.00000000 - Y)/2.00000000;
767
double tmp2 = tmp1*tmp1;
768
769
// Compute basisvalues.
770
basisvalues[0] = 1.00000000;
771
basisvalues[1] = tmp0;
772
for (unsigned int r = 1; r < 2; r++)
773
{
774
rr = (r + 1)*((r + 1) + 1)/2;
775
ss = r*(r + 1)/2;
776
tt = (r - 1)*((r - 1) + 1)/2;
777
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
778
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
779
}// end loop over 'r'
780
for (unsigned int r = 0; r < 2; r++)
781
{
782
rr = (r + 1)*(r + 1 + 1)/2 + 1;
783
ss = r*(r + 1)/2;
784
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
785
}// end loop over 'r'
786
for (unsigned int r = 0; r < 1; r++)
787
{
788
for (unsigned int s = 1; s < 2 - r; s++)
789
{
790
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
791
ss = (r + s)*(r + s + 1)/2 + s;
792
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
793
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
794
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
795
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
796
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
797
}// end loop over 's'
798
}// end loop over 'r'
799
for (unsigned int r = 0; r < 3; r++)
800
{
801
for (unsigned int s = 0; s < 3 - r; s++)
802
{
803
rr = (r + s)*(r + s + 1)/2 + s;
804
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
805
}// end loop over 's'
806
}// end loop over 'r'
807
808
// Table(s) of coefficients.
809
static const double coefficients0[6] = \
810
{0.47140452, 0.23094011, 0.13333333, 0.00000000, 0.18856181, -0.16329932};
811
812
// Compute value(s).
813
for (unsigned int r = 0; r < 6; r++)
814
{
815
*values += coefficients0[r]*basisvalues[r];
816
}// end loop over 'r'
817
break;
818
}
819
case 4:
820
{
821
822
// Array of basisvalues.
823
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
824
825
// Declare helper variables.
826
unsigned int rr = 0;
827
unsigned int ss = 0;
828
unsigned int tt = 0;
829
double tmp5 = 0.00000000;
830
double tmp6 = 0.00000000;
831
double tmp7 = 0.00000000;
832
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
833
double tmp1 = (1.00000000 - Y)/2.00000000;
834
double tmp2 = tmp1*tmp1;
835
836
// Compute basisvalues.
837
basisvalues[0] = 1.00000000;
838
basisvalues[1] = tmp0;
839
for (unsigned int r = 1; r < 2; r++)
840
{
841
rr = (r + 1)*((r + 1) + 1)/2;
842
ss = r*(r + 1)/2;
843
tt = (r - 1)*((r - 1) + 1)/2;
844
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
845
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
846
}// end loop over 'r'
847
for (unsigned int r = 0; r < 2; r++)
848
{
849
rr = (r + 1)*(r + 1 + 1)/2 + 1;
850
ss = r*(r + 1)/2;
851
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
852
}// end loop over 'r'
853
for (unsigned int r = 0; r < 1; r++)
854
{
855
for (unsigned int s = 1; s < 2 - r; s++)
856
{
857
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
858
ss = (r + s)*(r + s + 1)/2 + s;
859
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
860
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
861
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
862
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
863
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
864
}// end loop over 's'
865
}// end loop over 'r'
866
for (unsigned int r = 0; r < 3; r++)
867
{
868
for (unsigned int s = 0; s < 3 - r; s++)
869
{
870
rr = (r + s)*(r + s + 1)/2 + s;
871
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
872
}// end loop over 's'
873
}// end loop over 'r'
874
875
// Table(s) of coefficients.
876
static const double coefficients0[6] = \
877
{0.47140452, -0.23094011, 0.13333333, 0.00000000, -0.18856181, -0.16329932};
878
879
// Compute value(s).
880
for (unsigned int r = 0; r < 6; r++)
881
{
882
*values += coefficients0[r]*basisvalues[r];
883
}// end loop over 'r'
884
break;
885
}
886
case 5:
887
{
888
889
// Array of basisvalues.
890
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
891
892
// Declare helper variables.
893
unsigned int rr = 0;
894
unsigned int ss = 0;
895
unsigned int tt = 0;
896
double tmp5 = 0.00000000;
897
double tmp6 = 0.00000000;
898
double tmp7 = 0.00000000;
899
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
900
double tmp1 = (1.00000000 - Y)/2.00000000;
901
double tmp2 = tmp1*tmp1;
902
903
// Compute basisvalues.
904
basisvalues[0] = 1.00000000;
905
basisvalues[1] = tmp0;
906
for (unsigned int r = 1; r < 2; r++)
907
{
908
rr = (r + 1)*((r + 1) + 1)/2;
909
ss = r*(r + 1)/2;
910
tt = (r - 1)*((r - 1) + 1)/2;
911
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
912
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
913
}// end loop over 'r'
914
for (unsigned int r = 0; r < 2; r++)
915
{
916
rr = (r + 1)*(r + 1 + 1)/2 + 1;
917
ss = r*(r + 1)/2;
918
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
919
}// end loop over 'r'
920
for (unsigned int r = 0; r < 1; r++)
921
{
922
for (unsigned int s = 1; s < 2 - r; s++)
923
{
924
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
925
ss = (r + s)*(r + s + 1)/2 + s;
926
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
927
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
928
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
929
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
930
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
931
}// end loop over 's'
932
}// end loop over 'r'
933
for (unsigned int r = 0; r < 3; r++)
934
{
935
for (unsigned int s = 0; s < 3 - r; s++)
936
{
937
rr = (r + s)*(r + s + 1)/2 + s;
938
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
939
}// end loop over 's'
940
}// end loop over 'r'
941
942
// Table(s) of coefficients.
943
static const double coefficients0[6] = \
944
{0.47140452, 0.00000000, -0.26666667, -0.24343225, 0.00000000, 0.05443311};
945
946
// Compute value(s).
947
for (unsigned int r = 0; r < 6; r++)
948
{
949
*values += coefficients0[r]*basisvalues[r];
950
}// end loop over 'r'
951
break;
952
}
953
}
954
620
955
}
621
956
622
957
/// Evaluate all basis functions at given point in cell
732
1067
values[r] = 0.00000000;
733
1068
}// end loop over 'r'
734
1069
735
// Map degree of freedom to element degree of freedom
736
const unsigned int dof = i;
737
738
// Array of basisvalues.
739
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
740
741
// Declare helper variables.
742
unsigned int rr = 0;
743
unsigned int ss = 0;
744
unsigned int tt = 0;
745
double tmp5 = 0.00000000;
746
double tmp6 = 0.00000000;
747
double tmp7 = 0.00000000;
748
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
749
double tmp1 = (1.00000000 - Y)/2.00000000;
750
double tmp2 = tmp1*tmp1;
751
752
// Compute basisvalues.
753
basisvalues[0] = 1.00000000;
754
basisvalues[1] = tmp0;
755
for (unsigned int r = 1; r < 2; r++)
756
{
757
rr = (r + 1)*((r + 1) + 1)/2;
758
ss = r*(r + 1)/2;
759
tt = (r - 1)*((r - 1) + 1)/2;
760
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
761
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
762
}// end loop over 'r'
763
for (unsigned int r = 0; r < 2; r++)
764
{
765
rr = (r + 1)*(r + 1 + 1)/2 + 1;
766
ss = r*(r + 1)/2;
767
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
768
}// end loop over 'r'
769
for (unsigned int r = 0; r < 1; r++)
770
{
771
for (unsigned int s = 1; s < 2 - r; s++)
772
{
773
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
774
ss = (r + s)*(r + s + 1)/2 + s;
775
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
776
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
777
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
778
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
779
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
780
}// end loop over 's'
781
}// end loop over 'r'
782
for (unsigned int r = 0; r < 3; r++)
783
{
784
for (unsigned int s = 0; s < 3 - r; s++)
785
{
786
rr = (r + s)*(r + s + 1)/2 + s;
787
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
788
}// end loop over 's'
789
}// end loop over 'r'
790
791
// Table(s) of coefficients.
792
static const double coefficients0[6][6] = \
793
{{0.00000000, -0.17320508, -0.10000000, 0.12171612, 0.09428090, 0.05443311},
794
{0.00000000, 0.17320508, -0.10000000, 0.12171612, -0.09428090, 0.05443311},
795
{0.00000000, 0.00000000, 0.20000000, 0.00000000, 0.00000000, 0.16329932},
796
{0.47140452, 0.23094011, 0.13333333, 0.00000000, 0.18856181, -0.16329932},
797
{0.47140452, -0.23094011, 0.13333333, 0.00000000, -0.18856181, -0.16329932},
798
{0.47140452, 0.00000000, -0.26666667, -0.24343225, 0.00000000, 0.05443311}};
799
800
// Tables of derivatives of the polynomial base (transpose).
801
static const double dmats0[6][6] = \
802
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
803
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
804
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
805
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
806
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000},
807
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
808
809
static const double dmats1[6][6] = \
810
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
811
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
812
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
813
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000},
814
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000},
815
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000}};
816
817
// Compute reference derivatives.
818
// Declare pointer to array of derivatives on FIAT element.
819
double *derivatives = new double[num_derivatives];
820
for (unsigned int r = 0; r < num_derivatives; r++)
821
{
822
derivatives[r] = 0.00000000;
823
}// end loop over 'r'
824
825
// Declare derivative matrix (of polynomial basis).
826
double dmats[6][6] = \
827
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
828
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
829
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
830
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
831
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
832
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
833
834
// Declare (auxiliary) derivative matrix (of polynomial basis).
835
double dmats_old[6][6] = \
836
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
837
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
838
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
839
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
840
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
841
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
842
843
// Loop possible derivatives.
844
for (unsigned int r = 0; r < num_derivatives; r++)
845
{
846
// Resetting dmats values to compute next derivative.
847
for (unsigned int t = 0; t < 6; t++)
848
{
849
for (unsigned int u = 0; u < 6; u++)
850
{
851
dmats[t][u] = 0.00000000;
852
if (t == u)
853
{
854
dmats[t][u] = 1.00000000;
855
}
856
857
}// end loop over 'u'
858
}// end loop over 't'
859
860
// Looping derivative order to generate dmats.
861
for (unsigned int s = 0; s < n; s++)
862
{
863
// Updating dmats_old with new values and resetting dmats.
864
for (unsigned int t = 0; t < 6; t++)
865
{
866
for (unsigned int u = 0; u < 6; u++)
867
{
868
dmats_old[t][u] = dmats[t][u];
869
dmats[t][u] = 0.00000000;
870
}// end loop over 'u'
871
}// end loop over 't'
872
873
// Update dmats using an inner product.
874
if (combinations[r][s] == 0)
875
{
876
for (unsigned int t = 0; t < 6; t++)
877
{
878
for (unsigned int u = 0; u < 6; u++)
879
{
880
for (unsigned int tu = 0; tu < 6; tu++)
881
{
882
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
883
}// end loop over 'tu'
884
}// end loop over 'u'
885
}// end loop over 't'
886
}
887
888
if (combinations[r][s] == 1)
889
{
890
for (unsigned int t = 0; t < 6; t++)
891
{
892
for (unsigned int u = 0; u < 6; u++)
893
{
894
for (unsigned int tu = 0; tu < 6; tu++)
895
{
896
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
897
}// end loop over 'tu'
898
}// end loop over 'u'
899
}// end loop over 't'
900
}
901
902
}// end loop over 's'
903
for (unsigned int s = 0; s < 6; s++)
904
{
905
for (unsigned int t = 0; t < 6; t++)
906
{
907
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
908
}// end loop over 't'
909
}// end loop over 's'
910
}// end loop over 'r'
911
912
// Transform derivatives back to physical element
913
for (unsigned int r = 0; r < num_derivatives; r++)
914
{
915
for (unsigned int s = 0; s < num_derivatives; s++)
916
{
917
values[r] += transform[r][s]*derivatives[s];
918
}// end loop over 's'
919
}// end loop over 'r'
920
921
// Delete pointer to array of derivatives on FIAT element
922
delete [] derivatives;
923
924
// Delete pointer to array of combinations of derivatives and transform
925
for (unsigned int r = 0; r < num_derivatives; r++)
926
{
927
delete [] combinations[r];
928
}// end loop over 'r'
929
delete [] combinations;
930
for (unsigned int r = 0; r < num_derivatives; r++)
931
{
932
delete [] transform[r];
933
}// end loop over 'r'
934
delete [] transform;
1070
switch (i)
1071
{
1072
case 0:
1073
{
1074
1075
// Array of basisvalues.
1076
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
1077
1078
// Declare helper variables.
1079
unsigned int rr = 0;
1080
unsigned int ss = 0;
1081
unsigned int tt = 0;
1082
double tmp5 = 0.00000000;
1083
double tmp6 = 0.00000000;
1084
double tmp7 = 0.00000000;
1085
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1086
double tmp1 = (1.00000000 - Y)/2.00000000;
1087
double tmp2 = tmp1*tmp1;
1088
1089
// Compute basisvalues.
1090
basisvalues[0] = 1.00000000;
1091
basisvalues[1] = tmp0;
1092
for (unsigned int r = 1; r < 2; r++)
1093
{
1094
rr = (r + 1)*((r + 1) + 1)/2;
1095
ss = r*(r + 1)/2;
1096
tt = (r - 1)*((r - 1) + 1)/2;
1097
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
1098
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
1099
}// end loop over 'r'
1100
for (unsigned int r = 0; r < 2; r++)
1101
{
1102
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1103
ss = r*(r + 1)/2;
1104
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1105
}// end loop over 'r'
1106
for (unsigned int r = 0; r < 1; r++)
1107
{
1108
for (unsigned int s = 1; s < 2 - r; s++)
1109
{
1110
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
1111
ss = (r + s)*(r + s + 1)/2 + s;
1112
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
1113
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1114
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1115
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1116
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
1117
}// end loop over 's'
1118
}// end loop over 'r'
1119
for (unsigned int r = 0; r < 3; r++)
1120
{
1121
for (unsigned int s = 0; s < 3 - r; s++)
1122
{
1123
rr = (r + s)*(r + s + 1)/2 + s;
1124
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1125
}// end loop over 's'
1126
}// end loop over 'r'
1127
1128
// Table(s) of coefficients.
1129
static const double coefficients0[6] = \
1130
{0.00000000, -0.17320508, -0.10000000, 0.12171612, 0.09428090, 0.05443311};
1131
1132
// Tables of derivatives of the polynomial base (transpose).
1133
static const double dmats0[6][6] = \
1134
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1135
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1136
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1137
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1138
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000},
1139
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1140
1141
static const double dmats1[6][6] = \
1142
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1143
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1144
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1145
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000},
1146
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000},
1147
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000}};
1148
1149
// Compute reference derivatives.
1150
// Declare pointer to array of derivatives on FIAT element.
1151
double *derivatives = new double[num_derivatives];
1152
for (unsigned int r = 0; r < num_derivatives; r++)
1153
{
1154
derivatives[r] = 0.00000000;
1155
}// end loop over 'r'
1156
1157
// Declare derivative matrix (of polynomial basis).
1158
double dmats[6][6] = \
1159
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1160
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1161
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1162
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1163
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1164
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1165
1166
// Declare (auxiliary) derivative matrix (of polynomial basis).
1167
double dmats_old[6][6] = \
1168
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1169
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1170
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1171
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1172
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1173
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1174
1175
// Loop possible derivatives.
1176
for (unsigned int r = 0; r < num_derivatives; r++)
1177
{
1178
// Resetting dmats values to compute next derivative.
1179
for (unsigned int t = 0; t < 6; t++)
1180
{
1181
for (unsigned int u = 0; u < 6; u++)
1182
{
1183
dmats[t][u] = 0.00000000;
1184
if (t == u)
1185
{
1186
dmats[t][u] = 1.00000000;
1187
}
1188
1189
}// end loop over 'u'
1190
}// end loop over 't'
1191
1192
// Looping derivative order to generate dmats.
1193
for (unsigned int s = 0; s < n; s++)
1194
{
1195
// Updating dmats_old with new values and resetting dmats.
1196
for (unsigned int t = 0; t < 6; t++)
1197
{
1198
for (unsigned int u = 0; u < 6; u++)
1199
{
1200
dmats_old[t][u] = dmats[t][u];
1201
dmats[t][u] = 0.00000000;
1202
}// end loop over 'u'
1203
}// end loop over 't'
1204
1205
// Update dmats using an inner product.
1206
if (combinations[r][s] == 0)
1207
{
1208
for (unsigned int t = 0; t < 6; t++)
1209
{
1210
for (unsigned int u = 0; u < 6; u++)
1211
{
1212
for (unsigned int tu = 0; tu < 6; tu++)
1213
{
1214
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1215
}// end loop over 'tu'
1216
}// end loop over 'u'
1217
}// end loop over 't'
1218
}
1219
1220
if (combinations[r][s] == 1)
1221
{
1222
for (unsigned int t = 0; t < 6; t++)
1223
{
1224
for (unsigned int u = 0; u < 6; u++)
1225
{
1226
for (unsigned int tu = 0; tu < 6; tu++)
1227
{
1228
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1229
}// end loop over 'tu'
1230
}// end loop over 'u'
1231
}// end loop over 't'
1232
}
1233
1234
}// end loop over 's'
1235
for (unsigned int s = 0; s < 6; s++)
1236
{
1237
for (unsigned int t = 0; t < 6; t++)
1238
{
1239
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1240
}// end loop over 't'
1241
}// end loop over 's'
1242
}// end loop over 'r'
1243
1244
// Transform derivatives back to physical element
1245
for (unsigned int r = 0; r < num_derivatives; r++)
1246
{
1247
for (unsigned int s = 0; s < num_derivatives; s++)
1248
{
1249
values[r] += transform[r][s]*derivatives[s];
1250
}// end loop over 's'
1251
}// end loop over 'r'
1252
1253
// Delete pointer to array of derivatives on FIAT element
1254
delete [] derivatives;
1255
1256
// Delete pointer to array of combinations of derivatives and transform
1257
for (unsigned int r = 0; r < num_derivatives; r++)
1258
{
1259
delete [] combinations[r];
1260
}// end loop over 'r'
1261
delete [] combinations;
1262
for (unsigned int r = 0; r < num_derivatives; r++)
1263
{
1264
delete [] transform[r];
1265
}// end loop over 'r'
1266
delete [] transform;
1267
break;
1268
}
1269
case 1:
1270
{
1271
1272
// Array of basisvalues.
1273
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
1274
1275
// Declare helper variables.
1276
unsigned int rr = 0;
1277
unsigned int ss = 0;
1278
unsigned int tt = 0;
1279
double tmp5 = 0.00000000;
1280
double tmp6 = 0.00000000;
1281
double tmp7 = 0.00000000;
1282
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1283
double tmp1 = (1.00000000 - Y)/2.00000000;
1284
double tmp2 = tmp1*tmp1;
1285
1286
// Compute basisvalues.
1287
basisvalues[0] = 1.00000000;
1288
basisvalues[1] = tmp0;
1289
for (unsigned int r = 1; r < 2; r++)
1290
{
1291
rr = (r + 1)*((r + 1) + 1)/2;
1292
ss = r*(r + 1)/2;
1293
tt = (r - 1)*((r - 1) + 1)/2;
1294
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
1295
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
1296
}// end loop over 'r'
1297
for (unsigned int r = 0; r < 2; r++)
1298
{
1299
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1300
ss = r*(r + 1)/2;
1301
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1302
}// end loop over 'r'
1303
for (unsigned int r = 0; r < 1; r++)
1304
{
1305
for (unsigned int s = 1; s < 2 - r; s++)
1306
{
1307
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
1308
ss = (r + s)*(r + s + 1)/2 + s;
1309
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
1310
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1311
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1312
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1313
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
1314
}// end loop over 's'
1315
}// end loop over 'r'
1316
for (unsigned int r = 0; r < 3; r++)
1317
{
1318
for (unsigned int s = 0; s < 3 - r; s++)
1319
{
1320
rr = (r + s)*(r + s + 1)/2 + s;
1321
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1322
}// end loop over 's'
1323
}// end loop over 'r'
1324
1325
// Table(s) of coefficients.
1326
static const double coefficients0[6] = \
1327
{0.00000000, 0.17320508, -0.10000000, 0.12171612, -0.09428090, 0.05443311};
1328
1329
// Tables of derivatives of the polynomial base (transpose).
1330
static const double dmats0[6][6] = \
1331
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1332
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1333
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1334
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1335
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000},
1336
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1337
1338
static const double dmats1[6][6] = \
1339
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1340
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1341
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1342
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000},
1343
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000},
1344
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000}};
1345
1346
// Compute reference derivatives.
1347
// Declare pointer to array of derivatives on FIAT element.
1348
double *derivatives = new double[num_derivatives];
1349
for (unsigned int r = 0; r < num_derivatives; r++)
1350
{
1351
derivatives[r] = 0.00000000;
1352
}// end loop over 'r'
1353
1354
// Declare derivative matrix (of polynomial basis).
1355
double dmats[6][6] = \
1356
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1357
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1358
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1359
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1360
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1361
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1362
1363
// Declare (auxiliary) derivative matrix (of polynomial basis).
1364
double dmats_old[6][6] = \
1365
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1366
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1367
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1368
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1369
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1370
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1371
1372
// Loop possible derivatives.
1373
for (unsigned int r = 0; r < num_derivatives; r++)
1374
{
1375
// Resetting dmats values to compute next derivative.
1376
for (unsigned int t = 0; t < 6; t++)
1377
{
1378
for (unsigned int u = 0; u < 6; u++)
1379
{
1380
dmats[t][u] = 0.00000000;
1381
if (t == u)
1382
{
1383
dmats[t][u] = 1.00000000;
1384
}
1385
1386
}// end loop over 'u'
1387
}// end loop over 't'
1388
1389
// Looping derivative order to generate dmats.
1390
for (unsigned int s = 0; s < n; s++)
1391
{
1392
// Updating dmats_old with new values and resetting dmats.
1393
for (unsigned int t = 0; t < 6; t++)
1394
{
1395
for (unsigned int u = 0; u < 6; u++)
1396
{
1397
dmats_old[t][u] = dmats[t][u];
1398
dmats[t][u] = 0.00000000;
1399
}// end loop over 'u'
1400
}// end loop over 't'
1401
1402
// Update dmats using an inner product.
1403
if (combinations[r][s] == 0)
1404
{
1405
for (unsigned int t = 0; t < 6; t++)
1406
{
1407
for (unsigned int u = 0; u < 6; u++)
1408
{
1409
for (unsigned int tu = 0; tu < 6; tu++)
1410
{
1411
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1412
}// end loop over 'tu'
1413
}// end loop over 'u'
1414
}// end loop over 't'
1415
}
1416
1417
if (combinations[r][s] == 1)
1418
{
1419
for (unsigned int t = 0; t < 6; t++)
1420
{
1421
for (unsigned int u = 0; u < 6; u++)
1422
{
1423
for (unsigned int tu = 0; tu < 6; tu++)
1424
{
1425
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1426
}// end loop over 'tu'
1427
}// end loop over 'u'
1428
}// end loop over 't'
1429
}
1430
1431
}// end loop over 's'
1432
for (unsigned int s = 0; s < 6; s++)
1433
{
1434
for (unsigned int t = 0; t < 6; t++)
1435
{
1436
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1437
}// end loop over 't'
1438
}// end loop over 's'
1439
}// end loop over 'r'
1440
1441
// Transform derivatives back to physical element
1442
for (unsigned int r = 0; r < num_derivatives; r++)
1443
{
1444
for (unsigned int s = 0; s < num_derivatives; s++)
1445
{
1446
values[r] += transform[r][s]*derivatives[s];
1447
}// end loop over 's'
1448
}// end loop over 'r'
1449
1450
// Delete pointer to array of derivatives on FIAT element
1451
delete [] derivatives;
1452
1453
// Delete pointer to array of combinations of derivatives and transform
1454
for (unsigned int r = 0; r < num_derivatives; r++)
1455
{
1456
delete [] combinations[r];
1457
}// end loop over 'r'
1458
delete [] combinations;
1459
for (unsigned int r = 0; r < num_derivatives; r++)
1460
{
1461
delete [] transform[r];
1462
}// end loop over 'r'
1463
delete [] transform;
1464
break;
1465
}
1466
case 2:
1467
{
1468
1469
// Array of basisvalues.
1470
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
1471
1472
// Declare helper variables.
1473
unsigned int rr = 0;
1474
unsigned int ss = 0;
1475
unsigned int tt = 0;
1476
double tmp5 = 0.00000000;
1477
double tmp6 = 0.00000000;
1478
double tmp7 = 0.00000000;
1479
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1480
double tmp1 = (1.00000000 - Y)/2.00000000;
1481
double tmp2 = tmp1*tmp1;
1482
1483
// Compute basisvalues.
1484
basisvalues[0] = 1.00000000;
1485
basisvalues[1] = tmp0;
1486
for (unsigned int r = 1; r < 2; r++)
1487
{
1488
rr = (r + 1)*((r + 1) + 1)/2;
1489
ss = r*(r + 1)/2;
1490
tt = (r - 1)*((r - 1) + 1)/2;
1491
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
1492
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
1493
}// end loop over 'r'
1494
for (unsigned int r = 0; r < 2; r++)
1495
{
1496
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1497
ss = r*(r + 1)/2;
1498
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1499
}// end loop over 'r'
1500
for (unsigned int r = 0; r < 1; r++)
1501
{
1502
for (unsigned int s = 1; s < 2 - r; s++)
1503
{
1504
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
1505
ss = (r + s)*(r + s + 1)/2 + s;
1506
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
1507
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1508
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1509
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1510
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
1511
}// end loop over 's'
1512
}// end loop over 'r'
1513
for (unsigned int r = 0; r < 3; r++)
1514
{
1515
for (unsigned int s = 0; s < 3 - r; s++)
1516
{
1517
rr = (r + s)*(r + s + 1)/2 + s;
1518
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1519
}// end loop over 's'
1520
}// end loop over 'r'
1521
1522
// Table(s) of coefficients.
1523
static const double coefficients0[6] = \
1524
{0.00000000, 0.00000000, 0.20000000, 0.00000000, 0.00000000, 0.16329932};
1525
1526
// Tables of derivatives of the polynomial base (transpose).
1527
static const double dmats0[6][6] = \
1528
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1529
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1530
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1531
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1532
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000},
1533
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1534
1535
static const double dmats1[6][6] = \
1536
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1537
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1538
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1539
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000},
1540
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000},
1541
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000}};
1542
1543
// Compute reference derivatives.
1544
// Declare pointer to array of derivatives on FIAT element.
1545
double *derivatives = new double[num_derivatives];
1546
for (unsigned int r = 0; r < num_derivatives; r++)
1547
{
1548
derivatives[r] = 0.00000000;
1549
}// end loop over 'r'
1550
1551
// Declare derivative matrix (of polynomial basis).
1552
double dmats[6][6] = \
1553
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1554
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1555
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1556
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1557
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1558
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1559
1560
// Declare (auxiliary) derivative matrix (of polynomial basis).
1561
double dmats_old[6][6] = \
1562
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1563
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1564
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1565
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1566
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1567
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1568
1569
// Loop possible derivatives.
1570
for (unsigned int r = 0; r < num_derivatives; r++)
1571
{
1572
// Resetting dmats values to compute next derivative.
1573
for (unsigned int t = 0; t < 6; t++)
1574
{
1575
for (unsigned int u = 0; u < 6; u++)
1576
{
1577
dmats[t][u] = 0.00000000;
1578
if (t == u)
1579
{
1580
dmats[t][u] = 1.00000000;
1581
}
1582
1583
}// end loop over 'u'
1584
}// end loop over 't'
1585
1586
// Looping derivative order to generate dmats.
1587
for (unsigned int s = 0; s < n; s++)
1588
{
1589
// Updating dmats_old with new values and resetting dmats.
1590
for (unsigned int t = 0; t < 6; t++)
1591
{
1592
for (unsigned int u = 0; u < 6; u++)
1593
{
1594
dmats_old[t][u] = dmats[t][u];
1595
dmats[t][u] = 0.00000000;
1596
}// end loop over 'u'
1597
}// end loop over 't'
1598
1599
// Update dmats using an inner product.
1600
if (combinations[r][s] == 0)
1601
{
1602
for (unsigned int t = 0; t < 6; t++)
1603
{
1604
for (unsigned int u = 0; u < 6; u++)
1605
{
1606
for (unsigned int tu = 0; tu < 6; tu++)
1607
{
1608
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1609
}// end loop over 'tu'
1610
}// end loop over 'u'
1611
}// end loop over 't'
1612
}
1613
1614
if (combinations[r][s] == 1)
1615
{
1616
for (unsigned int t = 0; t < 6; t++)
1617
{
1618
for (unsigned int u = 0; u < 6; u++)
1619
{
1620
for (unsigned int tu = 0; tu < 6; tu++)
1621
{
1622
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1623
}// end loop over 'tu'
1624
}// end loop over 'u'
1625
}// end loop over 't'
1626
}
1627
1628
}// end loop over 's'
1629
for (unsigned int s = 0; s < 6; s++)
1630
{
1631
for (unsigned int t = 0; t < 6; t++)
1632
{
1633
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1634
}// end loop over 't'
1635
}// end loop over 's'
1636
}// end loop over 'r'
1637
1638
// Transform derivatives back to physical element
1639
for (unsigned int r = 0; r < num_derivatives; r++)
1640
{
1641
for (unsigned int s = 0; s < num_derivatives; s++)
1642
{
1643
values[r] += transform[r][s]*derivatives[s];
1644
}// end loop over 's'
1645
}// end loop over 'r'
1646
1647
// Delete pointer to array of derivatives on FIAT element
1648
delete [] derivatives;
1649
1650
// Delete pointer to array of combinations of derivatives and transform
1651
for (unsigned int r = 0; r < num_derivatives; r++)
1652
{
1653
delete [] combinations[r];
1654
}// end loop over 'r'
1655
delete [] combinations;
1656
for (unsigned int r = 0; r < num_derivatives; r++)
1657
{
1658
delete [] transform[r];
1659
}// end loop over 'r'
1660
delete [] transform;
1661
break;
1662
}
1663
case 3:
1664
{
1665
1666
// Array of basisvalues.
1667
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
1668
1669
// Declare helper variables.
1670
unsigned int rr = 0;
1671
unsigned int ss = 0;
1672
unsigned int tt = 0;
1673
double tmp5 = 0.00000000;
1674
double tmp6 = 0.00000000;
1675
double tmp7 = 0.00000000;
1676
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1677
double tmp1 = (1.00000000 - Y)/2.00000000;
1678
double tmp2 = tmp1*tmp1;
1679
1680
// Compute basisvalues.
1681
basisvalues[0] = 1.00000000;
1682
basisvalues[1] = tmp0;
1683
for (unsigned int r = 1; r < 2; r++)
1684
{
1685
rr = (r + 1)*((r + 1) + 1)/2;
1686
ss = r*(r + 1)/2;
1687
tt = (r - 1)*((r - 1) + 1)/2;
1688
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
1689
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
1690
}// end loop over 'r'
1691
for (unsigned int r = 0; r < 2; r++)
1692
{
1693
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1694
ss = r*(r + 1)/2;
1695
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1696
}// end loop over 'r'
1697
for (unsigned int r = 0; r < 1; r++)
1698
{
1699
for (unsigned int s = 1; s < 2 - r; s++)
1700
{
1701
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
1702
ss = (r + s)*(r + s + 1)/2 + s;
1703
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
1704
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1705
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1706
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1707
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
1708
}// end loop over 's'
1709
}// end loop over 'r'
1710
for (unsigned int r = 0; r < 3; r++)
1711
{
1712
for (unsigned int s = 0; s < 3 - r; s++)
1713
{
1714
rr = (r + s)*(r + s + 1)/2 + s;
1715
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1716
}// end loop over 's'
1717
}// end loop over 'r'
1718
1719
// Table(s) of coefficients.
1720
static const double coefficients0[6] = \
1721
{0.47140452, 0.23094011, 0.13333333, 0.00000000, 0.18856181, -0.16329932};
1722
1723
// Tables of derivatives of the polynomial base (transpose).
1724
static const double dmats0[6][6] = \
1725
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1726
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1727
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1728
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1729
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000},
1730
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1731
1732
static const double dmats1[6][6] = \
1733
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1734
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1735
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1736
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000},
1737
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000},
1738
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000}};
1739
1740
// Compute reference derivatives.
1741
// Declare pointer to array of derivatives on FIAT element.
1742
double *derivatives = new double[num_derivatives];
1743
for (unsigned int r = 0; r < num_derivatives; r++)
1744
{
1745
derivatives[r] = 0.00000000;
1746
}// end loop over 'r'
1747
1748
// Declare derivative matrix (of polynomial basis).
1749
double dmats[6][6] = \
1750
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1751
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1752
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1753
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1754
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1755
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1756
1757
// Declare (auxiliary) derivative matrix (of polynomial basis).
1758
double dmats_old[6][6] = \
1759
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1760
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1761
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1762
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1763
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1764
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1765
1766
// Loop possible derivatives.
1767
for (unsigned int r = 0; r < num_derivatives; r++)
1768
{
1769
// Resetting dmats values to compute next derivative.
1770
for (unsigned int t = 0; t < 6; t++)
1771
{
1772
for (unsigned int u = 0; u < 6; u++)
1773
{
1774
dmats[t][u] = 0.00000000;
1775
if (t == u)
1776
{
1777
dmats[t][u] = 1.00000000;
1778
}
1779
1780
}// end loop over 'u'
1781
}// end loop over 't'
1782
1783
// Looping derivative order to generate dmats.
1784
for (unsigned int s = 0; s < n; s++)
1785
{
1786
// Updating dmats_old with new values and resetting dmats.
1787
for (unsigned int t = 0; t < 6; t++)
1788
{
1789
for (unsigned int u = 0; u < 6; u++)
1790
{
1791
dmats_old[t][u] = dmats[t][u];
1792
dmats[t][u] = 0.00000000;
1793
}// end loop over 'u'
1794
}// end loop over 't'
1795
1796
// Update dmats using an inner product.
1797
if (combinations[r][s] == 0)
1798
{
1799
for (unsigned int t = 0; t < 6; t++)
1800
{
1801
for (unsigned int u = 0; u < 6; u++)
1802
{
1803
for (unsigned int tu = 0; tu < 6; tu++)
1804
{
1805
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1806
}// end loop over 'tu'
1807
}// end loop over 'u'
1808
}// end loop over 't'
1809
}
1810
1811
if (combinations[r][s] == 1)
1812
{
1813
for (unsigned int t = 0; t < 6; t++)
1814
{
1815
for (unsigned int u = 0; u < 6; u++)
1816
{
1817
for (unsigned int tu = 0; tu < 6; tu++)
1818
{
1819
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1820
}// end loop over 'tu'
1821
}// end loop over 'u'
1822
}// end loop over 't'
1823
}
1824
1825
}// end loop over 's'
1826
for (unsigned int s = 0; s < 6; s++)
1827
{
1828
for (unsigned int t = 0; t < 6; t++)
1829
{
1830
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1831
}// end loop over 't'
1832
}// end loop over 's'
1833
}// end loop over 'r'
1834
1835
// Transform derivatives back to physical element
1836
for (unsigned int r = 0; r < num_derivatives; r++)
1837
{
1838
for (unsigned int s = 0; s < num_derivatives; s++)
1839
{
1840
values[r] += transform[r][s]*derivatives[s];
1841
}// end loop over 's'
1842
}// end loop over 'r'
1843
1844
// Delete pointer to array of derivatives on FIAT element
1845
delete [] derivatives;
1846
1847
// Delete pointer to array of combinations of derivatives and transform
1848
for (unsigned int r = 0; r < num_derivatives; r++)
1849
{
1850
delete [] combinations[r];
1851
}// end loop over 'r'
1852
delete [] combinations;
1853
for (unsigned int r = 0; r < num_derivatives; r++)
1854
{
1855
delete [] transform[r];
1856
}// end loop over 'r'
1857
delete [] transform;
1858
break;
1859
}
1860
case 4:
1861
{
1862
1863
// Array of basisvalues.
1864
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
1865
1866
// Declare helper variables.
1867
unsigned int rr = 0;
1868
unsigned int ss = 0;
1869
unsigned int tt = 0;
1870
double tmp5 = 0.00000000;
1871
double tmp6 = 0.00000000;
1872
double tmp7 = 0.00000000;
1873
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1874
double tmp1 = (1.00000000 - Y)/2.00000000;
1875
double tmp2 = tmp1*tmp1;
1876
1877
// Compute basisvalues.
1878
basisvalues[0] = 1.00000000;
1879
basisvalues[1] = tmp0;
1880
for (unsigned int r = 1; r < 2; r++)
1881
{
1882
rr = (r + 1)*((r + 1) + 1)/2;
1883
ss = r*(r + 1)/2;
1884
tt = (r - 1)*((r - 1) + 1)/2;
1885
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
1886
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
1887
}// end loop over 'r'
1888
for (unsigned int r = 0; r < 2; r++)
1889
{
1890
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1891
ss = r*(r + 1)/2;
1892
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1893
}// end loop over 'r'
1894
for (unsigned int r = 0; r < 1; r++)
1895
{
1896
for (unsigned int s = 1; s < 2 - r; s++)
1897
{
1898
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
1899
ss = (r + s)*(r + s + 1)/2 + s;
1900
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
1901
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1902
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1903
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
1904
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
1905
}// end loop over 's'
1906
}// end loop over 'r'
1907
for (unsigned int r = 0; r < 3; r++)
1908
{
1909
for (unsigned int s = 0; s < 3 - r; s++)
1910
{
1911
rr = (r + s)*(r + s + 1)/2 + s;
1912
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1913
}// end loop over 's'
1914
}// end loop over 'r'
1915
1916
// Table(s) of coefficients.
1917
static const double coefficients0[6] = \
1918
{0.47140452, -0.23094011, 0.13333333, 0.00000000, -0.18856181, -0.16329932};
1919
1920
// Tables of derivatives of the polynomial base (transpose).
1921
static const double dmats0[6][6] = \
1922
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1923
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1924
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1925
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1926
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000},
1927
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1928
1929
static const double dmats1[6][6] = \
1930
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1931
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1932
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1933
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000},
1934
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000},
1935
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000}};
1936
1937
// Compute reference derivatives.
1938
// Declare pointer to array of derivatives on FIAT element.
1939
double *derivatives = new double[num_derivatives];
1940
for (unsigned int r = 0; r < num_derivatives; r++)
1941
{
1942
derivatives[r] = 0.00000000;
1943
}// end loop over 'r'
1944
1945
// Declare derivative matrix (of polynomial basis).
1946
double dmats[6][6] = \
1947
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1948
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1949
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1950
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1951
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1952
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1953
1954
// Declare (auxiliary) derivative matrix (of polynomial basis).
1955
double dmats_old[6][6] = \
1956
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1957
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
1958
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
1959
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
1960
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
1961
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1962
1963
// Loop possible derivatives.
1964
for (unsigned int r = 0; r < num_derivatives; r++)
1965
{
1966
// Resetting dmats values to compute next derivative.
1967
for (unsigned int t = 0; t < 6; t++)
1968
{
1969
for (unsigned int u = 0; u < 6; u++)
1970
{
1971
dmats[t][u] = 0.00000000;
1972
if (t == u)
1973
{
1974
dmats[t][u] = 1.00000000;
1975
}
1976
1977
}// end loop over 'u'
1978
}// end loop over 't'
1979
1980
// Looping derivative order to generate dmats.
1981
for (unsigned int s = 0; s < n; s++)
1982
{
1983
// Updating dmats_old with new values and resetting dmats.
1984
for (unsigned int t = 0; t < 6; t++)
1985
{
1986
for (unsigned int u = 0; u < 6; u++)
1987
{
1988
dmats_old[t][u] = dmats[t][u];
1989
dmats[t][u] = 0.00000000;
1990
}// end loop over 'u'
1991
}// end loop over 't'
1992
1993
// Update dmats using an inner product.
1994
if (combinations[r][s] == 0)
1995
{
1996
for (unsigned int t = 0; t < 6; t++)
1997
{
1998
for (unsigned int u = 0; u < 6; u++)
1999
{
2000
for (unsigned int tu = 0; tu < 6; tu++)
2001
{
2002
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2003
}// end loop over 'tu'
2004
}// end loop over 'u'
2005
}// end loop over 't'
2006
}
2007
2008
if (combinations[r][s] == 1)
2009
{
2010
for (unsigned int t = 0; t < 6; t++)
2011
{
2012
for (unsigned int u = 0; u < 6; u++)
2013
{
2014
for (unsigned int tu = 0; tu < 6; tu++)
2015
{
2016
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2017
}// end loop over 'tu'
2018
}// end loop over 'u'
2019
}// end loop over 't'
2020
}
2021
2022
}// end loop over 's'
2023
for (unsigned int s = 0; s < 6; s++)
2024
{
2025
for (unsigned int t = 0; t < 6; t++)
2026
{
2027
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2028
}// end loop over 't'
2029
}// end loop over 's'
2030
}// end loop over 'r'
2031
2032
// Transform derivatives back to physical element
2033
for (unsigned int r = 0; r < num_derivatives; r++)
2034
{
2035
for (unsigned int s = 0; s < num_derivatives; s++)
2036
{
2037
values[r] += transform[r][s]*derivatives[s];
2038
}// end loop over 's'
2039
}// end loop over 'r'
2040
2041
// Delete pointer to array of derivatives on FIAT element
2042
delete [] derivatives;
2043
2044
// Delete pointer to array of combinations of derivatives and transform
2045
for (unsigned int r = 0; r < num_derivatives; r++)
2046
{
2047
delete [] combinations[r];
2048
}// end loop over 'r'
2049
delete [] combinations;
2050
for (unsigned int r = 0; r < num_derivatives; r++)
2051
{
2052
delete [] transform[r];
2053
}// end loop over 'r'
2054
delete [] transform;
2055
break;
2056
}
2057
case 5:
2058
{
2059
2060
// Array of basisvalues.
2061
double basisvalues[6] = {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000};
2062
2063
// Declare helper variables.
2064
unsigned int rr = 0;
2065
unsigned int ss = 0;
2066
unsigned int tt = 0;
2067
double tmp5 = 0.00000000;
2068
double tmp6 = 0.00000000;
2069
double tmp7 = 0.00000000;
2070
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
2071
double tmp1 = (1.00000000 - Y)/2.00000000;
2072
double tmp2 = tmp1*tmp1;
2073
2074
// Compute basisvalues.
2075
basisvalues[0] = 1.00000000;
2076
basisvalues[1] = tmp0;
2077
for (unsigned int r = 1; r < 2; r++)
2078
{
2079
rr = (r + 1)*((r + 1) + 1)/2;
2080
ss = r*(r + 1)/2;
2081
tt = (r - 1)*((r - 1) + 1)/2;
2082
tmp5 = (1.00000000 + 2.00000000*r)/(1.00000000 + r);
2083
basisvalues[rr] = (basisvalues[ss]*tmp0*tmp5 - basisvalues[tt]*tmp2*r/(1.00000000 + r));
2084
}// end loop over 'r'
2085
for (unsigned int r = 0; r < 2; r++)
2086
{
2087
rr = (r + 1)*(r + 1 + 1)/2 + 1;
2088
ss = r*(r + 1)/2;
2089
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
2090
}// end loop over 'r'
2091
for (unsigned int r = 0; r < 1; r++)
2092
{
2093
for (unsigned int s = 1; s < 2 - r; s++)
2094
{
2095
rr = (r + s + 1)*(r + s + 1 + 1)/2 + s + 1;
2096
ss = (r + s)*(r + s + 1)/2 + s;
2097
tt = (r + s - 1)*(r + s - 1 + 1)/2 + s - 1;
2098
tmp5 = (2.00000000 + 2.00000000*r + 2.00000000*s)*(3.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2099
tmp6 = (1.00000000 + 4.00000000*r*r + 4.00000000*r)*(2.00000000 + 2.00000000*r + 2.00000000*s)/(2.00000000*(1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2100
tmp7 = (1.00000000 + s + 2.00000000*r)*(3.00000000 + 2.00000000*r + 2.00000000*s)*s/((1.00000000 + 2.00000000*r + 2.00000000*s)*(1.00000000 + s)*(2.00000000 + s + 2.00000000*r));
2101
basisvalues[rr] = (basisvalues[ss]*(tmp6 + Y*tmp5) - basisvalues[tt]*tmp7);
2102
}// end loop over 's'
2103
}// end loop over 'r'
2104
for (unsigned int r = 0; r < 3; r++)
2105
{
2106
for (unsigned int s = 0; s < 3 - r; s++)
2107
{
2108
rr = (r + s)*(r + s + 1)/2 + s;
2109
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
2110
}// end loop over 's'
2111
}// end loop over 'r'
2112
2113
// Table(s) of coefficients.
2114
static const double coefficients0[6] = \
2115
{0.47140452, 0.00000000, -0.26666667, -0.24343225, 0.00000000, 0.05443311};
2116
2117
// Tables of derivatives of the polynomial base (transpose).
2118
static const double dmats0[6][6] = \
2119
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2120
{4.89897949, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2121
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2122
{0.00000000, 9.48683298, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2123
{4.00000000, 0.00000000, 7.07106781, 0.00000000, 0.00000000, 0.00000000},
2124
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2125
2126
static const double dmats1[6][6] = \
2127
{{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2128
{2.44948974, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2129
{4.24264069, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2130
{2.58198890, 4.74341649, -0.91287093, 0.00000000, 0.00000000, 0.00000000},
2131
{2.00000000, 6.12372436, 3.53553391, 0.00000000, 0.00000000, 0.00000000},
2132
{-2.30940108, 0.00000000, 8.16496581, 0.00000000, 0.00000000, 0.00000000}};
2133
2134
// Compute reference derivatives.
2135
// Declare pointer to array of derivatives on FIAT element.
2136
double *derivatives = new double[num_derivatives];
2137
for (unsigned int r = 0; r < num_derivatives; r++)
2138
{
2139
derivatives[r] = 0.00000000;
2140
}// end loop over 'r'
2141
2142
// Declare derivative matrix (of polynomial basis).
2143
double dmats[6][6] = \
2144
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2145
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2146
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2147
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2148
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2149
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2150
2151
// Declare (auxiliary) derivative matrix (of polynomial basis).
2152
double dmats_old[6][6] = \
2153
{{1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2154
{0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000},
2155
{0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000},
2156
{0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000},
2157
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000},
2158
{0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2159
2160
// Loop possible derivatives.
2161
for (unsigned int r = 0; r < num_derivatives; r++)
2162
{
2163
// Resetting dmats values to compute next derivative.
2164
for (unsigned int t = 0; t < 6; t++)
2165
{
2166
for (unsigned int u = 0; u < 6; u++)
2167
{
2168
dmats[t][u] = 0.00000000;
2169
if (t == u)
2170
{
2171
dmats[t][u] = 1.00000000;
2172
}
2173
2174
}// end loop over 'u'
2175
}// end loop over 't'
2176
2177
// Looping derivative order to generate dmats.
2178
for (unsigned int s = 0; s < n; s++)
2179
{
2180
// Updating dmats_old with new values and resetting dmats.
2181
for (unsigned int t = 0; t < 6; t++)
2182
{
2183
for (unsigned int u = 0; u < 6; u++)
2184
{
2185
dmats_old[t][u] = dmats[t][u];
2186
dmats[t][u] = 0.00000000;
2187
}// end loop over 'u'
2188
}// end loop over 't'
2189
2190
// Update dmats using an inner product.
2191
if (combinations[r][s] == 0)
2192
{
2193
for (unsigned int t = 0; t < 6; t++)
2194
{
2195
for (unsigned int u = 0; u < 6; u++)
2196
{
2197
for (unsigned int tu = 0; tu < 6; tu++)
2198
{
2199
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2200
}// end loop over 'tu'
2201
}// end loop over 'u'
2202
}// end loop over 't'
2203
}
2204
2205
if (combinations[r][s] == 1)
2206
{
2207
for (unsigned int t = 0; t < 6; t++)
2208
{
2209
for (unsigned int u = 0; u < 6; u++)
2210
{
2211
for (unsigned int tu = 0; tu < 6; tu++)
2212
{
2213
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2214
}// end loop over 'tu'
2215
}// end loop over 'u'
2216
}// end loop over 't'
2217
}
2218
2219
}// end loop over 's'
2220
for (unsigned int s = 0; s < 6; s++)
2221
{
2222
for (unsigned int t = 0; t < 6; t++)
2223
{
2224
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2225
}// end loop over 't'
2226
}// end loop over 's'
2227
}// end loop over 'r'
2228
2229
// Transform derivatives back to physical element
2230
for (unsigned int r = 0; r < num_derivatives; r++)
2231
{
2232
for (unsigned int s = 0; s < num_derivatives; s++)
2233
{
2234
values[r] += transform[r][s]*derivatives[s];
2235
}// end loop over 's'
2236
}// end loop over 'r'
2237
2238
// Delete pointer to array of derivatives on FIAT element
2239
delete [] derivatives;
2240
2241
// Delete pointer to array of combinations of derivatives and transform
2242
for (unsigned int r = 0; r < num_derivatives; r++)
2243
{
2244
delete [] combinations[r];
2245
}// end loop over 'r'
2246
delete [] combinations;
2247
for (unsigned int r = 0; r < num_derivatives; r++)
2248
{
2249
delete [] transform[r];
2250
}// end loop over 'r'
2251
delete [] transform;
2252
break;
2253
}
2254
}
2255
935
2256
}
936
2257
937
2258
/// Evaluate order n derivatives of all basis functions at given point in cell
1954
3275
}// end loop over 'r'
1955
3276
1956
3277
// Compute element tensor using UFL quadrature representation
1957
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
3278
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
1958
3279
1959
3280
// Loop quadrature points for integral.
1960
3281
// Number of operations to compute element tensor for following IP loop = 6192
2064
3385
}// end loop over 'r'
2065
3386
2066
3387
// Compute element tensor using UFL quadrature representation
2067
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
3388
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
2068
3389
2069
3390
// Loop quadrature points for integral.
2070
3391
// Number of operations to compute element tensor for following IP loop = 408
Loggerhead is a web-based interface for Breezy
Version: 2.0.1