~ubuntu-branches/ubuntu/trusty/ffc/trusty
» Revision
8
: test/regression/references/HyperElasticity.h
« back to all changes in this revision
Viewing changes to test/regression/references/HyperElasticity.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/HyperElasticity.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
//
99
99
// Reset values.
100
100
*values = 0.00000000;
101
101
102
// Map degree of freedom to element degree of freedom
103
const unsigned int dof = i;
104
105
102
// Array of basisvalues.
106
103
double basisvalues[1] = {0.00000000};
107
104
111
108
basisvalues[0] = 1.00000000;
112
109
113
110
// Table(s) of coefficients.
114
static const double coefficients0[1][1] = \
115
{{1.00000000}};
111
static const double coefficients0[1] = \
112
{1.00000000};
116
113
117
114
// Compute value(s).
118
115
for (unsigned int r = 0; r < 1; r++)
119
116
{
120
*values += coefficients0[dof][r]*basisvalues[r];
117
*values += coefficients0[r]*basisvalues[r];
121
118
}// end loop over 'r'
122
119
}
123
120
245
242
values[r] = 0.00000000;
246
243
}// end loop over 'r'
247
244
248
// Map degree of freedom to element degree of freedom
249
const unsigned int dof = i;
250
245
251
246
// Array of basisvalues.
252
247
double basisvalues[1] = {0.00000000};
257
252
basisvalues[0] = 1.00000000;
258
253
259
254
// Table(s) of coefficients.
260
static const double coefficients0[1][1] = \
261
{{1.00000000}};
255
static const double coefficients0[1] = \
256
{1.00000000};
262
257
263
258
// Tables of derivatives of the polynomial base (transpose).
264
259
static const double dmats0[1][1] = \
364
359
{
365
360
for (unsigned int t = 0; t < 1; t++)
366
361
{
367
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
362
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
368
363
}// end loop over 't'
369
364
}// end loop over 's'
370
365
}// end loop over 'r'
572
567
573
568
// Reset values.
574
569
*values = 0.00000000;
575
576
// Map degree of freedom to element degree of freedom
577
const unsigned int dof = i;
578
579
// Array of basisvalues.
580
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
581
582
// Declare helper variables.
583
unsigned int rr = 0;
584
unsigned int ss = 0;
585
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
586
587
// Compute basisvalues.
588
basisvalues[0] = 1.00000000;
589
basisvalues[1] = tmp0;
590
for (unsigned int r = 0; r < 1; r++)
591
{
592
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
593
ss = r*(r + 1)*(r + 2)/6;
594
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
595
}// end loop over 'r'
596
for (unsigned int r = 0; r < 1; r++)
597
{
598
for (unsigned int s = 0; s < 1 - r; s++)
599
{
600
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
601
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
602
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
603
}// end loop over 's'
604
}// end loop over 'r'
605
for (unsigned int r = 0; r < 2; r++)
606
{
607
for (unsigned int s = 0; s < 2 - r; s++)
608
{
609
for (unsigned int t = 0; t < 2 - r - s; t++)
610
{
611
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
612
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
613
}// end loop over 't'
614
}// end loop over 's'
615
}// end loop over 'r'
616
617
// Table(s) of coefficients.
618
static const double coefficients0[4][4] = \
619
{{0.28867513, -0.18257419, -0.10540926, -0.07453560},
620
{0.28867513, 0.18257419, -0.10540926, -0.07453560},
621
{0.28867513, 0.00000000, 0.21081851, -0.07453560},
622
{0.28867513, 0.00000000, 0.00000000, 0.22360680}};
623
624
// Compute value(s).
625
for (unsigned int r = 0; r < 4; r++)
626
{
627
*values += coefficients0[dof][r]*basisvalues[r];
628
}// end loop over 'r'
570
switch (i)
571
{
572
case 0:
573
{
574
575
// Array of basisvalues.
576
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
577
578
// Declare helper variables.
579
unsigned int rr = 0;
580
unsigned int ss = 0;
581
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
582
583
// Compute basisvalues.
584
basisvalues[0] = 1.00000000;
585
basisvalues[1] = tmp0;
586
for (unsigned int r = 0; r < 1; r++)
587
{
588
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
589
ss = r*(r + 1)*(r + 2)/6;
590
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
591
}// end loop over 'r'
592
for (unsigned int r = 0; r < 1; r++)
593
{
594
for (unsigned int s = 0; s < 1 - r; s++)
595
{
596
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
597
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
598
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
599
}// end loop over 's'
600
}// end loop over 'r'
601
for (unsigned int r = 0; r < 2; r++)
602
{
603
for (unsigned int s = 0; s < 2 - r; s++)
604
{
605
for (unsigned int t = 0; t < 2 - r - s; t++)
606
{
607
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
608
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
609
}// end loop over 't'
610
}// end loop over 's'
611
}// end loop over 'r'
612
613
// Table(s) of coefficients.
614
static const double coefficients0[4] = \
615
{0.28867513, -0.18257419, -0.10540926, -0.07453560};
616
617
// Compute value(s).
618
for (unsigned int r = 0; r < 4; r++)
619
{
620
*values += coefficients0[r]*basisvalues[r];
621
}// end loop over 'r'
622
break;
623
}
624
case 1:
625
{
626
627
// Array of basisvalues.
628
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
629
630
// Declare helper variables.
631
unsigned int rr = 0;
632
unsigned int ss = 0;
633
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
634
635
// Compute basisvalues.
636
basisvalues[0] = 1.00000000;
637
basisvalues[1] = tmp0;
638
for (unsigned int r = 0; r < 1; r++)
639
{
640
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
641
ss = r*(r + 1)*(r + 2)/6;
642
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
643
}// end loop over 'r'
644
for (unsigned int r = 0; r < 1; r++)
645
{
646
for (unsigned int s = 0; s < 1 - r; s++)
647
{
648
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
649
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
650
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
651
}// end loop over 's'
652
}// end loop over 'r'
653
for (unsigned int r = 0; r < 2; r++)
654
{
655
for (unsigned int s = 0; s < 2 - r; s++)
656
{
657
for (unsigned int t = 0; t < 2 - r - s; t++)
658
{
659
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
660
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
661
}// end loop over 't'
662
}// end loop over 's'
663
}// end loop over 'r'
664
665
// Table(s) of coefficients.
666
static const double coefficients0[4] = \
667
{0.28867513, 0.18257419, -0.10540926, -0.07453560};
668
669
// Compute value(s).
670
for (unsigned int r = 0; r < 4; r++)
671
{
672
*values += coefficients0[r]*basisvalues[r];
673
}// end loop over 'r'
674
break;
675
}
676
case 2:
677
{
678
679
// Array of basisvalues.
680
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
681
682
// Declare helper variables.
683
unsigned int rr = 0;
684
unsigned int ss = 0;
685
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
686
687
// Compute basisvalues.
688
basisvalues[0] = 1.00000000;
689
basisvalues[1] = tmp0;
690
for (unsigned int r = 0; r < 1; r++)
691
{
692
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
693
ss = r*(r + 1)*(r + 2)/6;
694
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
695
}// end loop over 'r'
696
for (unsigned int r = 0; r < 1; r++)
697
{
698
for (unsigned int s = 0; s < 1 - r; s++)
699
{
700
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
701
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
702
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
703
}// end loop over 's'
704
}// end loop over 'r'
705
for (unsigned int r = 0; r < 2; r++)
706
{
707
for (unsigned int s = 0; s < 2 - r; s++)
708
{
709
for (unsigned int t = 0; t < 2 - r - s; t++)
710
{
711
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
712
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
713
}// end loop over 't'
714
}// end loop over 's'
715
}// end loop over 'r'
716
717
// Table(s) of coefficients.
718
static const double coefficients0[4] = \
719
{0.28867513, 0.00000000, 0.21081851, -0.07453560};
720
721
// Compute value(s).
722
for (unsigned int r = 0; r < 4; r++)
723
{
724
*values += coefficients0[r]*basisvalues[r];
725
}// end loop over 'r'
726
break;
727
}
728
case 3:
729
{
730
731
// Array of basisvalues.
732
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
733
734
// Declare helper variables.
735
unsigned int rr = 0;
736
unsigned int ss = 0;
737
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
738
739
// Compute basisvalues.
740
basisvalues[0] = 1.00000000;
741
basisvalues[1] = tmp0;
742
for (unsigned int r = 0; r < 1; r++)
743
{
744
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
745
ss = r*(r + 1)*(r + 2)/6;
746
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
747
}// end loop over 'r'
748
for (unsigned int r = 0; r < 1; r++)
749
{
750
for (unsigned int s = 0; s < 1 - r; s++)
751
{
752
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
753
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
754
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
755
}// end loop over 's'
756
}// end loop over 'r'
757
for (unsigned int r = 0; r < 2; r++)
758
{
759
for (unsigned int s = 0; s < 2 - r; s++)
760
{
761
for (unsigned int t = 0; t < 2 - r - s; t++)
762
{
763
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
764
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
765
}// end loop over 't'
766
}// end loop over 's'
767
}// end loop over 'r'
768
769
// Table(s) of coefficients.
770
static const double coefficients0[4] = \
771
{0.28867513, 0.00000000, 0.00000000, 0.22360680};
772
773
// Compute value(s).
774
for (unsigned int r = 0; r < 4; r++)
775
{
776
*values += coefficients0[r]*basisvalues[r];
777
}// end loop over 'r'
778
break;
779
}
780
}
781
629
782
}
630
783
631
784
/// Evaluate all basis functions at given point in cell
765
918
values[r] = 0.00000000;
766
919
}// end loop over 'r'
767
920
768
// Map degree of freedom to element degree of freedom
769
const unsigned int dof = i;
770
771
// Array of basisvalues.
772
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
773
774
// Declare helper variables.
775
unsigned int rr = 0;
776
unsigned int ss = 0;
777
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
778
779
// Compute basisvalues.
780
basisvalues[0] = 1.00000000;
781
basisvalues[1] = tmp0;
782
for (unsigned int r = 0; r < 1; r++)
783
{
784
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
785
ss = r*(r + 1)*(r + 2)/6;
786
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
787
}// end loop over 'r'
788
for (unsigned int r = 0; r < 1; r++)
789
{
790
for (unsigned int s = 0; s < 1 - r; s++)
791
{
792
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
793
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
794
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
795
}// end loop over 's'
796
}// end loop over 'r'
797
for (unsigned int r = 0; r < 2; r++)
798
{
799
for (unsigned int s = 0; s < 2 - r; s++)
800
{
801
for (unsigned int t = 0; t < 2 - r - s; t++)
802
{
803
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
804
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
805
}// end loop over 't'
806
}// end loop over 's'
807
}// end loop over 'r'
808
809
// Table(s) of coefficients.
810
static const double coefficients0[4][4] = \
811
{{0.28867513, -0.18257419, -0.10540926, -0.07453560},
812
{0.28867513, 0.18257419, -0.10540926, -0.07453560},
813
{0.28867513, 0.00000000, 0.21081851, -0.07453560},
814
{0.28867513, 0.00000000, 0.00000000, 0.22360680}};
815
816
// Tables of derivatives of the polynomial base (transpose).
817
static const double dmats0[4][4] = \
818
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
819
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
820
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
821
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
822
823
static const double dmats1[4][4] = \
824
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
825
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
826
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
827
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
828
829
static const double dmats2[4][4] = \
830
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
831
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
832
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
833
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
834
835
// Compute reference derivatives.
836
// Declare pointer to array of derivatives on FIAT element.
837
double *derivatives = new double[num_derivatives];
838
for (unsigned int r = 0; r < num_derivatives; r++)
839
{
840
derivatives[r] = 0.00000000;
841
}// end loop over 'r'
842
843
// Declare derivative matrix (of polynomial basis).
844
double dmats[4][4] = \
845
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
846
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
847
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
848
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
849
850
// Declare (auxiliary) derivative matrix (of polynomial basis).
851
double dmats_old[4][4] = \
852
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
853
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
854
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
855
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
856
857
// Loop possible derivatives.
858
for (unsigned int r = 0; r < num_derivatives; r++)
859
{
860
// Resetting dmats values to compute next derivative.
861
for (unsigned int t = 0; t < 4; t++)
862
{
863
for (unsigned int u = 0; u < 4; u++)
864
{
865
dmats[t][u] = 0.00000000;
866
if (t == u)
867
{
868
dmats[t][u] = 1.00000000;
869
}
870
871
}// end loop over 'u'
872
}// end loop over 't'
873
874
// Looping derivative order to generate dmats.
875
for (unsigned int s = 0; s < n; s++)
876
{
877
// Updating dmats_old with new values and resetting dmats.
878
for (unsigned int t = 0; t < 4; t++)
879
{
880
for (unsigned int u = 0; u < 4; u++)
881
{
882
dmats_old[t][u] = dmats[t][u];
883
dmats[t][u] = 0.00000000;
884
}// end loop over 'u'
885
}// end loop over 't'
886
887
// Update dmats using an inner product.
888
if (combinations[r][s] == 0)
889
{
890
for (unsigned int t = 0; t < 4; t++)
891
{
892
for (unsigned int u = 0; u < 4; u++)
893
{
894
for (unsigned int tu = 0; tu < 4; tu++)
895
{
896
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
897
}// end loop over 'tu'
898
}// end loop over 'u'
899
}// end loop over 't'
900
}
901
902
if (combinations[r][s] == 1)
903
{
904
for (unsigned int t = 0; t < 4; t++)
905
{
906
for (unsigned int u = 0; u < 4; u++)
907
{
908
for (unsigned int tu = 0; tu < 4; tu++)
909
{
910
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
911
}// end loop over 'tu'
912
}// end loop over 'u'
913
}// end loop over 't'
914
}
915
916
if (combinations[r][s] == 2)
917
{
918
for (unsigned int t = 0; t < 4; t++)
919
{
920
for (unsigned int u = 0; u < 4; u++)
921
{
922
for (unsigned int tu = 0; tu < 4; tu++)
923
{
924
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
925
}// end loop over 'tu'
926
}// end loop over 'u'
927
}// end loop over 't'
928
}
929
930
}// end loop over 's'
931
for (unsigned int s = 0; s < 4; s++)
932
{
933
for (unsigned int t = 0; t < 4; t++)
934
{
935
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
936
}// end loop over 't'
937
}// end loop over 's'
938
}// end loop over 'r'
939
940
// Transform derivatives back to physical element
941
for (unsigned int r = 0; r < num_derivatives; r++)
942
{
943
for (unsigned int s = 0; s < num_derivatives; s++)
944
{
945
values[r] += transform[r][s]*derivatives[s];
946
}// end loop over 's'
947
}// end loop over 'r'
948
949
// Delete pointer to array of derivatives on FIAT element
950
delete [] derivatives;
951
952
// Delete pointer to array of combinations of derivatives and transform
953
for (unsigned int r = 0; r < num_derivatives; r++)
954
{
955
delete [] combinations[r];
956
}// end loop over 'r'
957
delete [] combinations;
958
for (unsigned int r = 0; r < num_derivatives; r++)
959
{
960
delete [] transform[r];
961
}// end loop over 'r'
962
delete [] transform;
921
switch (i)
922
{
923
case 0:
924
{
925
926
// Array of basisvalues.
927
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
928
929
// Declare helper variables.
930
unsigned int rr = 0;
931
unsigned int ss = 0;
932
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
933
934
// Compute basisvalues.
935
basisvalues[0] = 1.00000000;
936
basisvalues[1] = tmp0;
937
for (unsigned int r = 0; r < 1; r++)
938
{
939
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
940
ss = r*(r + 1)*(r + 2)/6;
941
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
942
}// end loop over 'r'
943
for (unsigned int r = 0; r < 1; r++)
944
{
945
for (unsigned int s = 0; s < 1 - r; s++)
946
{
947
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
948
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
949
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
950
}// end loop over 's'
951
}// end loop over 'r'
952
for (unsigned int r = 0; r < 2; r++)
953
{
954
for (unsigned int s = 0; s < 2 - r; s++)
955
{
956
for (unsigned int t = 0; t < 2 - r - s; t++)
957
{
958
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
959
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
960
}// end loop over 't'
961
}// end loop over 's'
962
}// end loop over 'r'
963
964
// Table(s) of coefficients.
965
static const double coefficients0[4] = \
966
{0.28867513, -0.18257419, -0.10540926, -0.07453560};
967
968
// Tables of derivatives of the polynomial base (transpose).
969
static const double dmats0[4][4] = \
970
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
971
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
972
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
973
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
974
975
static const double dmats1[4][4] = \
976
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
977
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
978
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
979
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
980
981
static const double dmats2[4][4] = \
982
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
983
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
984
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
985
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
986
987
// Compute reference derivatives.
988
// Declare pointer to array of derivatives on FIAT element.
989
double *derivatives = new double[num_derivatives];
990
for (unsigned int r = 0; r < num_derivatives; r++)
991
{
992
derivatives[r] = 0.00000000;
993
}// end loop over 'r'
994
995
// Declare derivative matrix (of polynomial basis).
996
double dmats[4][4] = \
997
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
998
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
999
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1000
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1001
1002
// Declare (auxiliary) derivative matrix (of polynomial basis).
1003
double dmats_old[4][4] = \
1004
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1005
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1006
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1007
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1008
1009
// Loop possible derivatives.
1010
for (unsigned int r = 0; r < num_derivatives; r++)
1011
{
1012
// Resetting dmats values to compute next derivative.
1013
for (unsigned int t = 0; t < 4; t++)
1014
{
1015
for (unsigned int u = 0; u < 4; u++)
1016
{
1017
dmats[t][u] = 0.00000000;
1018
if (t == u)
1019
{
1020
dmats[t][u] = 1.00000000;
1021
}
1022
1023
}// end loop over 'u'
1024
}// end loop over 't'
1025
1026
// Looping derivative order to generate dmats.
1027
for (unsigned int s = 0; s < n; s++)
1028
{
1029
// Updating dmats_old with new values and resetting dmats.
1030
for (unsigned int t = 0; t < 4; t++)
1031
{
1032
for (unsigned int u = 0; u < 4; u++)
1033
{
1034
dmats_old[t][u] = dmats[t][u];
1035
dmats[t][u] = 0.00000000;
1036
}// end loop over 'u'
1037
}// end loop over 't'
1038
1039
// Update dmats using an inner product.
1040
if (combinations[r][s] == 0)
1041
{
1042
for (unsigned int t = 0; t < 4; t++)
1043
{
1044
for (unsigned int u = 0; u < 4; u++)
1045
{
1046
for (unsigned int tu = 0; tu < 4; tu++)
1047
{
1048
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1049
}// end loop over 'tu'
1050
}// end loop over 'u'
1051
}// end loop over 't'
1052
}
1053
1054
if (combinations[r][s] == 1)
1055
{
1056
for (unsigned int t = 0; t < 4; t++)
1057
{
1058
for (unsigned int u = 0; u < 4; u++)
1059
{
1060
for (unsigned int tu = 0; tu < 4; tu++)
1061
{
1062
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1063
}// end loop over 'tu'
1064
}// end loop over 'u'
1065
}// end loop over 't'
1066
}
1067
1068
if (combinations[r][s] == 2)
1069
{
1070
for (unsigned int t = 0; t < 4; t++)
1071
{
1072
for (unsigned int u = 0; u < 4; u++)
1073
{
1074
for (unsigned int tu = 0; tu < 4; tu++)
1075
{
1076
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1077
}// end loop over 'tu'
1078
}// end loop over 'u'
1079
}// end loop over 't'
1080
}
1081
1082
}// end loop over 's'
1083
for (unsigned int s = 0; s < 4; s++)
1084
{
1085
for (unsigned int t = 0; t < 4; t++)
1086
{
1087
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1088
}// end loop over 't'
1089
}// end loop over 's'
1090
}// end loop over 'r'
1091
1092
// Transform derivatives back to physical element
1093
for (unsigned int r = 0; r < num_derivatives; r++)
1094
{
1095
for (unsigned int s = 0; s < num_derivatives; s++)
1096
{
1097
values[r] += transform[r][s]*derivatives[s];
1098
}// end loop over 's'
1099
}// end loop over 'r'
1100
1101
// Delete pointer to array of derivatives on FIAT element
1102
delete [] derivatives;
1103
1104
// Delete pointer to array of combinations of derivatives and transform
1105
for (unsigned int r = 0; r < num_derivatives; r++)
1106
{
1107
delete [] combinations[r];
1108
}// end loop over 'r'
1109
delete [] combinations;
1110
for (unsigned int r = 0; r < num_derivatives; r++)
1111
{
1112
delete [] transform[r];
1113
}// end loop over 'r'
1114
delete [] transform;
1115
break;
1116
}
1117
case 1:
1118
{
1119
1120
// Array of basisvalues.
1121
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
1122
1123
// Declare helper variables.
1124
unsigned int rr = 0;
1125
unsigned int ss = 0;
1126
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
1127
1128
// Compute basisvalues.
1129
basisvalues[0] = 1.00000000;
1130
basisvalues[1] = tmp0;
1131
for (unsigned int r = 0; r < 1; r++)
1132
{
1133
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
1134
ss = r*(r + 1)*(r + 2)/6;
1135
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
1136
}// end loop over 'r'
1137
for (unsigned int r = 0; r < 1; r++)
1138
{
1139
for (unsigned int s = 0; s < 1 - r; s++)
1140
{
1141
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
1142
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
1143
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
1144
}// end loop over 's'
1145
}// end loop over 'r'
1146
for (unsigned int r = 0; r < 2; r++)
1147
{
1148
for (unsigned int s = 0; s < 2 - r; s++)
1149
{
1150
for (unsigned int t = 0; t < 2 - r - s; t++)
1151
{
1152
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
1153
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
1154
}// end loop over 't'
1155
}// end loop over 's'
1156
}// end loop over 'r'
1157
1158
// Table(s) of coefficients.
1159
static const double coefficients0[4] = \
1160
{0.28867513, 0.18257419, -0.10540926, -0.07453560};
1161
1162
// Tables of derivatives of the polynomial base (transpose).
1163
static const double dmats0[4][4] = \
1164
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1165
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
1166
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1167
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1168
1169
static const double dmats1[4][4] = \
1170
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1171
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1172
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
1173
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1174
1175
static const double dmats2[4][4] = \
1176
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1177
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1178
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
1179
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
1180
1181
// Compute reference derivatives.
1182
// Declare pointer to array of derivatives on FIAT element.
1183
double *derivatives = new double[num_derivatives];
1184
for (unsigned int r = 0; r < num_derivatives; r++)
1185
{
1186
derivatives[r] = 0.00000000;
1187
}// end loop over 'r'
1188
1189
// Declare derivative matrix (of polynomial basis).
1190
double dmats[4][4] = \
1191
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1192
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1193
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1194
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1195
1196
// Declare (auxiliary) derivative matrix (of polynomial basis).
1197
double dmats_old[4][4] = \
1198
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1199
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1200
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1201
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1202
1203
// Loop possible derivatives.
1204
for (unsigned int r = 0; r < num_derivatives; r++)
1205
{
1206
// Resetting dmats values to compute next derivative.
1207
for (unsigned int t = 0; t < 4; t++)
1208
{
1209
for (unsigned int u = 0; u < 4; u++)
1210
{
1211
dmats[t][u] = 0.00000000;
1212
if (t == u)
1213
{
1214
dmats[t][u] = 1.00000000;
1215
}
1216
1217
}// end loop over 'u'
1218
}// end loop over 't'
1219
1220
// Looping derivative order to generate dmats.
1221
for (unsigned int s = 0; s < n; s++)
1222
{
1223
// Updating dmats_old with new values and resetting dmats.
1224
for (unsigned int t = 0; t < 4; t++)
1225
{
1226
for (unsigned int u = 0; u < 4; u++)
1227
{
1228
dmats_old[t][u] = dmats[t][u];
1229
dmats[t][u] = 0.00000000;
1230
}// end loop over 'u'
1231
}// end loop over 't'
1232
1233
// Update dmats using an inner product.
1234
if (combinations[r][s] == 0)
1235
{
1236
for (unsigned int t = 0; t < 4; t++)
1237
{
1238
for (unsigned int u = 0; u < 4; u++)
1239
{
1240
for (unsigned int tu = 0; tu < 4; tu++)
1241
{
1242
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1243
}// end loop over 'tu'
1244
}// end loop over 'u'
1245
}// end loop over 't'
1246
}
1247
1248
if (combinations[r][s] == 1)
1249
{
1250
for (unsigned int t = 0; t < 4; t++)
1251
{
1252
for (unsigned int u = 0; u < 4; u++)
1253
{
1254
for (unsigned int tu = 0; tu < 4; tu++)
1255
{
1256
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1257
}// end loop over 'tu'
1258
}// end loop over 'u'
1259
}// end loop over 't'
1260
}
1261
1262
if (combinations[r][s] == 2)
1263
{
1264
for (unsigned int t = 0; t < 4; t++)
1265
{
1266
for (unsigned int u = 0; u < 4; u++)
1267
{
1268
for (unsigned int tu = 0; tu < 4; tu++)
1269
{
1270
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1271
}// end loop over 'tu'
1272
}// end loop over 'u'
1273
}// end loop over 't'
1274
}
1275
1276
}// end loop over 's'
1277
for (unsigned int s = 0; s < 4; s++)
1278
{
1279
for (unsigned int t = 0; t < 4; t++)
1280
{
1281
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1282
}// end loop over 't'
1283
}// end loop over 's'
1284
}// end loop over 'r'
1285
1286
// Transform derivatives back to physical element
1287
for (unsigned int r = 0; r < num_derivatives; r++)
1288
{
1289
for (unsigned int s = 0; s < num_derivatives; s++)
1290
{
1291
values[r] += transform[r][s]*derivatives[s];
1292
}// end loop over 's'
1293
}// end loop over 'r'
1294
1295
// Delete pointer to array of derivatives on FIAT element
1296
delete [] derivatives;
1297
1298
// Delete pointer to array of combinations of derivatives and transform
1299
for (unsigned int r = 0; r < num_derivatives; r++)
1300
{
1301
delete [] combinations[r];
1302
}// end loop over 'r'
1303
delete [] combinations;
1304
for (unsigned int r = 0; r < num_derivatives; r++)
1305
{
1306
delete [] transform[r];
1307
}// end loop over 'r'
1308
delete [] transform;
1309
break;
1310
}
1311
case 2:
1312
{
1313
1314
// Array of basisvalues.
1315
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
1316
1317
// Declare helper variables.
1318
unsigned int rr = 0;
1319
unsigned int ss = 0;
1320
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
1321
1322
// Compute basisvalues.
1323
basisvalues[0] = 1.00000000;
1324
basisvalues[1] = tmp0;
1325
for (unsigned int r = 0; r < 1; r++)
1326
{
1327
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
1328
ss = r*(r + 1)*(r + 2)/6;
1329
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
1330
}// end loop over 'r'
1331
for (unsigned int r = 0; r < 1; r++)
1332
{
1333
for (unsigned int s = 0; s < 1 - r; s++)
1334
{
1335
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
1336
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
1337
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
1338
}// end loop over 's'
1339
}// end loop over 'r'
1340
for (unsigned int r = 0; r < 2; r++)
1341
{
1342
for (unsigned int s = 0; s < 2 - r; s++)
1343
{
1344
for (unsigned int t = 0; t < 2 - r - s; t++)
1345
{
1346
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
1347
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
1348
}// end loop over 't'
1349
}// end loop over 's'
1350
}// end loop over 'r'
1351
1352
// Table(s) of coefficients.
1353
static const double coefficients0[4] = \
1354
{0.28867513, 0.00000000, 0.21081851, -0.07453560};
1355
1356
// Tables of derivatives of the polynomial base (transpose).
1357
static const double dmats0[4][4] = \
1358
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1359
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
1360
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1361
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1362
1363
static const double dmats1[4][4] = \
1364
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1365
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1366
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
1367
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1368
1369
static const double dmats2[4][4] = \
1370
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1371
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1372
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
1373
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
1374
1375
// Compute reference derivatives.
1376
// Declare pointer to array of derivatives on FIAT element.
1377
double *derivatives = new double[num_derivatives];
1378
for (unsigned int r = 0; r < num_derivatives; r++)
1379
{
1380
derivatives[r] = 0.00000000;
1381
}// end loop over 'r'
1382
1383
// Declare derivative matrix (of polynomial basis).
1384
double dmats[4][4] = \
1385
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1386
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1387
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1388
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1389
1390
// Declare (auxiliary) derivative matrix (of polynomial basis).
1391
double dmats_old[4][4] = \
1392
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1393
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1394
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1395
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1396
1397
// Loop possible derivatives.
1398
for (unsigned int r = 0; r < num_derivatives; r++)
1399
{
1400
// Resetting dmats values to compute next derivative.
1401
for (unsigned int t = 0; t < 4; t++)
1402
{
1403
for (unsigned int u = 0; u < 4; u++)
1404
{
1405
dmats[t][u] = 0.00000000;
1406
if (t == u)
1407
{
1408
dmats[t][u] = 1.00000000;
1409
}
1410
1411
}// end loop over 'u'
1412
}// end loop over 't'
1413
1414
// Looping derivative order to generate dmats.
1415
for (unsigned int s = 0; s < n; s++)
1416
{
1417
// Updating dmats_old with new values and resetting dmats.
1418
for (unsigned int t = 0; t < 4; t++)
1419
{
1420
for (unsigned int u = 0; u < 4; u++)
1421
{
1422
dmats_old[t][u] = dmats[t][u];
1423
dmats[t][u] = 0.00000000;
1424
}// end loop over 'u'
1425
}// end loop over 't'
1426
1427
// Update dmats using an inner product.
1428
if (combinations[r][s] == 0)
1429
{
1430
for (unsigned int t = 0; t < 4; t++)
1431
{
1432
for (unsigned int u = 0; u < 4; u++)
1433
{
1434
for (unsigned int tu = 0; tu < 4; tu++)
1435
{
1436
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1437
}// end loop over 'tu'
1438
}// end loop over 'u'
1439
}// end loop over 't'
1440
}
1441
1442
if (combinations[r][s] == 1)
1443
{
1444
for (unsigned int t = 0; t < 4; t++)
1445
{
1446
for (unsigned int u = 0; u < 4; u++)
1447
{
1448
for (unsigned int tu = 0; tu < 4; tu++)
1449
{
1450
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1451
}// end loop over 'tu'
1452
}// end loop over 'u'
1453
}// end loop over 't'
1454
}
1455
1456
if (combinations[r][s] == 2)
1457
{
1458
for (unsigned int t = 0; t < 4; t++)
1459
{
1460
for (unsigned int u = 0; u < 4; u++)
1461
{
1462
for (unsigned int tu = 0; tu < 4; tu++)
1463
{
1464
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1465
}// end loop over 'tu'
1466
}// end loop over 'u'
1467
}// end loop over 't'
1468
}
1469
1470
}// end loop over 's'
1471
for (unsigned int s = 0; s < 4; s++)
1472
{
1473
for (unsigned int t = 0; t < 4; t++)
1474
{
1475
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1476
}// end loop over 't'
1477
}// end loop over 's'
1478
}// end loop over 'r'
1479
1480
// Transform derivatives back to physical element
1481
for (unsigned int r = 0; r < num_derivatives; r++)
1482
{
1483
for (unsigned int s = 0; s < num_derivatives; s++)
1484
{
1485
values[r] += transform[r][s]*derivatives[s];
1486
}// end loop over 's'
1487
}// end loop over 'r'
1488
1489
// Delete pointer to array of derivatives on FIAT element
1490
delete [] derivatives;
1491
1492
// Delete pointer to array of combinations of derivatives and transform
1493
for (unsigned int r = 0; r < num_derivatives; r++)
1494
{
1495
delete [] combinations[r];
1496
}// end loop over 'r'
1497
delete [] combinations;
1498
for (unsigned int r = 0; r < num_derivatives; r++)
1499
{
1500
delete [] transform[r];
1501
}// end loop over 'r'
1502
delete [] transform;
1503
break;
1504
}
1505
case 3:
1506
{
1507
1508
// Array of basisvalues.
1509
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
1510
1511
// Declare helper variables.
1512
unsigned int rr = 0;
1513
unsigned int ss = 0;
1514
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
1515
1516
// Compute basisvalues.
1517
basisvalues[0] = 1.00000000;
1518
basisvalues[1] = tmp0;
1519
for (unsigned int r = 0; r < 1; r++)
1520
{
1521
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
1522
ss = r*(r + 1)*(r + 2)/6;
1523
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
1524
}// end loop over 'r'
1525
for (unsigned int r = 0; r < 1; r++)
1526
{
1527
for (unsigned int s = 0; s < 1 - r; s++)
1528
{
1529
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
1530
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
1531
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
1532
}// end loop over 's'
1533
}// end loop over 'r'
1534
for (unsigned int r = 0; r < 2; r++)
1535
{
1536
for (unsigned int s = 0; s < 2 - r; s++)
1537
{
1538
for (unsigned int t = 0; t < 2 - r - s; t++)
1539
{
1540
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
1541
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
1542
}// end loop over 't'
1543
}// end loop over 's'
1544
}// end loop over 'r'
1545
1546
// Table(s) of coefficients.
1547
static const double coefficients0[4] = \
1548
{0.28867513, 0.00000000, 0.00000000, 0.22360680};
1549
1550
// Tables of derivatives of the polynomial base (transpose).
1551
static const double dmats0[4][4] = \
1552
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1553
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
1554
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1555
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1556
1557
static const double dmats1[4][4] = \
1558
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1559
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1560
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
1561
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1562
1563
static const double dmats2[4][4] = \
1564
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1565
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1566
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
1567
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
1568
1569
// Compute reference derivatives.
1570
// Declare pointer to array of derivatives on FIAT element.
1571
double *derivatives = new double[num_derivatives];
1572
for (unsigned int r = 0; r < num_derivatives; r++)
1573
{
1574
derivatives[r] = 0.00000000;
1575
}// end loop over 'r'
1576
1577
// Declare derivative matrix (of polynomial basis).
1578
double dmats[4][4] = \
1579
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1580
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1581
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1582
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1583
1584
// Declare (auxiliary) derivative matrix (of polynomial basis).
1585
double dmats_old[4][4] = \
1586
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1587
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1588
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1589
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1590
1591
// Loop possible derivatives.
1592
for (unsigned int r = 0; r < num_derivatives; r++)
1593
{
1594
// Resetting dmats values to compute next derivative.
1595
for (unsigned int t = 0; t < 4; t++)
1596
{
1597
for (unsigned int u = 0; u < 4; u++)
1598
{
1599
dmats[t][u] = 0.00000000;
1600
if (t == u)
1601
{
1602
dmats[t][u] = 1.00000000;
1603
}
1604
1605
}// end loop over 'u'
1606
}// end loop over 't'
1607
1608
// Looping derivative order to generate dmats.
1609
for (unsigned int s = 0; s < n; s++)
1610
{
1611
// Updating dmats_old with new values and resetting dmats.
1612
for (unsigned int t = 0; t < 4; t++)
1613
{
1614
for (unsigned int u = 0; u < 4; u++)
1615
{
1616
dmats_old[t][u] = dmats[t][u];
1617
dmats[t][u] = 0.00000000;
1618
}// end loop over 'u'
1619
}// end loop over 't'
1620
1621
// Update dmats using an inner product.
1622
if (combinations[r][s] == 0)
1623
{
1624
for (unsigned int t = 0; t < 4; t++)
1625
{
1626
for (unsigned int u = 0; u < 4; u++)
1627
{
1628
for (unsigned int tu = 0; tu < 4; tu++)
1629
{
1630
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1631
}// end loop over 'tu'
1632
}// end loop over 'u'
1633
}// end loop over 't'
1634
}
1635
1636
if (combinations[r][s] == 1)
1637
{
1638
for (unsigned int t = 0; t < 4; t++)
1639
{
1640
for (unsigned int u = 0; u < 4; u++)
1641
{
1642
for (unsigned int tu = 0; tu < 4; tu++)
1643
{
1644
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1645
}// end loop over 'tu'
1646
}// end loop over 'u'
1647
}// end loop over 't'
1648
}
1649
1650
if (combinations[r][s] == 2)
1651
{
1652
for (unsigned int t = 0; t < 4; t++)
1653
{
1654
for (unsigned int u = 0; u < 4; u++)
1655
{
1656
for (unsigned int tu = 0; tu < 4; tu++)
1657
{
1658
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1659
}// end loop over 'tu'
1660
}// end loop over 'u'
1661
}// end loop over 't'
1662
}
1663
1664
}// end loop over 's'
1665
for (unsigned int s = 0; s < 4; s++)
1666
{
1667
for (unsigned int t = 0; t < 4; t++)
1668
{
1669
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1670
}// end loop over 't'
1671
}// end loop over 's'
1672
}// end loop over 'r'
1673
1674
// Transform derivatives back to physical element
1675
for (unsigned int r = 0; r < num_derivatives; r++)
1676
{
1677
for (unsigned int s = 0; s < num_derivatives; s++)
1678
{
1679
values[r] += transform[r][s]*derivatives[s];
1680
}// end loop over 's'
1681
}// end loop over 'r'
1682
1683
// Delete pointer to array of derivatives on FIAT element
1684
delete [] derivatives;
1685
1686
// Delete pointer to array of combinations of derivatives and transform
1687
for (unsigned int r = 0; r < num_derivatives; r++)
1688
{
1689
delete [] combinations[r];
1690
}// end loop over 'r'
1691
delete [] combinations;
1692
for (unsigned int r = 0; r < num_derivatives; r++)
1693
{
1694
delete [] transform[r];
1695
}// end loop over 'r'
1696
delete [] transform;
1697
break;
1698
}
1699
}
1700
963
1701
}
964
1702
965
1703
/// Evaluate order n derivatives of all basis functions at given point in cell
1217
1955
values[0] = 0.00000000;
1218
1956
values[1] = 0.00000000;
1219
1957
values[2] = 0.00000000;
1220
if (0 <= i && i <= 3)
1221
{
1222
// Map degree of freedom to element degree of freedom
1223
const unsigned int dof = i;
1224
1225
// Array of basisvalues.
1226
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
1227
1228
// Declare helper variables.
1229
unsigned int rr = 0;
1230
unsigned int ss = 0;
1231
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
1232
1233
// Compute basisvalues.
1234
basisvalues[0] = 1.00000000;
1235
basisvalues[1] = tmp0;
1236
for (unsigned int r = 0; r < 1; r++)
1237
{
1238
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
1239
ss = r*(r + 1)*(r + 2)/6;
1240
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
1241
}// end loop over 'r'
1242
for (unsigned int r = 0; r < 1; r++)
1243
{
1244
for (unsigned int s = 0; s < 1 - r; s++)
1245
{
1246
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
1247
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
1248
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
1249
}// end loop over 's'
1250
}// end loop over 'r'
1251
for (unsigned int r = 0; r < 2; r++)
1252
{
1253
for (unsigned int s = 0; s < 2 - r; s++)
1254
{
1255
for (unsigned int t = 0; t < 2 - r - s; t++)
1256
{
1257
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
1258
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
1259
}// end loop over 't'
1260
}// end loop over 's'
1261
}// end loop over 'r'
1262
1263
// Table(s) of coefficients.
1264
static const double coefficients0[4][4] = \
1265
{{0.28867513, -0.18257419, -0.10540926, -0.07453560},
1266
{0.28867513, 0.18257419, -0.10540926, -0.07453560},
1267
{0.28867513, 0.00000000, 0.21081851, -0.07453560},
1268
{0.28867513, 0.00000000, 0.00000000, 0.22360680}};
1269
1270
// Compute value(s).
1271
for (unsigned int r = 0; r < 4; r++)
1272
{
1273
values[0] += coefficients0[dof][r]*basisvalues[r];
1274
}// end loop over 'r'
1275
}
1276
1277
if (4 <= i && i <= 7)
1278
{
1279
// Map degree of freedom to element degree of freedom
1280
const unsigned int dof = i - 4;
1281
1282
// Array of basisvalues.
1283
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
1284
1285
// Declare helper variables.
1286
unsigned int rr = 0;
1287
unsigned int ss = 0;
1288
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
1289
1290
// Compute basisvalues.
1291
basisvalues[0] = 1.00000000;
1292
basisvalues[1] = tmp0;
1293
for (unsigned int r = 0; r < 1; r++)
1294
{
1295
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
1296
ss = r*(r + 1)*(r + 2)/6;
1297
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
1298
}// end loop over 'r'
1299
for (unsigned int r = 0; r < 1; r++)
1300
{
1301
for (unsigned int s = 0; s < 1 - r; s++)
1302
{
1303
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
1304
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
1305
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
1306
}// end loop over 's'
1307
}// end loop over 'r'
1308
for (unsigned int r = 0; r < 2; r++)
1309
{
1310
for (unsigned int s = 0; s < 2 - r; s++)
1311
{
1312
for (unsigned int t = 0; t < 2 - r - s; t++)
1313
{
1314
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
1315
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
1316
}// end loop over 't'
1317
}// end loop over 's'
1318
}// end loop over 'r'
1319
1320
// Table(s) of coefficients.
1321
static const double coefficients0[4][4] = \
1322
{{0.28867513, -0.18257419, -0.10540926, -0.07453560},
1323
{0.28867513, 0.18257419, -0.10540926, -0.07453560},
1324
{0.28867513, 0.00000000, 0.21081851, -0.07453560},
1325
{0.28867513, 0.00000000, 0.00000000, 0.22360680}};
1326
1327
// Compute value(s).
1328
for (unsigned int r = 0; r < 4; r++)
1329
{
1330
values[1] += coefficients0[dof][r]*basisvalues[r];
1331
}// end loop over 'r'
1332
}
1333
1334
if (8 <= i && i <= 11)
1335
{
1336
// Map degree of freedom to element degree of freedom
1337
const unsigned int dof = i - 8;
1338
1339
// Array of basisvalues.
1340
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
1341
1342
// Declare helper variables.
1343
unsigned int rr = 0;
1344
unsigned int ss = 0;
1345
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
1346
1347
// Compute basisvalues.
1348
basisvalues[0] = 1.00000000;
1349
basisvalues[1] = tmp0;
1350
for (unsigned int r = 0; r < 1; r++)
1351
{
1352
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
1353
ss = r*(r + 1)*(r + 2)/6;
1354
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
1355
}// end loop over 'r'
1356
for (unsigned int r = 0; r < 1; r++)
1357
{
1358
for (unsigned int s = 0; s < 1 - r; s++)
1359
{
1360
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
1361
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
1362
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
1363
}// end loop over 's'
1364
}// end loop over 'r'
1365
for (unsigned int r = 0; r < 2; r++)
1366
{
1367
for (unsigned int s = 0; s < 2 - r; s++)
1368
{
1369
for (unsigned int t = 0; t < 2 - r - s; t++)
1370
{
1371
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
1372
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
1373
}// end loop over 't'
1374
}// end loop over 's'
1375
}// end loop over 'r'
1376
1377
// Table(s) of coefficients.
1378
static const double coefficients0[4][4] = \
1379
{{0.28867513, -0.18257419, -0.10540926, -0.07453560},
1380
{0.28867513, 0.18257419, -0.10540926, -0.07453560},
1381
{0.28867513, 0.00000000, 0.21081851, -0.07453560},
1382
{0.28867513, 0.00000000, 0.00000000, 0.22360680}};
1383
1384
// Compute value(s).
1385
for (unsigned int r = 0; r < 4; r++)
1386
{
1387
values[2] += coefficients0[dof][r]*basisvalues[r];
1388
}// end loop over 'r'
1958
switch (i)
1959
{
1960
case 0:
1961
{
1962
1963
// Array of basisvalues.
1964
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
1965
1966
// Declare helper variables.
1967
unsigned int rr = 0;
1968
unsigned int ss = 0;
1969
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
1970
1971
// Compute basisvalues.
1972
basisvalues[0] = 1.00000000;
1973
basisvalues[1] = tmp0;
1974
for (unsigned int r = 0; r < 1; r++)
1975
{
1976
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
1977
ss = r*(r + 1)*(r + 2)/6;
1978
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
1979
}// end loop over 'r'
1980
for (unsigned int r = 0; r < 1; r++)
1981
{
1982
for (unsigned int s = 0; s < 1 - r; s++)
1983
{
1984
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
1985
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
1986
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
1987
}// end loop over 's'
1988
}// end loop over 'r'
1989
for (unsigned int r = 0; r < 2; r++)
1990
{
1991
for (unsigned int s = 0; s < 2 - r; s++)
1992
{
1993
for (unsigned int t = 0; t < 2 - r - s; t++)
1994
{
1995
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
1996
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
1997
}// end loop over 't'
1998
}// end loop over 's'
1999
}// end loop over 'r'
2000
2001
// Table(s) of coefficients.
2002
static const double coefficients0[4] = \
2003
{0.28867513, -0.18257419, -0.10540926, -0.07453560};
2004
2005
// Compute value(s).
2006
for (unsigned int r = 0; r < 4; r++)
2007
{
2008
values[0] += coefficients0[r]*basisvalues[r];
2009
}// end loop over 'r'
2010
break;
2011
}
2012
case 1:
2013
{
2014
2015
// Array of basisvalues.
2016
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2017
2018
// Declare helper variables.
2019
unsigned int rr = 0;
2020
unsigned int ss = 0;
2021
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2022
2023
// Compute basisvalues.
2024
basisvalues[0] = 1.00000000;
2025
basisvalues[1] = tmp0;
2026
for (unsigned int r = 0; r < 1; r++)
2027
{
2028
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2029
ss = r*(r + 1)*(r + 2)/6;
2030
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2031
}// end loop over 'r'
2032
for (unsigned int r = 0; r < 1; r++)
2033
{
2034
for (unsigned int s = 0; s < 1 - r; s++)
2035
{
2036
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2037
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2038
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2039
}// end loop over 's'
2040
}// end loop over 'r'
2041
for (unsigned int r = 0; r < 2; r++)
2042
{
2043
for (unsigned int s = 0; s < 2 - r; s++)
2044
{
2045
for (unsigned int t = 0; t < 2 - r - s; t++)
2046
{
2047
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2048
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2049
}// end loop over 't'
2050
}// end loop over 's'
2051
}// end loop over 'r'
2052
2053
// Table(s) of coefficients.
2054
static const double coefficients0[4] = \
2055
{0.28867513, 0.18257419, -0.10540926, -0.07453560};
2056
2057
// Compute value(s).
2058
for (unsigned int r = 0; r < 4; r++)
2059
{
2060
values[0] += coefficients0[r]*basisvalues[r];
2061
}// end loop over 'r'
2062
break;
2063
}
2064
case 2:
2065
{
2066
2067
// Array of basisvalues.
2068
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2069
2070
// Declare helper variables.
2071
unsigned int rr = 0;
2072
unsigned int ss = 0;
2073
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2074
2075
// Compute basisvalues.
2076
basisvalues[0] = 1.00000000;
2077
basisvalues[1] = tmp0;
2078
for (unsigned int r = 0; r < 1; r++)
2079
{
2080
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2081
ss = r*(r + 1)*(r + 2)/6;
2082
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2083
}// end loop over 'r'
2084
for (unsigned int r = 0; r < 1; r++)
2085
{
2086
for (unsigned int s = 0; s < 1 - r; s++)
2087
{
2088
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2089
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2090
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2091
}// end loop over 's'
2092
}// end loop over 'r'
2093
for (unsigned int r = 0; r < 2; r++)
2094
{
2095
for (unsigned int s = 0; s < 2 - r; s++)
2096
{
2097
for (unsigned int t = 0; t < 2 - r - s; t++)
2098
{
2099
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2100
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2101
}// end loop over 't'
2102
}// end loop over 's'
2103
}// end loop over 'r'
2104
2105
// Table(s) of coefficients.
2106
static const double coefficients0[4] = \
2107
{0.28867513, 0.00000000, 0.21081851, -0.07453560};
2108
2109
// Compute value(s).
2110
for (unsigned int r = 0; r < 4; r++)
2111
{
2112
values[0] += coefficients0[r]*basisvalues[r];
2113
}// end loop over 'r'
2114
break;
2115
}
2116
case 3:
2117
{
2118
2119
// Array of basisvalues.
2120
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2121
2122
// Declare helper variables.
2123
unsigned int rr = 0;
2124
unsigned int ss = 0;
2125
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2126
2127
// Compute basisvalues.
2128
basisvalues[0] = 1.00000000;
2129
basisvalues[1] = tmp0;
2130
for (unsigned int r = 0; r < 1; r++)
2131
{
2132
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2133
ss = r*(r + 1)*(r + 2)/6;
2134
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2135
}// end loop over 'r'
2136
for (unsigned int r = 0; r < 1; r++)
2137
{
2138
for (unsigned int s = 0; s < 1 - r; s++)
2139
{
2140
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2141
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2142
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2143
}// end loop over 's'
2144
}// end loop over 'r'
2145
for (unsigned int r = 0; r < 2; r++)
2146
{
2147
for (unsigned int s = 0; s < 2 - r; s++)
2148
{
2149
for (unsigned int t = 0; t < 2 - r - s; t++)
2150
{
2151
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2152
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2153
}// end loop over 't'
2154
}// end loop over 's'
2155
}// end loop over 'r'
2156
2157
// Table(s) of coefficients.
2158
static const double coefficients0[4] = \
2159
{0.28867513, 0.00000000, 0.00000000, 0.22360680};
2160
2161
// Compute value(s).
2162
for (unsigned int r = 0; r < 4; r++)
2163
{
2164
values[0] += coefficients0[r]*basisvalues[r];
2165
}// end loop over 'r'
2166
break;
2167
}
2168
case 4:
2169
{
2170
2171
// Array of basisvalues.
2172
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2173
2174
// Declare helper variables.
2175
unsigned int rr = 0;
2176
unsigned int ss = 0;
2177
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2178
2179
// Compute basisvalues.
2180
basisvalues[0] = 1.00000000;
2181
basisvalues[1] = tmp0;
2182
for (unsigned int r = 0; r < 1; r++)
2183
{
2184
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2185
ss = r*(r + 1)*(r + 2)/6;
2186
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2187
}// end loop over 'r'
2188
for (unsigned int r = 0; r < 1; r++)
2189
{
2190
for (unsigned int s = 0; s < 1 - r; s++)
2191
{
2192
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2193
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2194
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2195
}// end loop over 's'
2196
}// end loop over 'r'
2197
for (unsigned int r = 0; r < 2; r++)
2198
{
2199
for (unsigned int s = 0; s < 2 - r; s++)
2200
{
2201
for (unsigned int t = 0; t < 2 - r - s; t++)
2202
{
2203
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2204
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2205
}// end loop over 't'
2206
}// end loop over 's'
2207
}// end loop over 'r'
2208
2209
// Table(s) of coefficients.
2210
static const double coefficients0[4] = \
2211
{0.28867513, -0.18257419, -0.10540926, -0.07453560};
2212
2213
// Compute value(s).
2214
for (unsigned int r = 0; r < 4; r++)
2215
{
2216
values[1] += coefficients0[r]*basisvalues[r];
2217
}// end loop over 'r'
2218
break;
2219
}
2220
case 5:
2221
{
2222
2223
// Array of basisvalues.
2224
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2225
2226
// Declare helper variables.
2227
unsigned int rr = 0;
2228
unsigned int ss = 0;
2229
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2230
2231
// Compute basisvalues.
2232
basisvalues[0] = 1.00000000;
2233
basisvalues[1] = tmp0;
2234
for (unsigned int r = 0; r < 1; r++)
2235
{
2236
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2237
ss = r*(r + 1)*(r + 2)/6;
2238
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2239
}// end loop over 'r'
2240
for (unsigned int r = 0; r < 1; r++)
2241
{
2242
for (unsigned int s = 0; s < 1 - r; s++)
2243
{
2244
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2245
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2246
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2247
}// end loop over 's'
2248
}// end loop over 'r'
2249
for (unsigned int r = 0; r < 2; r++)
2250
{
2251
for (unsigned int s = 0; s < 2 - r; s++)
2252
{
2253
for (unsigned int t = 0; t < 2 - r - s; t++)
2254
{
2255
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2256
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2257
}// end loop over 't'
2258
}// end loop over 's'
2259
}// end loop over 'r'
2260
2261
// Table(s) of coefficients.
2262
static const double coefficients0[4] = \
2263
{0.28867513, 0.18257419, -0.10540926, -0.07453560};
2264
2265
// Compute value(s).
2266
for (unsigned int r = 0; r < 4; r++)
2267
{
2268
values[1] += coefficients0[r]*basisvalues[r];
2269
}// end loop over 'r'
2270
break;
2271
}
2272
case 6:
2273
{
2274
2275
// Array of basisvalues.
2276
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2277
2278
// Declare helper variables.
2279
unsigned int rr = 0;
2280
unsigned int ss = 0;
2281
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2282
2283
// Compute basisvalues.
2284
basisvalues[0] = 1.00000000;
2285
basisvalues[1] = tmp0;
2286
for (unsigned int r = 0; r < 1; r++)
2287
{
2288
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2289
ss = r*(r + 1)*(r + 2)/6;
2290
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2291
}// end loop over 'r'
2292
for (unsigned int r = 0; r < 1; r++)
2293
{
2294
for (unsigned int s = 0; s < 1 - r; s++)
2295
{
2296
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2297
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2298
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2299
}// end loop over 's'
2300
}// end loop over 'r'
2301
for (unsigned int r = 0; r < 2; r++)
2302
{
2303
for (unsigned int s = 0; s < 2 - r; s++)
2304
{
2305
for (unsigned int t = 0; t < 2 - r - s; t++)
2306
{
2307
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2308
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2309
}// end loop over 't'
2310
}// end loop over 's'
2311
}// end loop over 'r'
2312
2313
// Table(s) of coefficients.
2314
static const double coefficients0[4] = \
2315
{0.28867513, 0.00000000, 0.21081851, -0.07453560};
2316
2317
// Compute value(s).
2318
for (unsigned int r = 0; r < 4; r++)
2319
{
2320
values[1] += coefficients0[r]*basisvalues[r];
2321
}// end loop over 'r'
2322
break;
2323
}
2324
case 7:
2325
{
2326
2327
// Array of basisvalues.
2328
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2329
2330
// Declare helper variables.
2331
unsigned int rr = 0;
2332
unsigned int ss = 0;
2333
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2334
2335
// Compute basisvalues.
2336
basisvalues[0] = 1.00000000;
2337
basisvalues[1] = tmp0;
2338
for (unsigned int r = 0; r < 1; r++)
2339
{
2340
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2341
ss = r*(r + 1)*(r + 2)/6;
2342
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2343
}// end loop over 'r'
2344
for (unsigned int r = 0; r < 1; r++)
2345
{
2346
for (unsigned int s = 0; s < 1 - r; s++)
2347
{
2348
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2349
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2350
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2351
}// end loop over 's'
2352
}// end loop over 'r'
2353
for (unsigned int r = 0; r < 2; r++)
2354
{
2355
for (unsigned int s = 0; s < 2 - r; s++)
2356
{
2357
for (unsigned int t = 0; t < 2 - r - s; t++)
2358
{
2359
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2360
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2361
}// end loop over 't'
2362
}// end loop over 's'
2363
}// end loop over 'r'
2364
2365
// Table(s) of coefficients.
2366
static const double coefficients0[4] = \
2367
{0.28867513, 0.00000000, 0.00000000, 0.22360680};
2368
2369
// Compute value(s).
2370
for (unsigned int r = 0; r < 4; r++)
2371
{
2372
values[1] += coefficients0[r]*basisvalues[r];
2373
}// end loop over 'r'
2374
break;
2375
}
2376
case 8:
2377
{
2378
2379
// Array of basisvalues.
2380
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2381
2382
// Declare helper variables.
2383
unsigned int rr = 0;
2384
unsigned int ss = 0;
2385
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2386
2387
// Compute basisvalues.
2388
basisvalues[0] = 1.00000000;
2389
basisvalues[1] = tmp0;
2390
for (unsigned int r = 0; r < 1; r++)
2391
{
2392
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2393
ss = r*(r + 1)*(r + 2)/6;
2394
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2395
}// end loop over 'r'
2396
for (unsigned int r = 0; r < 1; r++)
2397
{
2398
for (unsigned int s = 0; s < 1 - r; s++)
2399
{
2400
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2401
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2402
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2403
}// end loop over 's'
2404
}// end loop over 'r'
2405
for (unsigned int r = 0; r < 2; r++)
2406
{
2407
for (unsigned int s = 0; s < 2 - r; s++)
2408
{
2409
for (unsigned int t = 0; t < 2 - r - s; t++)
2410
{
2411
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2412
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2413
}// end loop over 't'
2414
}// end loop over 's'
2415
}// end loop over 'r'
2416
2417
// Table(s) of coefficients.
2418
static const double coefficients0[4] = \
2419
{0.28867513, -0.18257419, -0.10540926, -0.07453560};
2420
2421
// Compute value(s).
2422
for (unsigned int r = 0; r < 4; r++)
2423
{
2424
values[2] += coefficients0[r]*basisvalues[r];
2425
}// end loop over 'r'
2426
break;
2427
}
2428
case 9:
2429
{
2430
2431
// Array of basisvalues.
2432
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2433
2434
// Declare helper variables.
2435
unsigned int rr = 0;
2436
unsigned int ss = 0;
2437
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2438
2439
// Compute basisvalues.
2440
basisvalues[0] = 1.00000000;
2441
basisvalues[1] = tmp0;
2442
for (unsigned int r = 0; r < 1; r++)
2443
{
2444
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2445
ss = r*(r + 1)*(r + 2)/6;
2446
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2447
}// end loop over 'r'
2448
for (unsigned int r = 0; r < 1; r++)
2449
{
2450
for (unsigned int s = 0; s < 1 - r; s++)
2451
{
2452
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2453
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2454
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2455
}// end loop over 's'
2456
}// end loop over 'r'
2457
for (unsigned int r = 0; r < 2; r++)
2458
{
2459
for (unsigned int s = 0; s < 2 - r; s++)
2460
{
2461
for (unsigned int t = 0; t < 2 - r - s; t++)
2462
{
2463
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2464
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2465
}// end loop over 't'
2466
}// end loop over 's'
2467
}// end loop over 'r'
2468
2469
// Table(s) of coefficients.
2470
static const double coefficients0[4] = \
2471
{0.28867513, 0.18257419, -0.10540926, -0.07453560};
2472
2473
// Compute value(s).
2474
for (unsigned int r = 0; r < 4; r++)
2475
{
2476
values[2] += coefficients0[r]*basisvalues[r];
2477
}// end loop over 'r'
2478
break;
2479
}
2480
case 10:
2481
{
2482
2483
// Array of basisvalues.
2484
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2485
2486
// Declare helper variables.
2487
unsigned int rr = 0;
2488
unsigned int ss = 0;
2489
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2490
2491
// Compute basisvalues.
2492
basisvalues[0] = 1.00000000;
2493
basisvalues[1] = tmp0;
2494
for (unsigned int r = 0; r < 1; r++)
2495
{
2496
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2497
ss = r*(r + 1)*(r + 2)/6;
2498
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2499
}// end loop over 'r'
2500
for (unsigned int r = 0; r < 1; r++)
2501
{
2502
for (unsigned int s = 0; s < 1 - r; s++)
2503
{
2504
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2505
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2506
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2507
}// end loop over 's'
2508
}// end loop over 'r'
2509
for (unsigned int r = 0; r < 2; r++)
2510
{
2511
for (unsigned int s = 0; s < 2 - r; s++)
2512
{
2513
for (unsigned int t = 0; t < 2 - r - s; t++)
2514
{
2515
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2516
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2517
}// end loop over 't'
2518
}// end loop over 's'
2519
}// end loop over 'r'
2520
2521
// Table(s) of coefficients.
2522
static const double coefficients0[4] = \
2523
{0.28867513, 0.00000000, 0.21081851, -0.07453560};
2524
2525
// Compute value(s).
2526
for (unsigned int r = 0; r < 4; r++)
2527
{
2528
values[2] += coefficients0[r]*basisvalues[r];
2529
}// end loop over 'r'
2530
break;
2531
}
2532
case 11:
2533
{
2534
2535
// Array of basisvalues.
2536
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2537
2538
// Declare helper variables.
2539
unsigned int rr = 0;
2540
unsigned int ss = 0;
2541
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2542
2543
// Compute basisvalues.
2544
basisvalues[0] = 1.00000000;
2545
basisvalues[1] = tmp0;
2546
for (unsigned int r = 0; r < 1; r++)
2547
{
2548
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2549
ss = r*(r + 1)*(r + 2)/6;
2550
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2551
}// end loop over 'r'
2552
for (unsigned int r = 0; r < 1; r++)
2553
{
2554
for (unsigned int s = 0; s < 1 - r; s++)
2555
{
2556
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2557
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2558
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2559
}// end loop over 's'
2560
}// end loop over 'r'
2561
for (unsigned int r = 0; r < 2; r++)
2562
{
2563
for (unsigned int s = 0; s < 2 - r; s++)
2564
{
2565
for (unsigned int t = 0; t < 2 - r - s; t++)
2566
{
2567
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2568
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2569
}// end loop over 't'
2570
}// end loop over 's'
2571
}// end loop over 'r'
2572
2573
// Table(s) of coefficients.
2574
static const double coefficients0[4] = \
2575
{0.28867513, 0.00000000, 0.00000000, 0.22360680};
2576
2577
// Compute value(s).
2578
for (unsigned int r = 0; r < 4; r++)
2579
{
2580
values[2] += coefficients0[r]*basisvalues[r];
2581
}// end loop over 'r'
2582
break;
2583
}
1389
2584
}
1390
2585
1391
2586
}
1530
2725
values[r] = 0.00000000;
1531
2726
}// end loop over 'r'
1532
2727
1533
if (0 <= i && i <= 3)
1534
{
1535
// Map degree of freedom to element degree of freedom
1536
const unsigned int dof = i;
1537
1538
// Array of basisvalues.
1539
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
1540
1541
// Declare helper variables.
1542
unsigned int rr = 0;
1543
unsigned int ss = 0;
1544
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
1545
1546
// Compute basisvalues.
1547
basisvalues[0] = 1.00000000;
1548
basisvalues[1] = tmp0;
1549
for (unsigned int r = 0; r < 1; r++)
1550
{
1551
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
1552
ss = r*(r + 1)*(r + 2)/6;
1553
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
1554
}// end loop over 'r'
1555
for (unsigned int r = 0; r < 1; r++)
1556
{
1557
for (unsigned int s = 0; s < 1 - r; s++)
1558
{
1559
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
1560
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
1561
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
1562
}// end loop over 's'
1563
}// end loop over 'r'
1564
for (unsigned int r = 0; r < 2; r++)
1565
{
1566
for (unsigned int s = 0; s < 2 - r; s++)
1567
{
1568
for (unsigned int t = 0; t < 2 - r - s; t++)
1569
{
1570
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
1571
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
1572
}// end loop over 't'
1573
}// end loop over 's'
1574
}// end loop over 'r'
1575
1576
// Table(s) of coefficients.
1577
static const double coefficients0[4][4] = \
1578
{{0.28867513, -0.18257419, -0.10540926, -0.07453560},
1579
{0.28867513, 0.18257419, -0.10540926, -0.07453560},
1580
{0.28867513, 0.00000000, 0.21081851, -0.07453560},
1581
{0.28867513, 0.00000000, 0.00000000, 0.22360680}};
1582
1583
// Tables of derivatives of the polynomial base (transpose).
1584
static const double dmats0[4][4] = \
1585
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1586
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
1587
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1588
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1589
1590
static const double dmats1[4][4] = \
1591
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1592
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1593
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
1594
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1595
1596
static const double dmats2[4][4] = \
1597
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1598
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1599
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
1600
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
1601
1602
// Compute reference derivatives.
1603
// Declare pointer to array of derivatives on FIAT element.
1604
double *derivatives = new double[num_derivatives];
1605
for (unsigned int r = 0; r < num_derivatives; r++)
1606
{
1607
derivatives[r] = 0.00000000;
1608
}// end loop over 'r'
1609
1610
// Declare derivative matrix (of polynomial basis).
1611
double dmats[4][4] = \
1612
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1613
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1614
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1615
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1616
1617
// Declare (auxiliary) derivative matrix (of polynomial basis).
1618
double dmats_old[4][4] = \
1619
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1620
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1621
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1622
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1623
1624
// Loop possible derivatives.
1625
for (unsigned int r = 0; r < num_derivatives; r++)
1626
{
1627
// Resetting dmats values to compute next derivative.
1628
for (unsigned int t = 0; t < 4; t++)
1629
{
1630
for (unsigned int u = 0; u < 4; u++)
1631
{
1632
dmats[t][u] = 0.00000000;
1633
if (t == u)
1634
{
1635
dmats[t][u] = 1.00000000;
1636
}
1637
1638
}// end loop over 'u'
1639
}// end loop over 't'
1640
1641
// Looping derivative order to generate dmats.
1642
for (unsigned int s = 0; s < n; s++)
1643
{
1644
// Updating dmats_old with new values and resetting dmats.
1645
for (unsigned int t = 0; t < 4; t++)
1646
{
1647
for (unsigned int u = 0; u < 4; u++)
1648
{
1649
dmats_old[t][u] = dmats[t][u];
1650
dmats[t][u] = 0.00000000;
1651
}// end loop over 'u'
1652
}// end loop over 't'
1653
1654
// Update dmats using an inner product.
1655
if (combinations[r][s] == 0)
1656
{
1657
for (unsigned int t = 0; t < 4; t++)
1658
{
1659
for (unsigned int u = 0; u < 4; u++)
1660
{
1661
for (unsigned int tu = 0; tu < 4; tu++)
1662
{
1663
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1664
}// end loop over 'tu'
1665
}// end loop over 'u'
1666
}// end loop over 't'
1667
}
1668
1669
if (combinations[r][s] == 1)
1670
{
1671
for (unsigned int t = 0; t < 4; t++)
1672
{
1673
for (unsigned int u = 0; u < 4; u++)
1674
{
1675
for (unsigned int tu = 0; tu < 4; tu++)
1676
{
1677
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1678
}// end loop over 'tu'
1679
}// end loop over 'u'
1680
}// end loop over 't'
1681
}
1682
1683
if (combinations[r][s] == 2)
1684
{
1685
for (unsigned int t = 0; t < 4; t++)
1686
{
1687
for (unsigned int u = 0; u < 4; u++)
1688
{
1689
for (unsigned int tu = 0; tu < 4; tu++)
1690
{
1691
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1692
}// end loop over 'tu'
1693
}// end loop over 'u'
1694
}// end loop over 't'
1695
}
1696
1697
}// end loop over 's'
1698
for (unsigned int s = 0; s < 4; s++)
1699
{
1700
for (unsigned int t = 0; t < 4; t++)
1701
{
1702
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
1703
}// end loop over 't'
1704
}// end loop over 's'
1705
}// end loop over 'r'
1706
1707
// Transform derivatives back to physical element
1708
for (unsigned int r = 0; r < num_derivatives; r++)
1709
{
1710
for (unsigned int s = 0; s < num_derivatives; s++)
1711
{
1712
values[r] += transform[r][s]*derivatives[s];
1713
}// end loop over 's'
1714
}// end loop over 'r'
1715
1716
// Delete pointer to array of derivatives on FIAT element
1717
delete [] derivatives;
1718
1719
// Delete pointer to array of combinations of derivatives and transform
1720
for (unsigned int r = 0; r < num_derivatives; r++)
1721
{
1722
delete [] combinations[r];
1723
}// end loop over 'r'
1724
delete [] combinations;
1725
for (unsigned int r = 0; r < num_derivatives; r++)
1726
{
1727
delete [] transform[r];
1728
}// end loop over 'r'
1729
delete [] transform;
1730
}
1731
1732
if (4 <= i && i <= 7)
1733
{
1734
// Map degree of freedom to element degree of freedom
1735
const unsigned int dof = i - 4;
1736
1737
// Array of basisvalues.
1738
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
1739
1740
// Declare helper variables.
1741
unsigned int rr = 0;
1742
unsigned int ss = 0;
1743
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
1744
1745
// Compute basisvalues.
1746
basisvalues[0] = 1.00000000;
1747
basisvalues[1] = tmp0;
1748
for (unsigned int r = 0; r < 1; r++)
1749
{
1750
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
1751
ss = r*(r + 1)*(r + 2)/6;
1752
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
1753
}// end loop over 'r'
1754
for (unsigned int r = 0; r < 1; r++)
1755
{
1756
for (unsigned int s = 0; s < 1 - r; s++)
1757
{
1758
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
1759
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
1760
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
1761
}// end loop over 's'
1762
}// end loop over 'r'
1763
for (unsigned int r = 0; r < 2; r++)
1764
{
1765
for (unsigned int s = 0; s < 2 - r; s++)
1766
{
1767
for (unsigned int t = 0; t < 2 - r - s; t++)
1768
{
1769
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
1770
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
1771
}// end loop over 't'
1772
}// end loop over 's'
1773
}// end loop over 'r'
1774
1775
// Table(s) of coefficients.
1776
static const double coefficients0[4][4] = \
1777
{{0.28867513, -0.18257419, -0.10540926, -0.07453560},
1778
{0.28867513, 0.18257419, -0.10540926, -0.07453560},
1779
{0.28867513, 0.00000000, 0.21081851, -0.07453560},
1780
{0.28867513, 0.00000000, 0.00000000, 0.22360680}};
1781
1782
// Tables of derivatives of the polynomial base (transpose).
1783
static const double dmats0[4][4] = \
1784
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1785
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
1786
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1787
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1788
1789
static const double dmats1[4][4] = \
1790
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1791
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1792
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
1793
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1794
1795
static const double dmats2[4][4] = \
1796
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1797
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1798
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
1799
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
1800
1801
// Compute reference derivatives.
1802
// Declare pointer to array of derivatives on FIAT element.
1803
double *derivatives = new double[num_derivatives];
1804
for (unsigned int r = 0; r < num_derivatives; r++)
1805
{
1806
derivatives[r] = 0.00000000;
1807
}// end loop over 'r'
1808
1809
// Declare derivative matrix (of polynomial basis).
1810
double dmats[4][4] = \
1811
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1812
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1813
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1814
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1815
1816
// Declare (auxiliary) derivative matrix (of polynomial basis).
1817
double dmats_old[4][4] = \
1818
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
1819
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
1820
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
1821
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
1822
1823
// Loop possible derivatives.
1824
for (unsigned int r = 0; r < num_derivatives; r++)
1825
{
1826
// Resetting dmats values to compute next derivative.
1827
for (unsigned int t = 0; t < 4; t++)
1828
{
1829
for (unsigned int u = 0; u < 4; u++)
1830
{
1831
dmats[t][u] = 0.00000000;
1832
if (t == u)
1833
{
1834
dmats[t][u] = 1.00000000;
1835
}
1836
1837
}// end loop over 'u'
1838
}// end loop over 't'
1839
1840
// Looping derivative order to generate dmats.
1841
for (unsigned int s = 0; s < n; s++)
1842
{
1843
// Updating dmats_old with new values and resetting dmats.
1844
for (unsigned int t = 0; t < 4; t++)
1845
{
1846
for (unsigned int u = 0; u < 4; u++)
1847
{
1848
dmats_old[t][u] = dmats[t][u];
1849
dmats[t][u] = 0.00000000;
1850
}// end loop over 'u'
1851
}// end loop over 't'
1852
1853
// Update dmats using an inner product.
1854
if (combinations[r][s] == 0)
1855
{
1856
for (unsigned int t = 0; t < 4; t++)
1857
{
1858
for (unsigned int u = 0; u < 4; u++)
1859
{
1860
for (unsigned int tu = 0; tu < 4; tu++)
1861
{
1862
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1863
}// end loop over 'tu'
1864
}// end loop over 'u'
1865
}// end loop over 't'
1866
}
1867
1868
if (combinations[r][s] == 1)
1869
{
1870
for (unsigned int t = 0; t < 4; t++)
1871
{
1872
for (unsigned int u = 0; u < 4; u++)
1873
{
1874
for (unsigned int tu = 0; tu < 4; tu++)
1875
{
1876
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1877
}// end loop over 'tu'
1878
}// end loop over 'u'
1879
}// end loop over 't'
1880
}
1881
1882
if (combinations[r][s] == 2)
1883
{
1884
for (unsigned int t = 0; t < 4; t++)
1885
{
1886
for (unsigned int u = 0; u < 4; u++)
1887
{
1888
for (unsigned int tu = 0; tu < 4; tu++)
1889
{
1890
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1891
}// end loop over 'tu'
1892
}// end loop over 'u'
1893
}// end loop over 't'
1894
}
1895
1896
}// end loop over 's'
1897
for (unsigned int s = 0; s < 4; s++)
1898
{
1899
for (unsigned int t = 0; t < 4; t++)
1900
{
1901
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
1902
}// end loop over 't'
1903
}// end loop over 's'
1904
}// end loop over 'r'
1905
1906
// Transform derivatives back to physical element
1907
for (unsigned int r = 0; r < num_derivatives; r++)
1908
{
1909
for (unsigned int s = 0; s < num_derivatives; s++)
1910
{
1911
values[num_derivatives + r] += transform[r][s]*derivatives[s];
1912
}// end loop over 's'
1913
}// end loop over 'r'
1914
1915
// Delete pointer to array of derivatives on FIAT element
1916
delete [] derivatives;
1917
1918
// Delete pointer to array of combinations of derivatives and transform
1919
for (unsigned int r = 0; r < num_derivatives; r++)
1920
{
1921
delete [] combinations[r];
1922
}// end loop over 'r'
1923
delete [] combinations;
1924
for (unsigned int r = 0; r < num_derivatives; r++)
1925
{
1926
delete [] transform[r];
1927
}// end loop over 'r'
1928
delete [] transform;
1929
}
1930
1931
if (8 <= i && i <= 11)
1932
{
1933
// Map degree of freedom to element degree of freedom
1934
const unsigned int dof = i - 8;
1935
1936
// Array of basisvalues.
1937
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
1938
1939
// Declare helper variables.
1940
unsigned int rr = 0;
1941
unsigned int ss = 0;
1942
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
1943
1944
// Compute basisvalues.
1945
basisvalues[0] = 1.00000000;
1946
basisvalues[1] = tmp0;
1947
for (unsigned int r = 0; r < 1; r++)
1948
{
1949
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
1950
ss = r*(r + 1)*(r + 2)/6;
1951
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
1952
}// end loop over 'r'
1953
for (unsigned int r = 0; r < 1; r++)
1954
{
1955
for (unsigned int s = 0; s < 1 - r; s++)
1956
{
1957
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
1958
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
1959
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
1960
}// end loop over 's'
1961
}// end loop over 'r'
1962
for (unsigned int r = 0; r < 2; r++)
1963
{
1964
for (unsigned int s = 0; s < 2 - r; s++)
1965
{
1966
for (unsigned int t = 0; t < 2 - r - s; t++)
1967
{
1968
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
1969
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
1970
}// end loop over 't'
1971
}// end loop over 's'
1972
}// end loop over 'r'
1973
1974
// Table(s) of coefficients.
1975
static const double coefficients0[4][4] = \
1976
{{0.28867513, -0.18257419, -0.10540926, -0.07453560},
1977
{0.28867513, 0.18257419, -0.10540926, -0.07453560},
1978
{0.28867513, 0.00000000, 0.21081851, -0.07453560},
1979
{0.28867513, 0.00000000, 0.00000000, 0.22360680}};
1980
1981
// Tables of derivatives of the polynomial base (transpose).
1982
static const double dmats0[4][4] = \
1983
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1984
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
1985
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1986
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1987
1988
static const double dmats1[4][4] = \
1989
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1990
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1991
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
1992
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
1993
1994
static const double dmats2[4][4] = \
1995
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
1996
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
1997
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
1998
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
1999
2000
// Compute reference derivatives.
2001
// Declare pointer to array of derivatives on FIAT element.
2002
double *derivatives = new double[num_derivatives];
2003
for (unsigned int r = 0; r < num_derivatives; r++)
2004
{
2005
derivatives[r] = 0.00000000;
2006
}// end loop over 'r'
2007
2008
// Declare derivative matrix (of polynomial basis).
2009
double dmats[4][4] = \
2010
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
2011
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
2012
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
2013
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2014
2015
// Declare (auxiliary) derivative matrix (of polynomial basis).
2016
double dmats_old[4][4] = \
2017
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
2018
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
2019
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
2020
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2021
2022
// Loop possible derivatives.
2023
for (unsigned int r = 0; r < num_derivatives; r++)
2024
{
2025
// Resetting dmats values to compute next derivative.
2026
for (unsigned int t = 0; t < 4; t++)
2027
{
2028
for (unsigned int u = 0; u < 4; u++)
2029
{
2030
dmats[t][u] = 0.00000000;
2031
if (t == u)
2032
{
2033
dmats[t][u] = 1.00000000;
2034
}
2035
2036
}// end loop over 'u'
2037
}// end loop over 't'
2038
2039
// Looping derivative order to generate dmats.
2040
for (unsigned int s = 0; s < n; s++)
2041
{
2042
// Updating dmats_old with new values and resetting dmats.
2043
for (unsigned int t = 0; t < 4; t++)
2044
{
2045
for (unsigned int u = 0; u < 4; u++)
2046
{
2047
dmats_old[t][u] = dmats[t][u];
2048
dmats[t][u] = 0.00000000;
2049
}// end loop over 'u'
2050
}// end loop over 't'
2051
2052
// Update dmats using an inner product.
2053
if (combinations[r][s] == 0)
2054
{
2055
for (unsigned int t = 0; t < 4; t++)
2056
{
2057
for (unsigned int u = 0; u < 4; u++)
2058
{
2059
for (unsigned int tu = 0; tu < 4; tu++)
2060
{
2061
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2062
}// end loop over 'tu'
2063
}// end loop over 'u'
2064
}// end loop over 't'
2065
}
2066
2067
if (combinations[r][s] == 1)
2068
{
2069
for (unsigned int t = 0; t < 4; t++)
2070
{
2071
for (unsigned int u = 0; u < 4; u++)
2072
{
2073
for (unsigned int tu = 0; tu < 4; tu++)
2074
{
2075
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2076
}// end loop over 'tu'
2077
}// end loop over 'u'
2078
}// end loop over 't'
2079
}
2080
2081
if (combinations[r][s] == 2)
2082
{
2083
for (unsigned int t = 0; t < 4; t++)
2084
{
2085
for (unsigned int u = 0; u < 4; u++)
2086
{
2087
for (unsigned int tu = 0; tu < 4; tu++)
2088
{
2089
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2090
}// end loop over 'tu'
2091
}// end loop over 'u'
2092
}// end loop over 't'
2093
}
2094
2095
}// end loop over 's'
2096
for (unsigned int s = 0; s < 4; s++)
2097
{
2098
for (unsigned int t = 0; t < 4; t++)
2099
{
2100
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
2101
}// end loop over 't'
2102
}// end loop over 's'
2103
}// end loop over 'r'
2104
2105
// Transform derivatives back to physical element
2106
for (unsigned int r = 0; r < num_derivatives; r++)
2107
{
2108
for (unsigned int s = 0; s < num_derivatives; s++)
2109
{
2110
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
2111
}// end loop over 's'
2112
}// end loop over 'r'
2113
2114
// Delete pointer to array of derivatives on FIAT element
2115
delete [] derivatives;
2116
2117
// Delete pointer to array of combinations of derivatives and transform
2118
for (unsigned int r = 0; r < num_derivatives; r++)
2119
{
2120
delete [] combinations[r];
2121
}// end loop over 'r'
2122
delete [] combinations;
2123
for (unsigned int r = 0; r < num_derivatives; r++)
2124
{
2125
delete [] transform[r];
2126
}// end loop over 'r'
2127
delete [] transform;
2728
switch (i)
2729
{
2730
case 0:
2731
{
2732
2733
// Array of basisvalues.
2734
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2735
2736
// Declare helper variables.
2737
unsigned int rr = 0;
2738
unsigned int ss = 0;
2739
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2740
2741
// Compute basisvalues.
2742
basisvalues[0] = 1.00000000;
2743
basisvalues[1] = tmp0;
2744
for (unsigned int r = 0; r < 1; r++)
2745
{
2746
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2747
ss = r*(r + 1)*(r + 2)/6;
2748
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2749
}// end loop over 'r'
2750
for (unsigned int r = 0; r < 1; r++)
2751
{
2752
for (unsigned int s = 0; s < 1 - r; s++)
2753
{
2754
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2755
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2756
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2757
}// end loop over 's'
2758
}// end loop over 'r'
2759
for (unsigned int r = 0; r < 2; r++)
2760
{
2761
for (unsigned int s = 0; s < 2 - r; s++)
2762
{
2763
for (unsigned int t = 0; t < 2 - r - s; t++)
2764
{
2765
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2766
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2767
}// end loop over 't'
2768
}// end loop over 's'
2769
}// end loop over 'r'
2770
2771
// Table(s) of coefficients.
2772
static const double coefficients0[4] = \
2773
{0.28867513, -0.18257419, -0.10540926, -0.07453560};
2774
2775
// Tables of derivatives of the polynomial base (transpose).
2776
static const double dmats0[4][4] = \
2777
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
2778
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
2779
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
2780
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2781
2782
static const double dmats1[4][4] = \
2783
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
2784
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
2785
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
2786
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2787
2788
static const double dmats2[4][4] = \
2789
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
2790
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
2791
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
2792
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
2793
2794
// Compute reference derivatives.
2795
// Declare pointer to array of derivatives on FIAT element.
2796
double *derivatives = new double[num_derivatives];
2797
for (unsigned int r = 0; r < num_derivatives; r++)
2798
{
2799
derivatives[r] = 0.00000000;
2800
}// end loop over 'r'
2801
2802
// Declare derivative matrix (of polynomial basis).
2803
double dmats[4][4] = \
2804
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
2805
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
2806
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
2807
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2808
2809
// Declare (auxiliary) derivative matrix (of polynomial basis).
2810
double dmats_old[4][4] = \
2811
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
2812
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
2813
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
2814
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2815
2816
// Loop possible derivatives.
2817
for (unsigned int r = 0; r < num_derivatives; r++)
2818
{
2819
// Resetting dmats values to compute next derivative.
2820
for (unsigned int t = 0; t < 4; t++)
2821
{
2822
for (unsigned int u = 0; u < 4; u++)
2823
{
2824
dmats[t][u] = 0.00000000;
2825
if (t == u)
2826
{
2827
dmats[t][u] = 1.00000000;
2828
}
2829
2830
}// end loop over 'u'
2831
}// end loop over 't'
2832
2833
// Looping derivative order to generate dmats.
2834
for (unsigned int s = 0; s < n; s++)
2835
{
2836
// Updating dmats_old with new values and resetting dmats.
2837
for (unsigned int t = 0; t < 4; t++)
2838
{
2839
for (unsigned int u = 0; u < 4; u++)
2840
{
2841
dmats_old[t][u] = dmats[t][u];
2842
dmats[t][u] = 0.00000000;
2843
}// end loop over 'u'
2844
}// end loop over 't'
2845
2846
// Update dmats using an inner product.
2847
if (combinations[r][s] == 0)
2848
{
2849
for (unsigned int t = 0; t < 4; t++)
2850
{
2851
for (unsigned int u = 0; u < 4; u++)
2852
{
2853
for (unsigned int tu = 0; tu < 4; tu++)
2854
{
2855
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2856
}// end loop over 'tu'
2857
}// end loop over 'u'
2858
}// end loop over 't'
2859
}
2860
2861
if (combinations[r][s] == 1)
2862
{
2863
for (unsigned int t = 0; t < 4; t++)
2864
{
2865
for (unsigned int u = 0; u < 4; u++)
2866
{
2867
for (unsigned int tu = 0; tu < 4; tu++)
2868
{
2869
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2870
}// end loop over 'tu'
2871
}// end loop over 'u'
2872
}// end loop over 't'
2873
}
2874
2875
if (combinations[r][s] == 2)
2876
{
2877
for (unsigned int t = 0; t < 4; t++)
2878
{
2879
for (unsigned int u = 0; u < 4; u++)
2880
{
2881
for (unsigned int tu = 0; tu < 4; tu++)
2882
{
2883
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2884
}// end loop over 'tu'
2885
}// end loop over 'u'
2886
}// end loop over 't'
2887
}
2888
2889
}// end loop over 's'
2890
for (unsigned int s = 0; s < 4; s++)
2891
{
2892
for (unsigned int t = 0; t < 4; t++)
2893
{
2894
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2895
}// end loop over 't'
2896
}// end loop over 's'
2897
}// end loop over 'r'
2898
2899
// Transform derivatives back to physical element
2900
for (unsigned int r = 0; r < num_derivatives; r++)
2901
{
2902
for (unsigned int s = 0; s < num_derivatives; s++)
2903
{
2904
values[r] += transform[r][s]*derivatives[s];
2905
}// end loop over 's'
2906
}// end loop over 'r'
2907
2908
// Delete pointer to array of derivatives on FIAT element
2909
delete [] derivatives;
2910
2911
// Delete pointer to array of combinations of derivatives and transform
2912
for (unsigned int r = 0; r < num_derivatives; r++)
2913
{
2914
delete [] combinations[r];
2915
}// end loop over 'r'
2916
delete [] combinations;
2917
for (unsigned int r = 0; r < num_derivatives; r++)
2918
{
2919
delete [] transform[r];
2920
}// end loop over 'r'
2921
delete [] transform;
2922
break;
2923
}
2924
case 1:
2925
{
2926
2927
// Array of basisvalues.
2928
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
2929
2930
// Declare helper variables.
2931
unsigned int rr = 0;
2932
unsigned int ss = 0;
2933
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
2934
2935
// Compute basisvalues.
2936
basisvalues[0] = 1.00000000;
2937
basisvalues[1] = tmp0;
2938
for (unsigned int r = 0; r < 1; r++)
2939
{
2940
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
2941
ss = r*(r + 1)*(r + 2)/6;
2942
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
2943
}// end loop over 'r'
2944
for (unsigned int r = 0; r < 1; r++)
2945
{
2946
for (unsigned int s = 0; s < 1 - r; s++)
2947
{
2948
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
2949
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
2950
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
2951
}// end loop over 's'
2952
}// end loop over 'r'
2953
for (unsigned int r = 0; r < 2; r++)
2954
{
2955
for (unsigned int s = 0; s < 2 - r; s++)
2956
{
2957
for (unsigned int t = 0; t < 2 - r - s; t++)
2958
{
2959
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
2960
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
2961
}// end loop over 't'
2962
}// end loop over 's'
2963
}// end loop over 'r'
2964
2965
// Table(s) of coefficients.
2966
static const double coefficients0[4] = \
2967
{0.28867513, 0.18257419, -0.10540926, -0.07453560};
2968
2969
// Tables of derivatives of the polynomial base (transpose).
2970
static const double dmats0[4][4] = \
2971
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
2972
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
2973
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
2974
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2975
2976
static const double dmats1[4][4] = \
2977
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
2978
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
2979
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
2980
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
2981
2982
static const double dmats2[4][4] = \
2983
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
2984
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
2985
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
2986
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
2987
2988
// Compute reference derivatives.
2989
// Declare pointer to array of derivatives on FIAT element.
2990
double *derivatives = new double[num_derivatives];
2991
for (unsigned int r = 0; r < num_derivatives; r++)
2992
{
2993
derivatives[r] = 0.00000000;
2994
}// end loop over 'r'
2995
2996
// Declare derivative matrix (of polynomial basis).
2997
double dmats[4][4] = \
2998
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
2999
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3000
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3001
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3002
3003
// Declare (auxiliary) derivative matrix (of polynomial basis).
3004
double dmats_old[4][4] = \
3005
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
3006
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3007
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3008
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3009
3010
// Loop possible derivatives.
3011
for (unsigned int r = 0; r < num_derivatives; r++)
3012
{
3013
// Resetting dmats values to compute next derivative.
3014
for (unsigned int t = 0; t < 4; t++)
3015
{
3016
for (unsigned int u = 0; u < 4; u++)
3017
{
3018
dmats[t][u] = 0.00000000;
3019
if (t == u)
3020
{
3021
dmats[t][u] = 1.00000000;
3022
}
3023
3024
}// end loop over 'u'
3025
}// end loop over 't'
3026
3027
// Looping derivative order to generate dmats.
3028
for (unsigned int s = 0; s < n; s++)
3029
{
3030
// Updating dmats_old with new values and resetting dmats.
3031
for (unsigned int t = 0; t < 4; t++)
3032
{
3033
for (unsigned int u = 0; u < 4; u++)
3034
{
3035
dmats_old[t][u] = dmats[t][u];
3036
dmats[t][u] = 0.00000000;
3037
}// end loop over 'u'
3038
}// end loop over 't'
3039
3040
// Update dmats using an inner product.
3041
if (combinations[r][s] == 0)
3042
{
3043
for (unsigned int t = 0; t < 4; t++)
3044
{
3045
for (unsigned int u = 0; u < 4; u++)
3046
{
3047
for (unsigned int tu = 0; tu < 4; tu++)
3048
{
3049
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3050
}// end loop over 'tu'
3051
}// end loop over 'u'
3052
}// end loop over 't'
3053
}
3054
3055
if (combinations[r][s] == 1)
3056
{
3057
for (unsigned int t = 0; t < 4; t++)
3058
{
3059
for (unsigned int u = 0; u < 4; u++)
3060
{
3061
for (unsigned int tu = 0; tu < 4; tu++)
3062
{
3063
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3064
}// end loop over 'tu'
3065
}// end loop over 'u'
3066
}// end loop over 't'
3067
}
3068
3069
if (combinations[r][s] == 2)
3070
{
3071
for (unsigned int t = 0; t < 4; t++)
3072
{
3073
for (unsigned int u = 0; u < 4; u++)
3074
{
3075
for (unsigned int tu = 0; tu < 4; tu++)
3076
{
3077
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3078
}// end loop over 'tu'
3079
}// end loop over 'u'
3080
}// end loop over 't'
3081
}
3082
3083
}// end loop over 's'
3084
for (unsigned int s = 0; s < 4; s++)
3085
{
3086
for (unsigned int t = 0; t < 4; t++)
3087
{
3088
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3089
}// end loop over 't'
3090
}// end loop over 's'
3091
}// end loop over 'r'
3092
3093
// Transform derivatives back to physical element
3094
for (unsigned int r = 0; r < num_derivatives; r++)
3095
{
3096
for (unsigned int s = 0; s < num_derivatives; s++)
3097
{
3098
values[r] += transform[r][s]*derivatives[s];
3099
}// end loop over 's'
3100
}// end loop over 'r'
3101
3102
// Delete pointer to array of derivatives on FIAT element
3103
delete [] derivatives;
3104
3105
// Delete pointer to array of combinations of derivatives and transform
3106
for (unsigned int r = 0; r < num_derivatives; r++)
3107
{
3108
delete [] combinations[r];
3109
}// end loop over 'r'
3110
delete [] combinations;
3111
for (unsigned int r = 0; r < num_derivatives; r++)
3112
{
3113
delete [] transform[r];
3114
}// end loop over 'r'
3115
delete [] transform;
3116
break;
3117
}
3118
case 2:
3119
{
3120
3121
// Array of basisvalues.
3122
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
3123
3124
// Declare helper variables.
3125
unsigned int rr = 0;
3126
unsigned int ss = 0;
3127
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
3128
3129
// Compute basisvalues.
3130
basisvalues[0] = 1.00000000;
3131
basisvalues[1] = tmp0;
3132
for (unsigned int r = 0; r < 1; r++)
3133
{
3134
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
3135
ss = r*(r + 1)*(r + 2)/6;
3136
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
3137
}// end loop over 'r'
3138
for (unsigned int r = 0; r < 1; r++)
3139
{
3140
for (unsigned int s = 0; s < 1 - r; s++)
3141
{
3142
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
3143
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
3144
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
3145
}// end loop over 's'
3146
}// end loop over 'r'
3147
for (unsigned int r = 0; r < 2; r++)
3148
{
3149
for (unsigned int s = 0; s < 2 - r; s++)
3150
{
3151
for (unsigned int t = 0; t < 2 - r - s; t++)
3152
{
3153
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
3154
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
3155
}// end loop over 't'
3156
}// end loop over 's'
3157
}// end loop over 'r'
3158
3159
// Table(s) of coefficients.
3160
static const double coefficients0[4] = \
3161
{0.28867513, 0.00000000, 0.21081851, -0.07453560};
3162
3163
// Tables of derivatives of the polynomial base (transpose).
3164
static const double dmats0[4][4] = \
3165
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3166
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
3167
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3168
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3169
3170
static const double dmats1[4][4] = \
3171
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3172
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
3173
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
3174
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3175
3176
static const double dmats2[4][4] = \
3177
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3178
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
3179
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
3180
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
3181
3182
// Compute reference derivatives.
3183
// Declare pointer to array of derivatives on FIAT element.
3184
double *derivatives = new double[num_derivatives];
3185
for (unsigned int r = 0; r < num_derivatives; r++)
3186
{
3187
derivatives[r] = 0.00000000;
3188
}// end loop over 'r'
3189
3190
// Declare derivative matrix (of polynomial basis).
3191
double dmats[4][4] = \
3192
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
3193
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3194
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3195
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3196
3197
// Declare (auxiliary) derivative matrix (of polynomial basis).
3198
double dmats_old[4][4] = \
3199
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
3200
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3201
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3202
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3203
3204
// Loop possible derivatives.
3205
for (unsigned int r = 0; r < num_derivatives; r++)
3206
{
3207
// Resetting dmats values to compute next derivative.
3208
for (unsigned int t = 0; t < 4; t++)
3209
{
3210
for (unsigned int u = 0; u < 4; u++)
3211
{
3212
dmats[t][u] = 0.00000000;
3213
if (t == u)
3214
{
3215
dmats[t][u] = 1.00000000;
3216
}
3217
3218
}// end loop over 'u'
3219
}// end loop over 't'
3220
3221
// Looping derivative order to generate dmats.
3222
for (unsigned int s = 0; s < n; s++)
3223
{
3224
// Updating dmats_old with new values and resetting dmats.
3225
for (unsigned int t = 0; t < 4; t++)
3226
{
3227
for (unsigned int u = 0; u < 4; u++)
3228
{
3229
dmats_old[t][u] = dmats[t][u];
3230
dmats[t][u] = 0.00000000;
3231
}// end loop over 'u'
3232
}// end loop over 't'
3233
3234
// Update dmats using an inner product.
3235
if (combinations[r][s] == 0)
3236
{
3237
for (unsigned int t = 0; t < 4; t++)
3238
{
3239
for (unsigned int u = 0; u < 4; u++)
3240
{
3241
for (unsigned int tu = 0; tu < 4; tu++)
3242
{
3243
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3244
}// end loop over 'tu'
3245
}// end loop over 'u'
3246
}// end loop over 't'
3247
}
3248
3249
if (combinations[r][s] == 1)
3250
{
3251
for (unsigned int t = 0; t < 4; t++)
3252
{
3253
for (unsigned int u = 0; u < 4; u++)
3254
{
3255
for (unsigned int tu = 0; tu < 4; tu++)
3256
{
3257
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3258
}// end loop over 'tu'
3259
}// end loop over 'u'
3260
}// end loop over 't'
3261
}
3262
3263
if (combinations[r][s] == 2)
3264
{
3265
for (unsigned int t = 0; t < 4; t++)
3266
{
3267
for (unsigned int u = 0; u < 4; u++)
3268
{
3269
for (unsigned int tu = 0; tu < 4; tu++)
3270
{
3271
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3272
}// end loop over 'tu'
3273
}// end loop over 'u'
3274
}// end loop over 't'
3275
}
3276
3277
}// end loop over 's'
3278
for (unsigned int s = 0; s < 4; s++)
3279
{
3280
for (unsigned int t = 0; t < 4; t++)
3281
{
3282
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3283
}// end loop over 't'
3284
}// end loop over 's'
3285
}// end loop over 'r'
3286
3287
// Transform derivatives back to physical element
3288
for (unsigned int r = 0; r < num_derivatives; r++)
3289
{
3290
for (unsigned int s = 0; s < num_derivatives; s++)
3291
{
3292
values[r] += transform[r][s]*derivatives[s];
3293
}// end loop over 's'
3294
}// end loop over 'r'
3295
3296
// Delete pointer to array of derivatives on FIAT element
3297
delete [] derivatives;
3298
3299
// Delete pointer to array of combinations of derivatives and transform
3300
for (unsigned int r = 0; r < num_derivatives; r++)
3301
{
3302
delete [] combinations[r];
3303
}// end loop over 'r'
3304
delete [] combinations;
3305
for (unsigned int r = 0; r < num_derivatives; r++)
3306
{
3307
delete [] transform[r];
3308
}// end loop over 'r'
3309
delete [] transform;
3310
break;
3311
}
3312
case 3:
3313
{
3314
3315
// Array of basisvalues.
3316
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
3317
3318
// Declare helper variables.
3319
unsigned int rr = 0;
3320
unsigned int ss = 0;
3321
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
3322
3323
// Compute basisvalues.
3324
basisvalues[0] = 1.00000000;
3325
basisvalues[1] = tmp0;
3326
for (unsigned int r = 0; r < 1; r++)
3327
{
3328
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
3329
ss = r*(r + 1)*(r + 2)/6;
3330
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
3331
}// end loop over 'r'
3332
for (unsigned int r = 0; r < 1; r++)
3333
{
3334
for (unsigned int s = 0; s < 1 - r; s++)
3335
{
3336
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
3337
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
3338
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
3339
}// end loop over 's'
3340
}// end loop over 'r'
3341
for (unsigned int r = 0; r < 2; r++)
3342
{
3343
for (unsigned int s = 0; s < 2 - r; s++)
3344
{
3345
for (unsigned int t = 0; t < 2 - r - s; t++)
3346
{
3347
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
3348
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
3349
}// end loop over 't'
3350
}// end loop over 's'
3351
}// end loop over 'r'
3352
3353
// Table(s) of coefficients.
3354
static const double coefficients0[4] = \
3355
{0.28867513, 0.00000000, 0.00000000, 0.22360680};
3356
3357
// Tables of derivatives of the polynomial base (transpose).
3358
static const double dmats0[4][4] = \
3359
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3360
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
3361
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3362
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3363
3364
static const double dmats1[4][4] = \
3365
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3366
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
3367
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
3368
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3369
3370
static const double dmats2[4][4] = \
3371
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3372
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
3373
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
3374
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
3375
3376
// Compute reference derivatives.
3377
// Declare pointer to array of derivatives on FIAT element.
3378
double *derivatives = new double[num_derivatives];
3379
for (unsigned int r = 0; r < num_derivatives; r++)
3380
{
3381
derivatives[r] = 0.00000000;
3382
}// end loop over 'r'
3383
3384
// Declare derivative matrix (of polynomial basis).
3385
double dmats[4][4] = \
3386
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
3387
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3388
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3389
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3390
3391
// Declare (auxiliary) derivative matrix (of polynomial basis).
3392
double dmats_old[4][4] = \
3393
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
3394
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3395
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3396
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3397
3398
// Loop possible derivatives.
3399
for (unsigned int r = 0; r < num_derivatives; r++)
3400
{
3401
// Resetting dmats values to compute next derivative.
3402
for (unsigned int t = 0; t < 4; t++)
3403
{
3404
for (unsigned int u = 0; u < 4; u++)
3405
{
3406
dmats[t][u] = 0.00000000;
3407
if (t == u)
3408
{
3409
dmats[t][u] = 1.00000000;
3410
}
3411
3412
}// end loop over 'u'
3413
}// end loop over 't'
3414
3415
// Looping derivative order to generate dmats.
3416
for (unsigned int s = 0; s < n; s++)
3417
{
3418
// Updating dmats_old with new values and resetting dmats.
3419
for (unsigned int t = 0; t < 4; t++)
3420
{
3421
for (unsigned int u = 0; u < 4; u++)
3422
{
3423
dmats_old[t][u] = dmats[t][u];
3424
dmats[t][u] = 0.00000000;
3425
}// end loop over 'u'
3426
}// end loop over 't'
3427
3428
// Update dmats using an inner product.
3429
if (combinations[r][s] == 0)
3430
{
3431
for (unsigned int t = 0; t < 4; t++)
3432
{
3433
for (unsigned int u = 0; u < 4; u++)
3434
{
3435
for (unsigned int tu = 0; tu < 4; tu++)
3436
{
3437
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3438
}// end loop over 'tu'
3439
}// end loop over 'u'
3440
}// end loop over 't'
3441
}
3442
3443
if (combinations[r][s] == 1)
3444
{
3445
for (unsigned int t = 0; t < 4; t++)
3446
{
3447
for (unsigned int u = 0; u < 4; u++)
3448
{
3449
for (unsigned int tu = 0; tu < 4; tu++)
3450
{
3451
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3452
}// end loop over 'tu'
3453
}// end loop over 'u'
3454
}// end loop over 't'
3455
}
3456
3457
if (combinations[r][s] == 2)
3458
{
3459
for (unsigned int t = 0; t < 4; t++)
3460
{
3461
for (unsigned int u = 0; u < 4; u++)
3462
{
3463
for (unsigned int tu = 0; tu < 4; tu++)
3464
{
3465
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3466
}// end loop over 'tu'
3467
}// end loop over 'u'
3468
}// end loop over 't'
3469
}
3470
3471
}// end loop over 's'
3472
for (unsigned int s = 0; s < 4; s++)
3473
{
3474
for (unsigned int t = 0; t < 4; t++)
3475
{
3476
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3477
}// end loop over 't'
3478
}// end loop over 's'
3479
}// end loop over 'r'
3480
3481
// Transform derivatives back to physical element
3482
for (unsigned int r = 0; r < num_derivatives; r++)
3483
{
3484
for (unsigned int s = 0; s < num_derivatives; s++)
3485
{
3486
values[r] += transform[r][s]*derivatives[s];
3487
}// end loop over 's'
3488
}// end loop over 'r'
3489
3490
// Delete pointer to array of derivatives on FIAT element
3491
delete [] derivatives;
3492
3493
// Delete pointer to array of combinations of derivatives and transform
3494
for (unsigned int r = 0; r < num_derivatives; r++)
3495
{
3496
delete [] combinations[r];
3497
}// end loop over 'r'
3498
delete [] combinations;
3499
for (unsigned int r = 0; r < num_derivatives; r++)
3500
{
3501
delete [] transform[r];
3502
}// end loop over 'r'
3503
delete [] transform;
3504
break;
3505
}
3506
case 4:
3507
{
3508
3509
// Array of basisvalues.
3510
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
3511
3512
// Declare helper variables.
3513
unsigned int rr = 0;
3514
unsigned int ss = 0;
3515
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
3516
3517
// Compute basisvalues.
3518
basisvalues[0] = 1.00000000;
3519
basisvalues[1] = tmp0;
3520
for (unsigned int r = 0; r < 1; r++)
3521
{
3522
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
3523
ss = r*(r + 1)*(r + 2)/6;
3524
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
3525
}// end loop over 'r'
3526
for (unsigned int r = 0; r < 1; r++)
3527
{
3528
for (unsigned int s = 0; s < 1 - r; s++)
3529
{
3530
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
3531
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
3532
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
3533
}// end loop over 's'
3534
}// end loop over 'r'
3535
for (unsigned int r = 0; r < 2; r++)
3536
{
3537
for (unsigned int s = 0; s < 2 - r; s++)
3538
{
3539
for (unsigned int t = 0; t < 2 - r - s; t++)
3540
{
3541
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
3542
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
3543
}// end loop over 't'
3544
}// end loop over 's'
3545
}// end loop over 'r'
3546
3547
// Table(s) of coefficients.
3548
static const double coefficients0[4] = \
3549
{0.28867513, -0.18257419, -0.10540926, -0.07453560};
3550
3551
// Tables of derivatives of the polynomial base (transpose).
3552
static const double dmats0[4][4] = \
3553
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3554
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
3555
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3556
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3557
3558
static const double dmats1[4][4] = \
3559
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3560
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
3561
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
3562
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3563
3564
static const double dmats2[4][4] = \
3565
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3566
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
3567
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
3568
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
3569
3570
// Compute reference derivatives.
3571
// Declare pointer to array of derivatives on FIAT element.
3572
double *derivatives = new double[num_derivatives];
3573
for (unsigned int r = 0; r < num_derivatives; r++)
3574
{
3575
derivatives[r] = 0.00000000;
3576
}// end loop over 'r'
3577
3578
// Declare derivative matrix (of polynomial basis).
3579
double dmats[4][4] = \
3580
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
3581
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3582
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3583
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3584
3585
// Declare (auxiliary) derivative matrix (of polynomial basis).
3586
double dmats_old[4][4] = \
3587
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
3588
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3589
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3590
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3591
3592
// Loop possible derivatives.
3593
for (unsigned int r = 0; r < num_derivatives; r++)
3594
{
3595
// Resetting dmats values to compute next derivative.
3596
for (unsigned int t = 0; t < 4; t++)
3597
{
3598
for (unsigned int u = 0; u < 4; u++)
3599
{
3600
dmats[t][u] = 0.00000000;
3601
if (t == u)
3602
{
3603
dmats[t][u] = 1.00000000;
3604
}
3605
3606
}// end loop over 'u'
3607
}// end loop over 't'
3608
3609
// Looping derivative order to generate dmats.
3610
for (unsigned int s = 0; s < n; s++)
3611
{
3612
// Updating dmats_old with new values and resetting dmats.
3613
for (unsigned int t = 0; t < 4; t++)
3614
{
3615
for (unsigned int u = 0; u < 4; u++)
3616
{
3617
dmats_old[t][u] = dmats[t][u];
3618
dmats[t][u] = 0.00000000;
3619
}// end loop over 'u'
3620
}// end loop over 't'
3621
3622
// Update dmats using an inner product.
3623
if (combinations[r][s] == 0)
3624
{
3625
for (unsigned int t = 0; t < 4; t++)
3626
{
3627
for (unsigned int u = 0; u < 4; u++)
3628
{
3629
for (unsigned int tu = 0; tu < 4; tu++)
3630
{
3631
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3632
}// end loop over 'tu'
3633
}// end loop over 'u'
3634
}// end loop over 't'
3635
}
3636
3637
if (combinations[r][s] == 1)
3638
{
3639
for (unsigned int t = 0; t < 4; t++)
3640
{
3641
for (unsigned int u = 0; u < 4; u++)
3642
{
3643
for (unsigned int tu = 0; tu < 4; tu++)
3644
{
3645
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3646
}// end loop over 'tu'
3647
}// end loop over 'u'
3648
}// end loop over 't'
3649
}
3650
3651
if (combinations[r][s] == 2)
3652
{
3653
for (unsigned int t = 0; t < 4; t++)
3654
{
3655
for (unsigned int u = 0; u < 4; u++)
3656
{
3657
for (unsigned int tu = 0; tu < 4; tu++)
3658
{
3659
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3660
}// end loop over 'tu'
3661
}// end loop over 'u'
3662
}// end loop over 't'
3663
}
3664
3665
}// end loop over 's'
3666
for (unsigned int s = 0; s < 4; s++)
3667
{
3668
for (unsigned int t = 0; t < 4; t++)
3669
{
3670
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3671
}// end loop over 't'
3672
}// end loop over 's'
3673
}// end loop over 'r'
3674
3675
// Transform derivatives back to physical element
3676
for (unsigned int r = 0; r < num_derivatives; r++)
3677
{
3678
for (unsigned int s = 0; s < num_derivatives; s++)
3679
{
3680
values[num_derivatives + r] += transform[r][s]*derivatives[s];
3681
}// end loop over 's'
3682
}// end loop over 'r'
3683
3684
// Delete pointer to array of derivatives on FIAT element
3685
delete [] derivatives;
3686
3687
// Delete pointer to array of combinations of derivatives and transform
3688
for (unsigned int r = 0; r < num_derivatives; r++)
3689
{
3690
delete [] combinations[r];
3691
}// end loop over 'r'
3692
delete [] combinations;
3693
for (unsigned int r = 0; r < num_derivatives; r++)
3694
{
3695
delete [] transform[r];
3696
}// end loop over 'r'
3697
delete [] transform;
3698
break;
3699
}
3700
case 5:
3701
{
3702
3703
// Array of basisvalues.
3704
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
3705
3706
// Declare helper variables.
3707
unsigned int rr = 0;
3708
unsigned int ss = 0;
3709
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
3710
3711
// Compute basisvalues.
3712
basisvalues[0] = 1.00000000;
3713
basisvalues[1] = tmp0;
3714
for (unsigned int r = 0; r < 1; r++)
3715
{
3716
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
3717
ss = r*(r + 1)*(r + 2)/6;
3718
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
3719
}// end loop over 'r'
3720
for (unsigned int r = 0; r < 1; r++)
3721
{
3722
for (unsigned int s = 0; s < 1 - r; s++)
3723
{
3724
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
3725
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
3726
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
3727
}// end loop over 's'
3728
}// end loop over 'r'
3729
for (unsigned int r = 0; r < 2; r++)
3730
{
3731
for (unsigned int s = 0; s < 2 - r; s++)
3732
{
3733
for (unsigned int t = 0; t < 2 - r - s; t++)
3734
{
3735
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
3736
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
3737
}// end loop over 't'
3738
}// end loop over 's'
3739
}// end loop over 'r'
3740
3741
// Table(s) of coefficients.
3742
static const double coefficients0[4] = \
3743
{0.28867513, 0.18257419, -0.10540926, -0.07453560};
3744
3745
// Tables of derivatives of the polynomial base (transpose).
3746
static const double dmats0[4][4] = \
3747
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3748
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
3749
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3750
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3751
3752
static const double dmats1[4][4] = \
3753
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3754
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
3755
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
3756
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3757
3758
static const double dmats2[4][4] = \
3759
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3760
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
3761
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
3762
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
3763
3764
// Compute reference derivatives.
3765
// Declare pointer to array of derivatives on FIAT element.
3766
double *derivatives = new double[num_derivatives];
3767
for (unsigned int r = 0; r < num_derivatives; r++)
3768
{
3769
derivatives[r] = 0.00000000;
3770
}// end loop over 'r'
3771
3772
// Declare derivative matrix (of polynomial basis).
3773
double dmats[4][4] = \
3774
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
3775
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3776
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3777
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3778
3779
// Declare (auxiliary) derivative matrix (of polynomial basis).
3780
double dmats_old[4][4] = \
3781
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
3782
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3783
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3784
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3785
3786
// Loop possible derivatives.
3787
for (unsigned int r = 0; r < num_derivatives; r++)
3788
{
3789
// Resetting dmats values to compute next derivative.
3790
for (unsigned int t = 0; t < 4; t++)
3791
{
3792
for (unsigned int u = 0; u < 4; u++)
3793
{
3794
dmats[t][u] = 0.00000000;
3795
if (t == u)
3796
{
3797
dmats[t][u] = 1.00000000;
3798
}
3799
3800
}// end loop over 'u'
3801
}// end loop over 't'
3802
3803
// Looping derivative order to generate dmats.
3804
for (unsigned int s = 0; s < n; s++)
3805
{
3806
// Updating dmats_old with new values and resetting dmats.
3807
for (unsigned int t = 0; t < 4; t++)
3808
{
3809
for (unsigned int u = 0; u < 4; u++)
3810
{
3811
dmats_old[t][u] = dmats[t][u];
3812
dmats[t][u] = 0.00000000;
3813
}// end loop over 'u'
3814
}// end loop over 't'
3815
3816
// Update dmats using an inner product.
3817
if (combinations[r][s] == 0)
3818
{
3819
for (unsigned int t = 0; t < 4; t++)
3820
{
3821
for (unsigned int u = 0; u < 4; u++)
3822
{
3823
for (unsigned int tu = 0; tu < 4; tu++)
3824
{
3825
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
3826
}// end loop over 'tu'
3827
}// end loop over 'u'
3828
}// end loop over 't'
3829
}
3830
3831
if (combinations[r][s] == 1)
3832
{
3833
for (unsigned int t = 0; t < 4; t++)
3834
{
3835
for (unsigned int u = 0; u < 4; u++)
3836
{
3837
for (unsigned int tu = 0; tu < 4; tu++)
3838
{
3839
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
3840
}// end loop over 'tu'
3841
}// end loop over 'u'
3842
}// end loop over 't'
3843
}
3844
3845
if (combinations[r][s] == 2)
3846
{
3847
for (unsigned int t = 0; t < 4; t++)
3848
{
3849
for (unsigned int u = 0; u < 4; u++)
3850
{
3851
for (unsigned int tu = 0; tu < 4; tu++)
3852
{
3853
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
3854
}// end loop over 'tu'
3855
}// end loop over 'u'
3856
}// end loop over 't'
3857
}
3858
3859
}// end loop over 's'
3860
for (unsigned int s = 0; s < 4; s++)
3861
{
3862
for (unsigned int t = 0; t < 4; t++)
3863
{
3864
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
3865
}// end loop over 't'
3866
}// end loop over 's'
3867
}// end loop over 'r'
3868
3869
// Transform derivatives back to physical element
3870
for (unsigned int r = 0; r < num_derivatives; r++)
3871
{
3872
for (unsigned int s = 0; s < num_derivatives; s++)
3873
{
3874
values[num_derivatives + r] += transform[r][s]*derivatives[s];
3875
}// end loop over 's'
3876
}// end loop over 'r'
3877
3878
// Delete pointer to array of derivatives on FIAT element
3879
delete [] derivatives;
3880
3881
// Delete pointer to array of combinations of derivatives and transform
3882
for (unsigned int r = 0; r < num_derivatives; r++)
3883
{
3884
delete [] combinations[r];
3885
}// end loop over 'r'
3886
delete [] combinations;
3887
for (unsigned int r = 0; r < num_derivatives; r++)
3888
{
3889
delete [] transform[r];
3890
}// end loop over 'r'
3891
delete [] transform;
3892
break;
3893
}
3894
case 6:
3895
{
3896
3897
// Array of basisvalues.
3898
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
3899
3900
// Declare helper variables.
3901
unsigned int rr = 0;
3902
unsigned int ss = 0;
3903
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
3904
3905
// Compute basisvalues.
3906
basisvalues[0] = 1.00000000;
3907
basisvalues[1] = tmp0;
3908
for (unsigned int r = 0; r < 1; r++)
3909
{
3910
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
3911
ss = r*(r + 1)*(r + 2)/6;
3912
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
3913
}// end loop over 'r'
3914
for (unsigned int r = 0; r < 1; r++)
3915
{
3916
for (unsigned int s = 0; s < 1 - r; s++)
3917
{
3918
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
3919
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
3920
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
3921
}// end loop over 's'
3922
}// end loop over 'r'
3923
for (unsigned int r = 0; r < 2; r++)
3924
{
3925
for (unsigned int s = 0; s < 2 - r; s++)
3926
{
3927
for (unsigned int t = 0; t < 2 - r - s; t++)
3928
{
3929
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
3930
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
3931
}// end loop over 't'
3932
}// end loop over 's'
3933
}// end loop over 'r'
3934
3935
// Table(s) of coefficients.
3936
static const double coefficients0[4] = \
3937
{0.28867513, 0.00000000, 0.21081851, -0.07453560};
3938
3939
// Tables of derivatives of the polynomial base (transpose).
3940
static const double dmats0[4][4] = \
3941
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3942
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
3943
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3944
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3945
3946
static const double dmats1[4][4] = \
3947
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3948
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
3949
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
3950
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
3951
3952
static const double dmats2[4][4] = \
3953
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
3954
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
3955
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
3956
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
3957
3958
// Compute reference derivatives.
3959
// Declare pointer to array of derivatives on FIAT element.
3960
double *derivatives = new double[num_derivatives];
3961
for (unsigned int r = 0; r < num_derivatives; r++)
3962
{
3963
derivatives[r] = 0.00000000;
3964
}// end loop over 'r'
3965
3966
// Declare derivative matrix (of polynomial basis).
3967
double dmats[4][4] = \
3968
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
3969
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3970
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3971
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3972
3973
// Declare (auxiliary) derivative matrix (of polynomial basis).
3974
double dmats_old[4][4] = \
3975
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
3976
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
3977
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
3978
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
3979
3980
// Loop possible derivatives.
3981
for (unsigned int r = 0; r < num_derivatives; r++)
3982
{
3983
// Resetting dmats values to compute next derivative.
3984
for (unsigned int t = 0; t < 4; t++)
3985
{
3986
for (unsigned int u = 0; u < 4; u++)
3987
{
3988
dmats[t][u] = 0.00000000;
3989
if (t == u)
3990
{
3991
dmats[t][u] = 1.00000000;
3992
}
3993
3994
}// end loop over 'u'
3995
}// end loop over 't'
3996
3997
// Looping derivative order to generate dmats.
3998
for (unsigned int s = 0; s < n; s++)
3999
{
4000
// Updating dmats_old with new values and resetting dmats.
4001
for (unsigned int t = 0; t < 4; t++)
4002
{
4003
for (unsigned int u = 0; u < 4; u++)
4004
{
4005
dmats_old[t][u] = dmats[t][u];
4006
dmats[t][u] = 0.00000000;
4007
}// end loop over 'u'
4008
}// end loop over 't'
4009
4010
// Update dmats using an inner product.
4011
if (combinations[r][s] == 0)
4012
{
4013
for (unsigned int t = 0; t < 4; t++)
4014
{
4015
for (unsigned int u = 0; u < 4; u++)
4016
{
4017
for (unsigned int tu = 0; tu < 4; tu++)
4018
{
4019
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
4020
}// end loop over 'tu'
4021
}// end loop over 'u'
4022
}// end loop over 't'
4023
}
4024
4025
if (combinations[r][s] == 1)
4026
{
4027
for (unsigned int t = 0; t < 4; t++)
4028
{
4029
for (unsigned int u = 0; u < 4; u++)
4030
{
4031
for (unsigned int tu = 0; tu < 4; tu++)
4032
{
4033
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
4034
}// end loop over 'tu'
4035
}// end loop over 'u'
4036
}// end loop over 't'
4037
}
4038
4039
if (combinations[r][s] == 2)
4040
{
4041
for (unsigned int t = 0; t < 4; t++)
4042
{
4043
for (unsigned int u = 0; u < 4; u++)
4044
{
4045
for (unsigned int tu = 0; tu < 4; tu++)
4046
{
4047
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
4048
}// end loop over 'tu'
4049
}// end loop over 'u'
4050
}// end loop over 't'
4051
}
4052
4053
}// end loop over 's'
4054
for (unsigned int s = 0; s < 4; s++)
4055
{
4056
for (unsigned int t = 0; t < 4; t++)
4057
{
4058
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
4059
}// end loop over 't'
4060
}// end loop over 's'
4061
}// end loop over 'r'
4062
4063
// Transform derivatives back to physical element
4064
for (unsigned int r = 0; r < num_derivatives; r++)
4065
{
4066
for (unsigned int s = 0; s < num_derivatives; s++)
4067
{
4068
values[num_derivatives + r] += transform[r][s]*derivatives[s];
4069
}// end loop over 's'
4070
}// end loop over 'r'
4071
4072
// Delete pointer to array of derivatives on FIAT element
4073
delete [] derivatives;
4074
4075
// Delete pointer to array of combinations of derivatives and transform
4076
for (unsigned int r = 0; r < num_derivatives; r++)
4077
{
4078
delete [] combinations[r];
4079
}// end loop over 'r'
4080
delete [] combinations;
4081
for (unsigned int r = 0; r < num_derivatives; r++)
4082
{
4083
delete [] transform[r];
4084
}// end loop over 'r'
4085
delete [] transform;
4086
break;
4087
}
4088
case 7:
4089
{
4090
4091
// Array of basisvalues.
4092
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
4093
4094
// Declare helper variables.
4095
unsigned int rr = 0;
4096
unsigned int ss = 0;
4097
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
4098
4099
// Compute basisvalues.
4100
basisvalues[0] = 1.00000000;
4101
basisvalues[1] = tmp0;
4102
for (unsigned int r = 0; r < 1; r++)
4103
{
4104
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
4105
ss = r*(r + 1)*(r + 2)/6;
4106
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
4107
}// end loop over 'r'
4108
for (unsigned int r = 0; r < 1; r++)
4109
{
4110
for (unsigned int s = 0; s < 1 - r; s++)
4111
{
4112
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
4113
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
4114
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
4115
}// end loop over 's'
4116
}// end loop over 'r'
4117
for (unsigned int r = 0; r < 2; r++)
4118
{
4119
for (unsigned int s = 0; s < 2 - r; s++)
4120
{
4121
for (unsigned int t = 0; t < 2 - r - s; t++)
4122
{
4123
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
4124
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
4125
}// end loop over 't'
4126
}// end loop over 's'
4127
}// end loop over 'r'
4128
4129
// Table(s) of coefficients.
4130
static const double coefficients0[4] = \
4131
{0.28867513, 0.00000000, 0.00000000, 0.22360680};
4132
4133
// Tables of derivatives of the polynomial base (transpose).
4134
static const double dmats0[4][4] = \
4135
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4136
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
4137
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4138
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
4139
4140
static const double dmats1[4][4] = \
4141
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4142
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
4143
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
4144
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
4145
4146
static const double dmats2[4][4] = \
4147
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4148
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
4149
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
4150
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
4151
4152
// Compute reference derivatives.
4153
// Declare pointer to array of derivatives on FIAT element.
4154
double *derivatives = new double[num_derivatives];
4155
for (unsigned int r = 0; r < num_derivatives; r++)
4156
{
4157
derivatives[r] = 0.00000000;
4158
}// end loop over 'r'
4159
4160
// Declare derivative matrix (of polynomial basis).
4161
double dmats[4][4] = \
4162
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
4163
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
4164
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
4165
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
4166
4167
// Declare (auxiliary) derivative matrix (of polynomial basis).
4168
double dmats_old[4][4] = \
4169
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
4170
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
4171
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
4172
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
4173
4174
// Loop possible derivatives.
4175
for (unsigned int r = 0; r < num_derivatives; r++)
4176
{
4177
// Resetting dmats values to compute next derivative.
4178
for (unsigned int t = 0; t < 4; t++)
4179
{
4180
for (unsigned int u = 0; u < 4; u++)
4181
{
4182
dmats[t][u] = 0.00000000;
4183
if (t == u)
4184
{
4185
dmats[t][u] = 1.00000000;
4186
}
4187
4188
}// end loop over 'u'
4189
}// end loop over 't'
4190
4191
// Looping derivative order to generate dmats.
4192
for (unsigned int s = 0; s < n; s++)
4193
{
4194
// Updating dmats_old with new values and resetting dmats.
4195
for (unsigned int t = 0; t < 4; t++)
4196
{
4197
for (unsigned int u = 0; u < 4; u++)
4198
{
4199
dmats_old[t][u] = dmats[t][u];
4200
dmats[t][u] = 0.00000000;
4201
}// end loop over 'u'
4202
}// end loop over 't'
4203
4204
// Update dmats using an inner product.
4205
if (combinations[r][s] == 0)
4206
{
4207
for (unsigned int t = 0; t < 4; t++)
4208
{
4209
for (unsigned int u = 0; u < 4; u++)
4210
{
4211
for (unsigned int tu = 0; tu < 4; tu++)
4212
{
4213
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
4214
}// end loop over 'tu'
4215
}// end loop over 'u'
4216
}// end loop over 't'
4217
}
4218
4219
if (combinations[r][s] == 1)
4220
{
4221
for (unsigned int t = 0; t < 4; t++)
4222
{
4223
for (unsigned int u = 0; u < 4; u++)
4224
{
4225
for (unsigned int tu = 0; tu < 4; tu++)
4226
{
4227
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
4228
}// end loop over 'tu'
4229
}// end loop over 'u'
4230
}// end loop over 't'
4231
}
4232
4233
if (combinations[r][s] == 2)
4234
{
4235
for (unsigned int t = 0; t < 4; t++)
4236
{
4237
for (unsigned int u = 0; u < 4; u++)
4238
{
4239
for (unsigned int tu = 0; tu < 4; tu++)
4240
{
4241
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
4242
}// end loop over 'tu'
4243
}// end loop over 'u'
4244
}// end loop over 't'
4245
}
4246
4247
}// end loop over 's'
4248
for (unsigned int s = 0; s < 4; s++)
4249
{
4250
for (unsigned int t = 0; t < 4; t++)
4251
{
4252
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
4253
}// end loop over 't'
4254
}// end loop over 's'
4255
}// end loop over 'r'
4256
4257
// Transform derivatives back to physical element
4258
for (unsigned int r = 0; r < num_derivatives; r++)
4259
{
4260
for (unsigned int s = 0; s < num_derivatives; s++)
4261
{
4262
values[num_derivatives + r] += transform[r][s]*derivatives[s];
4263
}// end loop over 's'
4264
}// end loop over 'r'
4265
4266
// Delete pointer to array of derivatives on FIAT element
4267
delete [] derivatives;
4268
4269
// Delete pointer to array of combinations of derivatives and transform
4270
for (unsigned int r = 0; r < num_derivatives; r++)
4271
{
4272
delete [] combinations[r];
4273
}// end loop over 'r'
4274
delete [] combinations;
4275
for (unsigned int r = 0; r < num_derivatives; r++)
4276
{
4277
delete [] transform[r];
4278
}// end loop over 'r'
4279
delete [] transform;
4280
break;
4281
}
4282
case 8:
4283
{
4284
4285
// Array of basisvalues.
4286
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
4287
4288
// Declare helper variables.
4289
unsigned int rr = 0;
4290
unsigned int ss = 0;
4291
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
4292
4293
// Compute basisvalues.
4294
basisvalues[0] = 1.00000000;
4295
basisvalues[1] = tmp0;
4296
for (unsigned int r = 0; r < 1; r++)
4297
{
4298
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
4299
ss = r*(r + 1)*(r + 2)/6;
4300
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
4301
}// end loop over 'r'
4302
for (unsigned int r = 0; r < 1; r++)
4303
{
4304
for (unsigned int s = 0; s < 1 - r; s++)
4305
{
4306
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
4307
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
4308
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
4309
}// end loop over 's'
4310
}// end loop over 'r'
4311
for (unsigned int r = 0; r < 2; r++)
4312
{
4313
for (unsigned int s = 0; s < 2 - r; s++)
4314
{
4315
for (unsigned int t = 0; t < 2 - r - s; t++)
4316
{
4317
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
4318
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
4319
}// end loop over 't'
4320
}// end loop over 's'
4321
}// end loop over 'r'
4322
4323
// Table(s) of coefficients.
4324
static const double coefficients0[4] = \
4325
{0.28867513, -0.18257419, -0.10540926, -0.07453560};
4326
4327
// Tables of derivatives of the polynomial base (transpose).
4328
static const double dmats0[4][4] = \
4329
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4330
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
4331
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4332
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
4333
4334
static const double dmats1[4][4] = \
4335
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4336
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
4337
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
4338
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
4339
4340
static const double dmats2[4][4] = \
4341
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4342
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
4343
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
4344
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
4345
4346
// Compute reference derivatives.
4347
// Declare pointer to array of derivatives on FIAT element.
4348
double *derivatives = new double[num_derivatives];
4349
for (unsigned int r = 0; r < num_derivatives; r++)
4350
{
4351
derivatives[r] = 0.00000000;
4352
}// end loop over 'r'
4353
4354
// Declare derivative matrix (of polynomial basis).
4355
double dmats[4][4] = \
4356
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
4357
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
4358
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
4359
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
4360
4361
// Declare (auxiliary) derivative matrix (of polynomial basis).
4362
double dmats_old[4][4] = \
4363
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
4364
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
4365
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
4366
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
4367
4368
// Loop possible derivatives.
4369
for (unsigned int r = 0; r < num_derivatives; r++)
4370
{
4371
// Resetting dmats values to compute next derivative.
4372
for (unsigned int t = 0; t < 4; t++)
4373
{
4374
for (unsigned int u = 0; u < 4; u++)
4375
{
4376
dmats[t][u] = 0.00000000;
4377
if (t == u)
4378
{
4379
dmats[t][u] = 1.00000000;
4380
}
4381
4382
}// end loop over 'u'
4383
}// end loop over 't'
4384
4385
// Looping derivative order to generate dmats.
4386
for (unsigned int s = 0; s < n; s++)
4387
{
4388
// Updating dmats_old with new values and resetting dmats.
4389
for (unsigned int t = 0; t < 4; t++)
4390
{
4391
for (unsigned int u = 0; u < 4; u++)
4392
{
4393
dmats_old[t][u] = dmats[t][u];
4394
dmats[t][u] = 0.00000000;
4395
}// end loop over 'u'
4396
}// end loop over 't'
4397
4398
// Update dmats using an inner product.
4399
if (combinations[r][s] == 0)
4400
{
4401
for (unsigned int t = 0; t < 4; t++)
4402
{
4403
for (unsigned int u = 0; u < 4; u++)
4404
{
4405
for (unsigned int tu = 0; tu < 4; tu++)
4406
{
4407
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
4408
}// end loop over 'tu'
4409
}// end loop over 'u'
4410
}// end loop over 't'
4411
}
4412
4413
if (combinations[r][s] == 1)
4414
{
4415
for (unsigned int t = 0; t < 4; t++)
4416
{
4417
for (unsigned int u = 0; u < 4; u++)
4418
{
4419
for (unsigned int tu = 0; tu < 4; tu++)
4420
{
4421
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
4422
}// end loop over 'tu'
4423
}// end loop over 'u'
4424
}// end loop over 't'
4425
}
4426
4427
if (combinations[r][s] == 2)
4428
{
4429
for (unsigned int t = 0; t < 4; t++)
4430
{
4431
for (unsigned int u = 0; u < 4; u++)
4432
{
4433
for (unsigned int tu = 0; tu < 4; tu++)
4434
{
4435
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
4436
}// end loop over 'tu'
4437
}// end loop over 'u'
4438
}// end loop over 't'
4439
}
4440
4441
}// end loop over 's'
4442
for (unsigned int s = 0; s < 4; s++)
4443
{
4444
for (unsigned int t = 0; t < 4; t++)
4445
{
4446
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
4447
}// end loop over 't'
4448
}// end loop over 's'
4449
}// end loop over 'r'
4450
4451
// Transform derivatives back to physical element
4452
for (unsigned int r = 0; r < num_derivatives; r++)
4453
{
4454
for (unsigned int s = 0; s < num_derivatives; s++)
4455
{
4456
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
4457
}// end loop over 's'
4458
}// end loop over 'r'
4459
4460
// Delete pointer to array of derivatives on FIAT element
4461
delete [] derivatives;
4462
4463
// Delete pointer to array of combinations of derivatives and transform
4464
for (unsigned int r = 0; r < num_derivatives; r++)
4465
{
4466
delete [] combinations[r];
4467
}// end loop over 'r'
4468
delete [] combinations;
4469
for (unsigned int r = 0; r < num_derivatives; r++)
4470
{
4471
delete [] transform[r];
4472
}// end loop over 'r'
4473
delete [] transform;
4474
break;
4475
}
4476
case 9:
4477
{
4478
4479
// Array of basisvalues.
4480
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
4481
4482
// Declare helper variables.
4483
unsigned int rr = 0;
4484
unsigned int ss = 0;
4485
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
4486
4487
// Compute basisvalues.
4488
basisvalues[0] = 1.00000000;
4489
basisvalues[1] = tmp0;
4490
for (unsigned int r = 0; r < 1; r++)
4491
{
4492
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
4493
ss = r*(r + 1)*(r + 2)/6;
4494
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
4495
}// end loop over 'r'
4496
for (unsigned int r = 0; r < 1; r++)
4497
{
4498
for (unsigned int s = 0; s < 1 - r; s++)
4499
{
4500
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
4501
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
4502
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
4503
}// end loop over 's'
4504
}// end loop over 'r'
4505
for (unsigned int r = 0; r < 2; r++)
4506
{
4507
for (unsigned int s = 0; s < 2 - r; s++)
4508
{
4509
for (unsigned int t = 0; t < 2 - r - s; t++)
4510
{
4511
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
4512
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
4513
}// end loop over 't'
4514
}// end loop over 's'
4515
}// end loop over 'r'
4516
4517
// Table(s) of coefficients.
4518
static const double coefficients0[4] = \
4519
{0.28867513, 0.18257419, -0.10540926, -0.07453560};
4520
4521
// Tables of derivatives of the polynomial base (transpose).
4522
static const double dmats0[4][4] = \
4523
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4524
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
4525
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4526
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
4527
4528
static const double dmats1[4][4] = \
4529
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4530
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
4531
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
4532
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
4533
4534
static const double dmats2[4][4] = \
4535
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4536
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
4537
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
4538
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
4539
4540
// Compute reference derivatives.
4541
// Declare pointer to array of derivatives on FIAT element.
4542
double *derivatives = new double[num_derivatives];
4543
for (unsigned int r = 0; r < num_derivatives; r++)
4544
{
4545
derivatives[r] = 0.00000000;
4546
}// end loop over 'r'
4547
4548
// Declare derivative matrix (of polynomial basis).
4549
double dmats[4][4] = \
4550
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
4551
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
4552
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
4553
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
4554
4555
// Declare (auxiliary) derivative matrix (of polynomial basis).
4556
double dmats_old[4][4] = \
4557
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
4558
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
4559
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
4560
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
4561
4562
// Loop possible derivatives.
4563
for (unsigned int r = 0; r < num_derivatives; r++)
4564
{
4565
// Resetting dmats values to compute next derivative.
4566
for (unsigned int t = 0; t < 4; t++)
4567
{
4568
for (unsigned int u = 0; u < 4; u++)
4569
{
4570
dmats[t][u] = 0.00000000;
4571
if (t == u)
4572
{
4573
dmats[t][u] = 1.00000000;
4574
}
4575
4576
}// end loop over 'u'
4577
}// end loop over 't'
4578
4579
// Looping derivative order to generate dmats.
4580
for (unsigned int s = 0; s < n; s++)
4581
{
4582
// Updating dmats_old with new values and resetting dmats.
4583
for (unsigned int t = 0; t < 4; t++)
4584
{
4585
for (unsigned int u = 0; u < 4; u++)
4586
{
4587
dmats_old[t][u] = dmats[t][u];
4588
dmats[t][u] = 0.00000000;
4589
}// end loop over 'u'
4590
}// end loop over 't'
4591
4592
// Update dmats using an inner product.
4593
if (combinations[r][s] == 0)
4594
{
4595
for (unsigned int t = 0; t < 4; t++)
4596
{
4597
for (unsigned int u = 0; u < 4; u++)
4598
{
4599
for (unsigned int tu = 0; tu < 4; tu++)
4600
{
4601
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
4602
}// end loop over 'tu'
4603
}// end loop over 'u'
4604
}// end loop over 't'
4605
}
4606
4607
if (combinations[r][s] == 1)
4608
{
4609
for (unsigned int t = 0; t < 4; t++)
4610
{
4611
for (unsigned int u = 0; u < 4; u++)
4612
{
4613
for (unsigned int tu = 0; tu < 4; tu++)
4614
{
4615
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
4616
}// end loop over 'tu'
4617
}// end loop over 'u'
4618
}// end loop over 't'
4619
}
4620
4621
if (combinations[r][s] == 2)
4622
{
4623
for (unsigned int t = 0; t < 4; t++)
4624
{
4625
for (unsigned int u = 0; u < 4; u++)
4626
{
4627
for (unsigned int tu = 0; tu < 4; tu++)
4628
{
4629
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
4630
}// end loop over 'tu'
4631
}// end loop over 'u'
4632
}// end loop over 't'
4633
}
4634
4635
}// end loop over 's'
4636
for (unsigned int s = 0; s < 4; s++)
4637
{
4638
for (unsigned int t = 0; t < 4; t++)
4639
{
4640
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
4641
}// end loop over 't'
4642
}// end loop over 's'
4643
}// end loop over 'r'
4644
4645
// Transform derivatives back to physical element
4646
for (unsigned int r = 0; r < num_derivatives; r++)
4647
{
4648
for (unsigned int s = 0; s < num_derivatives; s++)
4649
{
4650
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
4651
}// end loop over 's'
4652
}// end loop over 'r'
4653
4654
// Delete pointer to array of derivatives on FIAT element
4655
delete [] derivatives;
4656
4657
// Delete pointer to array of combinations of derivatives and transform
4658
for (unsigned int r = 0; r < num_derivatives; r++)
4659
{
4660
delete [] combinations[r];
4661
}// end loop over 'r'
4662
delete [] combinations;
4663
for (unsigned int r = 0; r < num_derivatives; r++)
4664
{
4665
delete [] transform[r];
4666
}// end loop over 'r'
4667
delete [] transform;
4668
break;
4669
}
4670
case 10:
4671
{
4672
4673
// Array of basisvalues.
4674
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
4675
4676
// Declare helper variables.
4677
unsigned int rr = 0;
4678
unsigned int ss = 0;
4679
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
4680
4681
// Compute basisvalues.
4682
basisvalues[0] = 1.00000000;
4683
basisvalues[1] = tmp0;
4684
for (unsigned int r = 0; r < 1; r++)
4685
{
4686
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
4687
ss = r*(r + 1)*(r + 2)/6;
4688
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
4689
}// end loop over 'r'
4690
for (unsigned int r = 0; r < 1; r++)
4691
{
4692
for (unsigned int s = 0; s < 1 - r; s++)
4693
{
4694
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
4695
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
4696
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
4697
}// end loop over 's'
4698
}// end loop over 'r'
4699
for (unsigned int r = 0; r < 2; r++)
4700
{
4701
for (unsigned int s = 0; s < 2 - r; s++)
4702
{
4703
for (unsigned int t = 0; t < 2 - r - s; t++)
4704
{
4705
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
4706
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
4707
}// end loop over 't'
4708
}// end loop over 's'
4709
}// end loop over 'r'
4710
4711
// Table(s) of coefficients.
4712
static const double coefficients0[4] = \
4713
{0.28867513, 0.00000000, 0.21081851, -0.07453560};
4714
4715
// Tables of derivatives of the polynomial base (transpose).
4716
static const double dmats0[4][4] = \
4717
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4718
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
4719
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4720
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
4721
4722
static const double dmats1[4][4] = \
4723
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4724
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
4725
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
4726
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
4727
4728
static const double dmats2[4][4] = \
4729
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4730
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
4731
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
4732
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
4733
4734
// Compute reference derivatives.
4735
// Declare pointer to array of derivatives on FIAT element.
4736
double *derivatives = new double[num_derivatives];
4737
for (unsigned int r = 0; r < num_derivatives; r++)
4738
{
4739
derivatives[r] = 0.00000000;
4740
}// end loop over 'r'
4741
4742
// Declare derivative matrix (of polynomial basis).
4743
double dmats[4][4] = \
4744
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
4745
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
4746
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
4747
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
4748
4749
// Declare (auxiliary) derivative matrix (of polynomial basis).
4750
double dmats_old[4][4] = \
4751
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
4752
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
4753
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
4754
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
4755
4756
// Loop possible derivatives.
4757
for (unsigned int r = 0; r < num_derivatives; r++)
4758
{
4759
// Resetting dmats values to compute next derivative.
4760
for (unsigned int t = 0; t < 4; t++)
4761
{
4762
for (unsigned int u = 0; u < 4; u++)
4763
{
4764
dmats[t][u] = 0.00000000;
4765
if (t == u)
4766
{
4767
dmats[t][u] = 1.00000000;
4768
}
4769
4770
}// end loop over 'u'
4771
}// end loop over 't'
4772
4773
// Looping derivative order to generate dmats.
4774
for (unsigned int s = 0; s < n; s++)
4775
{
4776
// Updating dmats_old with new values and resetting dmats.
4777
for (unsigned int t = 0; t < 4; t++)
4778
{
4779
for (unsigned int u = 0; u < 4; u++)
4780
{
4781
dmats_old[t][u] = dmats[t][u];
4782
dmats[t][u] = 0.00000000;
4783
}// end loop over 'u'
4784
}// end loop over 't'
4785
4786
// Update dmats using an inner product.
4787
if (combinations[r][s] == 0)
4788
{
4789
for (unsigned int t = 0; t < 4; t++)
4790
{
4791
for (unsigned int u = 0; u < 4; u++)
4792
{
4793
for (unsigned int tu = 0; tu < 4; tu++)
4794
{
4795
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
4796
}// end loop over 'tu'
4797
}// end loop over 'u'
4798
}// end loop over 't'
4799
}
4800
4801
if (combinations[r][s] == 1)
4802
{
4803
for (unsigned int t = 0; t < 4; t++)
4804
{
4805
for (unsigned int u = 0; u < 4; u++)
4806
{
4807
for (unsigned int tu = 0; tu < 4; tu++)
4808
{
4809
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
4810
}// end loop over 'tu'
4811
}// end loop over 'u'
4812
}// end loop over 't'
4813
}
4814
4815
if (combinations[r][s] == 2)
4816
{
4817
for (unsigned int t = 0; t < 4; t++)
4818
{
4819
for (unsigned int u = 0; u < 4; u++)
4820
{
4821
for (unsigned int tu = 0; tu < 4; tu++)
4822
{
4823
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
4824
}// end loop over 'tu'
4825
}// end loop over 'u'
4826
}// end loop over 't'
4827
}
4828
4829
}// end loop over 's'
4830
for (unsigned int s = 0; s < 4; s++)
4831
{
4832
for (unsigned int t = 0; t < 4; t++)
4833
{
4834
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
4835
}// end loop over 't'
4836
}// end loop over 's'
4837
}// end loop over 'r'
4838
4839
// Transform derivatives back to physical element
4840
for (unsigned int r = 0; r < num_derivatives; r++)
4841
{
4842
for (unsigned int s = 0; s < num_derivatives; s++)
4843
{
4844
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
4845
}// end loop over 's'
4846
}// end loop over 'r'
4847
4848
// Delete pointer to array of derivatives on FIAT element
4849
delete [] derivatives;
4850
4851
// Delete pointer to array of combinations of derivatives and transform
4852
for (unsigned int r = 0; r < num_derivatives; r++)
4853
{
4854
delete [] combinations[r];
4855
}// end loop over 'r'
4856
delete [] combinations;
4857
for (unsigned int r = 0; r < num_derivatives; r++)
4858
{
4859
delete [] transform[r];
4860
}// end loop over 'r'
4861
delete [] transform;
4862
break;
4863
}
4864
case 11:
4865
{
4866
4867
// Array of basisvalues.
4868
double basisvalues[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
4869
4870
// Declare helper variables.
4871
unsigned int rr = 0;
4872
unsigned int ss = 0;
4873
double tmp0 = 0.50000000*(2.00000000 + Y + Z + 2.00000000*X);
4874
4875
// Compute basisvalues.
4876
basisvalues[0] = 1.00000000;
4877
basisvalues[1] = tmp0;
4878
for (unsigned int r = 0; r < 1; r++)
4879
{
4880
rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
4881
ss = r*(r + 1)*(r + 2)/6;
4882
basisvalues[rr] = basisvalues[ss]*(r*(1.00000000 + Y) + (2.00000000 + Z + 3.00000000*Y)/2.00000000);
4883
}// end loop over 'r'
4884
for (unsigned int r = 0; r < 1; r++)
4885
{
4886
for (unsigned int s = 0; s < 1 - r; s++)
4887
{
4888
rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
4889
ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
4890
basisvalues[rr] = basisvalues[ss]*(1.00000000 + r + s + Z*(2.00000000 + r + s));
4891
}// end loop over 's'
4892
}// end loop over 'r'
4893
for (unsigned int r = 0; r < 2; r++)
4894
{
4895
for (unsigned int s = 0; s < 2 - r; s++)
4896
{
4897
for (unsigned int t = 0; t < 2 - r - s; t++)
4898
{
4899
rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
4900
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s)*(1.50000000 + r + s + t));
4901
}// end loop over 't'
4902
}// end loop over 's'
4903
}// end loop over 'r'
4904
4905
// Table(s) of coefficients.
4906
static const double coefficients0[4] = \
4907
{0.28867513, 0.00000000, 0.00000000, 0.22360680};
4908
4909
// Tables of derivatives of the polynomial base (transpose).
4910
static const double dmats0[4][4] = \
4911
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4912
{6.32455532, 0.00000000, 0.00000000, 0.00000000},
4913
{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4914
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
4915
4916
static const double dmats1[4][4] = \
4917
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4918
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
4919
{5.47722558, 0.00000000, 0.00000000, 0.00000000},
4920
{0.00000000, 0.00000000, 0.00000000, 0.00000000}};
4921
4922
static const double dmats2[4][4] = \
4923
{{0.00000000, 0.00000000, 0.00000000, 0.00000000},
4924
{3.16227766, 0.00000000, 0.00000000, 0.00000000},
4925
{1.82574186, 0.00000000, 0.00000000, 0.00000000},
4926
{5.16397779, 0.00000000, 0.00000000, 0.00000000}};
4927
4928
// Compute reference derivatives.
4929
// Declare pointer to array of derivatives on FIAT element.
4930
double *derivatives = new double[num_derivatives];
4931
for (unsigned int r = 0; r < num_derivatives; r++)
4932
{
4933
derivatives[r] = 0.00000000;
4934
}// end loop over 'r'
4935
4936
// Declare derivative matrix (of polynomial basis).
4937
double dmats[4][4] = \
4938
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
4939
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
4940
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
4941
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
4942
4943
// Declare (auxiliary) derivative matrix (of polynomial basis).
4944
double dmats_old[4][4] = \
4945
{{1.00000000, 0.00000000, 0.00000000, 0.00000000},
4946
{0.00000000, 1.00000000, 0.00000000, 0.00000000},
4947
{0.00000000, 0.00000000, 1.00000000, 0.00000000},
4948
{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
4949
4950
// Loop possible derivatives.
4951
for (unsigned int r = 0; r < num_derivatives; r++)
4952
{
4953
// Resetting dmats values to compute next derivative.
4954
for (unsigned int t = 0; t < 4; t++)
4955
{
4956
for (unsigned int u = 0; u < 4; u++)
4957
{
4958
dmats[t][u] = 0.00000000;
4959
if (t == u)
4960
{
4961
dmats[t][u] = 1.00000000;
4962
}
4963
4964
}// end loop over 'u'
4965
}// end loop over 't'
4966
4967
// Looping derivative order to generate dmats.
4968
for (unsigned int s = 0; s < n; s++)
4969
{
4970
// Updating dmats_old with new values and resetting dmats.
4971
for (unsigned int t = 0; t < 4; t++)
4972
{
4973
for (unsigned int u = 0; u < 4; u++)
4974
{
4975
dmats_old[t][u] = dmats[t][u];
4976
dmats[t][u] = 0.00000000;
4977
}// end loop over 'u'
4978
}// end loop over 't'
4979
4980
// Update dmats using an inner product.
4981
if (combinations[r][s] == 0)
4982
{
4983
for (unsigned int t = 0; t < 4; t++)
4984
{
4985
for (unsigned int u = 0; u < 4; u++)
4986
{
4987
for (unsigned int tu = 0; tu < 4; tu++)
4988
{
4989
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
4990
}// end loop over 'tu'
4991
}// end loop over 'u'
4992
}// end loop over 't'
4993
}
4994
4995
if (combinations[r][s] == 1)
4996
{
4997
for (unsigned int t = 0; t < 4; t++)
4998
{
4999
for (unsigned int u = 0; u < 4; u++)
5000
{
5001
for (unsigned int tu = 0; tu < 4; tu++)
5002
{
5003
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
5004
}// end loop over 'tu'
5005
}// end loop over 'u'
5006
}// end loop over 't'
5007
}
5008
5009
if (combinations[r][s] == 2)
5010
{
5011
for (unsigned int t = 0; t < 4; t++)
5012
{
5013
for (unsigned int u = 0; u < 4; u++)
5014
{
5015
for (unsigned int tu = 0; tu < 4; tu++)
5016
{
5017
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
5018
}// end loop over 'tu'
5019
}// end loop over 'u'
5020
}// end loop over 't'
5021
}
5022
5023
}// end loop over 's'
5024
for (unsigned int s = 0; s < 4; s++)
5025
{
5026
for (unsigned int t = 0; t < 4; t++)
5027
{
5028
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
5029
}// end loop over 't'
5030
}// end loop over 's'
5031
}// end loop over 'r'
5032
5033
// Transform derivatives back to physical element
5034
for (unsigned int r = 0; r < num_derivatives; r++)
5035
{
5036
for (unsigned int s = 0; s < num_derivatives; s++)
5037
{
5038
values[2*num_derivatives + r] += transform[r][s]*derivatives[s];
5039
}// end loop over 's'
5040
}// end loop over 'r'
5041
5042
// Delete pointer to array of derivatives on FIAT element
5043
delete [] derivatives;
5044
5045
// Delete pointer to array of combinations of derivatives and transform
5046
for (unsigned int r = 0; r < num_derivatives; r++)
5047
{
5048
delete [] combinations[r];
5049
}// end loop over 'r'
5050
delete [] combinations;
5051
for (unsigned int r = 0; r < num_derivatives; r++)
5052
{
5053
delete [] transform[r];
5054
}// end loop over 'r'
5055
delete [] transform;
5056
break;
5057
}
2128
5058
}
2129
5059
2130
5060
}
3419
6349
}// end loop over 'r'
3420
6350
3421
6351
// Compute element tensor using UFL quadrature representation
3422
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
6352
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
3423
6353
3424
6354
// Loop quadrature points for integral.
3425
6355
// Number of operations to compute element tensor for following IP loop = 802728
3662
6592
}// end loop over 'r'
3663
6593
3664
6594
// Compute element tensor using UFL quadrature representation
3665
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
6595
// Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
3666
6596
3667
6597
// Loop quadrature points for integral.
3668
6598
// Number of operations to compute element tensor for following IP loop = 242496
Loggerhead is a web-based interface for Breezy
Version: 2.0.1