~ubuntu-branches/ubuntu/raring/ffc/raring

« back to all changes in this revision

Viewing changes to test/regression/references/PoissonDG.h

  • 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.

Show diffs side-by-side

added added

removed removed

Lines of Context:
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
// 
96
96
    // Reset values.
97
97
    *values = 0.00000000;
98
98
    
99
 
    // Map degree of freedom to element degree of freedom
100
 
    const unsigned int dof = i;
101
 
    
102
99
    // Array of basisvalues.
103
100
    double basisvalues[1] = {0.00000000};
104
101
    
108
105
    basisvalues[0] = 1.00000000;
109
106
    
110
107
    // Table(s) of coefficients.
111
 
    static const double coefficients0[1][1] = \
112
 
    {{1.00000000}};
 
108
    static const double coefficients0[1] = \
 
109
    {1.00000000};
113
110
    
114
111
    // Compute value(s).
115
112
    for (unsigned int r = 0; r < 1; r++)
116
113
    {
117
 
      *values += coefficients0[dof][r]*basisvalues[r];
 
114
      *values += coefficients0[r]*basisvalues[r];
118
115
    }// end loop over 'r'
119
116
  }
120
117
 
220
217
      values[r] = 0.00000000;
221
218
    }// end loop over 'r'
222
219
    
223
 
    // Map degree of freedom to element degree of freedom
224
 
    const unsigned int dof = i;
225
220
    
226
221
    // Array of basisvalues.
227
222
    double basisvalues[1] = {0.00000000};
232
227
    basisvalues[0] = 1.00000000;
233
228
    
234
229
    // Table(s) of coefficients.
235
 
    static const double coefficients0[1][1] = \
236
 
    {{1.00000000}};
 
230
    static const double coefficients0[1] = \
 
231
    {1.00000000};
237
232
    
238
233
    // Tables of derivatives of the polynomial base (transpose).
239
234
    static const double dmats0[1][1] = \
322
317
      {
323
318
        for (unsigned int t = 0; t < 1; t++)
324
319
        {
325
 
          derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
 
320
          derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
326
321
        }// end loop over 't'
327
322
      }// end loop over 's'
328
323
    }// end loop over 'r'
508
503
    
509
504
    // Reset values.
510
505
    *values = 0.00000000;
511
 
    
512
 
    // Map degree of freedom to element degree of freedom
513
 
    const unsigned int dof = i;
514
 
    
515
 
    // Array of basisvalues.
516
 
    double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
517
 
    
518
 
    // Declare helper variables.
519
 
    unsigned int rr = 0;
520
 
    unsigned int ss = 0;
521
 
    double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
522
 
    
523
 
    // Compute basisvalues.
524
 
    basisvalues[0] = 1.00000000;
525
 
    basisvalues[1] = tmp0;
526
 
    for (unsigned int r = 0; r < 1; r++)
527
 
    {
528
 
      rr = (r + 1)*(r + 1 + 1)/2 + 1;
529
 
      ss = r*(r + 1)/2;
530
 
      basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
531
 
    }// end loop over 'r'
532
 
    for (unsigned int r = 0; r < 2; r++)
533
 
    {
534
 
      for (unsigned int s = 0; s < 2 - r; s++)
535
 
      {
536
 
        rr = (r + s)*(r + s + 1)/2 + s;
537
 
        basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
538
 
      }// end loop over 's'
539
 
    }// end loop over 'r'
540
 
    
541
 
    // Table(s) of coefficients.
542
 
    static const double coefficients0[3][3] = \
543
 
    {{0.47140452, -0.28867513, -0.16666667},
544
 
    {0.47140452, 0.28867513, -0.16666667},
545
 
    {0.47140452, 0.00000000, 0.33333333}};
546
 
    
547
 
    // Compute value(s).
548
 
    for (unsigned int r = 0; r < 3; r++)
549
 
    {
550
 
      *values += coefficients0[dof][r]*basisvalues[r];
551
 
    }// end loop over 'r'
 
506
    switch (i)
 
507
    {
 
508
    case 0:
 
509
      {
 
510
        
 
511
      // Array of basisvalues.
 
512
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
513
      
 
514
      // Declare helper variables.
 
515
      unsigned int rr = 0;
 
516
      unsigned int ss = 0;
 
517
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
518
      
 
519
      // Compute basisvalues.
 
520
      basisvalues[0] = 1.00000000;
 
521
      basisvalues[1] = tmp0;
 
522
      for (unsigned int r = 0; r < 1; r++)
 
523
      {
 
524
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
525
        ss = r*(r + 1)/2;
 
526
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
527
      }// end loop over 'r'
 
528
      for (unsigned int r = 0; r < 2; r++)
 
529
      {
 
530
        for (unsigned int s = 0; s < 2 - r; s++)
 
531
        {
 
532
          rr = (r + s)*(r + s + 1)/2 + s;
 
533
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
534
        }// end loop over 's'
 
535
      }// end loop over 'r'
 
536
      
 
537
      // Table(s) of coefficients.
 
538
      static const double coefficients0[3] = \
 
539
      {0.47140452, -0.28867513, -0.16666667};
 
540
      
 
541
      // Compute value(s).
 
542
      for (unsigned int r = 0; r < 3; r++)
 
543
      {
 
544
        *values += coefficients0[r]*basisvalues[r];
 
545
      }// end loop over 'r'
 
546
        break;
 
547
      }
 
548
    case 1:
 
549
      {
 
550
        
 
551
      // Array of basisvalues.
 
552
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
553
      
 
554
      // Declare helper variables.
 
555
      unsigned int rr = 0;
 
556
      unsigned int ss = 0;
 
557
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
558
      
 
559
      // Compute basisvalues.
 
560
      basisvalues[0] = 1.00000000;
 
561
      basisvalues[1] = tmp0;
 
562
      for (unsigned int r = 0; r < 1; r++)
 
563
      {
 
564
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
565
        ss = r*(r + 1)/2;
 
566
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
567
      }// end loop over 'r'
 
568
      for (unsigned int r = 0; r < 2; r++)
 
569
      {
 
570
        for (unsigned int s = 0; s < 2 - r; s++)
 
571
        {
 
572
          rr = (r + s)*(r + s + 1)/2 + s;
 
573
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
574
        }// end loop over 's'
 
575
      }// end loop over 'r'
 
576
      
 
577
      // Table(s) of coefficients.
 
578
      static const double coefficients0[3] = \
 
579
      {0.47140452, 0.28867513, -0.16666667};
 
580
      
 
581
      // Compute value(s).
 
582
      for (unsigned int r = 0; r < 3; r++)
 
583
      {
 
584
        *values += coefficients0[r]*basisvalues[r];
 
585
      }// end loop over 'r'
 
586
        break;
 
587
      }
 
588
    case 2:
 
589
      {
 
590
        
 
591
      // Array of basisvalues.
 
592
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
593
      
 
594
      // Declare helper variables.
 
595
      unsigned int rr = 0;
 
596
      unsigned int ss = 0;
 
597
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
598
      
 
599
      // Compute basisvalues.
 
600
      basisvalues[0] = 1.00000000;
 
601
      basisvalues[1] = tmp0;
 
602
      for (unsigned int r = 0; r < 1; r++)
 
603
      {
 
604
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
605
        ss = r*(r + 1)/2;
 
606
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
607
      }// end loop over 'r'
 
608
      for (unsigned int r = 0; r < 2; r++)
 
609
      {
 
610
        for (unsigned int s = 0; s < 2 - r; s++)
 
611
        {
 
612
          rr = (r + s)*(r + s + 1)/2 + s;
 
613
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
614
        }// end loop over 's'
 
615
      }// end loop over 'r'
 
616
      
 
617
      // Table(s) of coefficients.
 
618
      static const double coefficients0[3] = \
 
619
      {0.47140452, 0.00000000, 0.33333333};
 
620
      
 
621
      // Compute value(s).
 
622
      for (unsigned int r = 0; r < 3; r++)
 
623
      {
 
624
        *values += coefficients0[r]*basisvalues[r];
 
625
      }// end loop over 'r'
 
626
        break;
 
627
      }
 
628
    }
 
629
    
552
630
  }
553
631
 
554
632
  /// Evaluate all basis functions at given point in cell
664
742
      values[r] = 0.00000000;
665
743
    }// end loop over 'r'
666
744
    
667
 
    // Map degree of freedom to element degree of freedom
668
 
    const unsigned int dof = i;
669
 
    
670
 
    // Array of basisvalues.
671
 
    double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
672
 
    
673
 
    // Declare helper variables.
674
 
    unsigned int rr = 0;
675
 
    unsigned int ss = 0;
676
 
    double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
677
 
    
678
 
    // Compute basisvalues.
679
 
    basisvalues[0] = 1.00000000;
680
 
    basisvalues[1] = tmp0;
681
 
    for (unsigned int r = 0; r < 1; r++)
682
 
    {
683
 
      rr = (r + 1)*(r + 1 + 1)/2 + 1;
684
 
      ss = r*(r + 1)/2;
685
 
      basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
686
 
    }// end loop over 'r'
687
 
    for (unsigned int r = 0; r < 2; r++)
688
 
    {
689
 
      for (unsigned int s = 0; s < 2 - r; s++)
690
 
      {
691
 
        rr = (r + s)*(r + s + 1)/2 + s;
692
 
        basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
693
 
      }// end loop over 's'
694
 
    }// end loop over 'r'
695
 
    
696
 
    // Table(s) of coefficients.
697
 
    static const double coefficients0[3][3] = \
698
 
    {{0.47140452, -0.28867513, -0.16666667},
699
 
    {0.47140452, 0.28867513, -0.16666667},
700
 
    {0.47140452, 0.00000000, 0.33333333}};
701
 
    
702
 
    // Tables of derivatives of the polynomial base (transpose).
703
 
    static const double dmats0[3][3] = \
704
 
    {{0.00000000, 0.00000000, 0.00000000},
705
 
    {4.89897949, 0.00000000, 0.00000000},
706
 
    {0.00000000, 0.00000000, 0.00000000}};
707
 
    
708
 
    static const double dmats1[3][3] = \
709
 
    {{0.00000000, 0.00000000, 0.00000000},
710
 
    {2.44948974, 0.00000000, 0.00000000},
711
 
    {4.24264069, 0.00000000, 0.00000000}};
712
 
    
713
 
    // Compute reference derivatives.
714
 
    // Declare pointer to array of derivatives on FIAT element.
715
 
    double *derivatives = new double[num_derivatives];
716
 
    for (unsigned int r = 0; r < num_derivatives; r++)
717
 
    {
718
 
      derivatives[r] = 0.00000000;
719
 
    }// end loop over 'r'
720
 
    
721
 
    // Declare derivative matrix (of polynomial basis).
722
 
    double dmats[3][3] = \
723
 
    {{1.00000000, 0.00000000, 0.00000000},
724
 
    {0.00000000, 1.00000000, 0.00000000},
725
 
    {0.00000000, 0.00000000, 1.00000000}};
726
 
    
727
 
    // Declare (auxiliary) derivative matrix (of polynomial basis).
728
 
    double dmats_old[3][3] = \
729
 
    {{1.00000000, 0.00000000, 0.00000000},
730
 
    {0.00000000, 1.00000000, 0.00000000},
731
 
    {0.00000000, 0.00000000, 1.00000000}};
732
 
    
733
 
    // Loop possible derivatives.
734
 
    for (unsigned int r = 0; r < num_derivatives; r++)
735
 
    {
736
 
      // Resetting dmats values to compute next derivative.
737
 
      for (unsigned int t = 0; t < 3; t++)
738
 
      {
739
 
        for (unsigned int u = 0; u < 3; u++)
740
 
        {
741
 
          dmats[t][u] = 0.00000000;
742
 
          if (t == u)
743
 
          {
744
 
          dmats[t][u] = 1.00000000;
745
 
          }
746
 
          
747
 
        }// end loop over 'u'
748
 
      }// end loop over 't'
749
 
      
750
 
      // Looping derivative order to generate dmats.
751
 
      for (unsigned int s = 0; s < n; s++)
752
 
      {
753
 
        // Updating dmats_old with new values and resetting dmats.
754
 
        for (unsigned int t = 0; t < 3; t++)
755
 
        {
756
 
          for (unsigned int u = 0; u < 3; u++)
757
 
          {
758
 
            dmats_old[t][u] = dmats[t][u];
759
 
            dmats[t][u] = 0.00000000;
760
 
          }// end loop over 'u'
761
 
        }// end loop over 't'
762
 
        
763
 
        // Update dmats using an inner product.
764
 
        if (combinations[r][s] == 0)
765
 
        {
766
 
        for (unsigned int t = 0; t < 3; t++)
767
 
        {
768
 
          for (unsigned int u = 0; u < 3; u++)
769
 
          {
770
 
            for (unsigned int tu = 0; tu < 3; tu++)
771
 
            {
772
 
              dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
773
 
            }// end loop over 'tu'
774
 
          }// end loop over 'u'
775
 
        }// end loop over 't'
776
 
        }
777
 
        
778
 
        if (combinations[r][s] == 1)
779
 
        {
780
 
        for (unsigned int t = 0; t < 3; t++)
781
 
        {
782
 
          for (unsigned int u = 0; u < 3; u++)
783
 
          {
784
 
            for (unsigned int tu = 0; tu < 3; tu++)
785
 
            {
786
 
              dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
787
 
            }// end loop over 'tu'
788
 
          }// end loop over 'u'
789
 
        }// end loop over 't'
790
 
        }
791
 
        
792
 
      }// end loop over 's'
793
 
      for (unsigned int s = 0; s < 3; s++)
794
 
      {
795
 
        for (unsigned int t = 0; t < 3; t++)
796
 
        {
797
 
          derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
798
 
        }// end loop over 't'
799
 
      }// end loop over 's'
800
 
    }// end loop over 'r'
801
 
    
802
 
    // Transform derivatives back to physical element
803
 
    for (unsigned int r = 0; r < num_derivatives; r++)
804
 
    {
805
 
      for (unsigned int s = 0; s < num_derivatives; s++)
806
 
      {
807
 
        values[r] += transform[r][s]*derivatives[s];
808
 
      }// end loop over 's'
809
 
    }// end loop over 'r'
810
 
    
811
 
    // Delete pointer to array of derivatives on FIAT element
812
 
    delete [] derivatives;
813
 
    
814
 
    // Delete pointer to array of combinations of derivatives and transform
815
 
    for (unsigned int r = 0; r < num_derivatives; r++)
816
 
    {
817
 
      delete [] combinations[r];
818
 
    }// end loop over 'r'
819
 
    delete [] combinations;
820
 
    for (unsigned int r = 0; r < num_derivatives; r++)
821
 
    {
822
 
      delete [] transform[r];
823
 
    }// end loop over 'r'
824
 
    delete [] transform;
 
745
    switch (i)
 
746
    {
 
747
    case 0:
 
748
      {
 
749
        
 
750
      // Array of basisvalues.
 
751
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
752
      
 
753
      // Declare helper variables.
 
754
      unsigned int rr = 0;
 
755
      unsigned int ss = 0;
 
756
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
757
      
 
758
      // Compute basisvalues.
 
759
      basisvalues[0] = 1.00000000;
 
760
      basisvalues[1] = tmp0;
 
761
      for (unsigned int r = 0; r < 1; r++)
 
762
      {
 
763
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
764
        ss = r*(r + 1)/2;
 
765
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
766
      }// end loop over 'r'
 
767
      for (unsigned int r = 0; r < 2; r++)
 
768
      {
 
769
        for (unsigned int s = 0; s < 2 - r; s++)
 
770
        {
 
771
          rr = (r + s)*(r + s + 1)/2 + s;
 
772
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
773
        }// end loop over 's'
 
774
      }// end loop over 'r'
 
775
      
 
776
      // Table(s) of coefficients.
 
777
      static const double coefficients0[3] = \
 
778
      {0.47140452, -0.28867513, -0.16666667};
 
779
      
 
780
      // Tables of derivatives of the polynomial base (transpose).
 
781
      static const double dmats0[3][3] = \
 
782
      {{0.00000000, 0.00000000, 0.00000000},
 
783
      {4.89897949, 0.00000000, 0.00000000},
 
784
      {0.00000000, 0.00000000, 0.00000000}};
 
785
      
 
786
      static const double dmats1[3][3] = \
 
787
      {{0.00000000, 0.00000000, 0.00000000},
 
788
      {2.44948974, 0.00000000, 0.00000000},
 
789
      {4.24264069, 0.00000000, 0.00000000}};
 
790
      
 
791
      // Compute reference derivatives.
 
792
      // Declare pointer to array of derivatives on FIAT element.
 
793
      double *derivatives = new double[num_derivatives];
 
794
      for (unsigned int r = 0; r < num_derivatives; r++)
 
795
      {
 
796
        derivatives[r] = 0.00000000;
 
797
      }// end loop over 'r'
 
798
      
 
799
      // Declare derivative matrix (of polynomial basis).
 
800
      double dmats[3][3] = \
 
801
      {{1.00000000, 0.00000000, 0.00000000},
 
802
      {0.00000000, 1.00000000, 0.00000000},
 
803
      {0.00000000, 0.00000000, 1.00000000}};
 
804
      
 
805
      // Declare (auxiliary) derivative matrix (of polynomial basis).
 
806
      double dmats_old[3][3] = \
 
807
      {{1.00000000, 0.00000000, 0.00000000},
 
808
      {0.00000000, 1.00000000, 0.00000000},
 
809
      {0.00000000, 0.00000000, 1.00000000}};
 
810
      
 
811
      // Loop possible derivatives.
 
812
      for (unsigned int r = 0; r < num_derivatives; r++)
 
813
      {
 
814
        // Resetting dmats values to compute next derivative.
 
815
        for (unsigned int t = 0; t < 3; t++)
 
816
        {
 
817
          for (unsigned int u = 0; u < 3; u++)
 
818
          {
 
819
            dmats[t][u] = 0.00000000;
 
820
            if (t == u)
 
821
            {
 
822
            dmats[t][u] = 1.00000000;
 
823
            }
 
824
            
 
825
          }// end loop over 'u'
 
826
        }// end loop over 't'
 
827
        
 
828
        // Looping derivative order to generate dmats.
 
829
        for (unsigned int s = 0; s < n; s++)
 
830
        {
 
831
          // Updating dmats_old with new values and resetting dmats.
 
832
          for (unsigned int t = 0; t < 3; t++)
 
833
          {
 
834
            for (unsigned int u = 0; u < 3; u++)
 
835
            {
 
836
              dmats_old[t][u] = dmats[t][u];
 
837
              dmats[t][u] = 0.00000000;
 
838
            }// end loop over 'u'
 
839
          }// end loop over 't'
 
840
          
 
841
          // Update dmats using an inner product.
 
842
          if (combinations[r][s] == 0)
 
843
          {
 
844
          for (unsigned int t = 0; t < 3; t++)
 
845
          {
 
846
            for (unsigned int u = 0; u < 3; u++)
 
847
            {
 
848
              for (unsigned int tu = 0; tu < 3; tu++)
 
849
              {
 
850
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
851
              }// end loop over 'tu'
 
852
            }// end loop over 'u'
 
853
          }// end loop over 't'
 
854
          }
 
855
          
 
856
          if (combinations[r][s] == 1)
 
857
          {
 
858
          for (unsigned int t = 0; t < 3; t++)
 
859
          {
 
860
            for (unsigned int u = 0; u < 3; u++)
 
861
            {
 
862
              for (unsigned int tu = 0; tu < 3; tu++)
 
863
              {
 
864
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
865
              }// end loop over 'tu'
 
866
            }// end loop over 'u'
 
867
          }// end loop over 't'
 
868
          }
 
869
          
 
870
        }// end loop over 's'
 
871
        for (unsigned int s = 0; s < 3; s++)
 
872
        {
 
873
          for (unsigned int t = 0; t < 3; t++)
 
874
          {
 
875
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
 
876
          }// end loop over 't'
 
877
        }// end loop over 's'
 
878
      }// end loop over 'r'
 
879
      
 
880
      // Transform derivatives back to physical element
 
881
      for (unsigned int r = 0; r < num_derivatives; r++)
 
882
      {
 
883
        for (unsigned int s = 0; s < num_derivatives; s++)
 
884
        {
 
885
          values[r] += transform[r][s]*derivatives[s];
 
886
        }// end loop over 's'
 
887
      }// end loop over 'r'
 
888
      
 
889
      // Delete pointer to array of derivatives on FIAT element
 
890
      delete [] derivatives;
 
891
      
 
892
      // Delete pointer to array of combinations of derivatives and transform
 
893
      for (unsigned int r = 0; r < num_derivatives; r++)
 
894
      {
 
895
        delete [] combinations[r];
 
896
      }// end loop over 'r'
 
897
      delete [] combinations;
 
898
      for (unsigned int r = 0; r < num_derivatives; r++)
 
899
      {
 
900
        delete [] transform[r];
 
901
      }// end loop over 'r'
 
902
      delete [] transform;
 
903
        break;
 
904
      }
 
905
    case 1:
 
906
      {
 
907
        
 
908
      // Array of basisvalues.
 
909
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
910
      
 
911
      // Declare helper variables.
 
912
      unsigned int rr = 0;
 
913
      unsigned int ss = 0;
 
914
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
915
      
 
916
      // Compute basisvalues.
 
917
      basisvalues[0] = 1.00000000;
 
918
      basisvalues[1] = tmp0;
 
919
      for (unsigned int r = 0; r < 1; r++)
 
920
      {
 
921
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
922
        ss = r*(r + 1)/2;
 
923
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
924
      }// end loop over 'r'
 
925
      for (unsigned int r = 0; r < 2; r++)
 
926
      {
 
927
        for (unsigned int s = 0; s < 2 - r; s++)
 
928
        {
 
929
          rr = (r + s)*(r + s + 1)/2 + s;
 
930
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
931
        }// end loop over 's'
 
932
      }// end loop over 'r'
 
933
      
 
934
      // Table(s) of coefficients.
 
935
      static const double coefficients0[3] = \
 
936
      {0.47140452, 0.28867513, -0.16666667};
 
937
      
 
938
      // Tables of derivatives of the polynomial base (transpose).
 
939
      static const double dmats0[3][3] = \
 
940
      {{0.00000000, 0.00000000, 0.00000000},
 
941
      {4.89897949, 0.00000000, 0.00000000},
 
942
      {0.00000000, 0.00000000, 0.00000000}};
 
943
      
 
944
      static const double dmats1[3][3] = \
 
945
      {{0.00000000, 0.00000000, 0.00000000},
 
946
      {2.44948974, 0.00000000, 0.00000000},
 
947
      {4.24264069, 0.00000000, 0.00000000}};
 
948
      
 
949
      // Compute reference derivatives.
 
950
      // Declare pointer to array of derivatives on FIAT element.
 
951
      double *derivatives = new double[num_derivatives];
 
952
      for (unsigned int r = 0; r < num_derivatives; r++)
 
953
      {
 
954
        derivatives[r] = 0.00000000;
 
955
      }// end loop over 'r'
 
956
      
 
957
      // Declare derivative matrix (of polynomial basis).
 
958
      double dmats[3][3] = \
 
959
      {{1.00000000, 0.00000000, 0.00000000},
 
960
      {0.00000000, 1.00000000, 0.00000000},
 
961
      {0.00000000, 0.00000000, 1.00000000}};
 
962
      
 
963
      // Declare (auxiliary) derivative matrix (of polynomial basis).
 
964
      double dmats_old[3][3] = \
 
965
      {{1.00000000, 0.00000000, 0.00000000},
 
966
      {0.00000000, 1.00000000, 0.00000000},
 
967
      {0.00000000, 0.00000000, 1.00000000}};
 
968
      
 
969
      // Loop possible derivatives.
 
970
      for (unsigned int r = 0; r < num_derivatives; r++)
 
971
      {
 
972
        // Resetting dmats values to compute next derivative.
 
973
        for (unsigned int t = 0; t < 3; t++)
 
974
        {
 
975
          for (unsigned int u = 0; u < 3; u++)
 
976
          {
 
977
            dmats[t][u] = 0.00000000;
 
978
            if (t == u)
 
979
            {
 
980
            dmats[t][u] = 1.00000000;
 
981
            }
 
982
            
 
983
          }// end loop over 'u'
 
984
        }// end loop over 't'
 
985
        
 
986
        // Looping derivative order to generate dmats.
 
987
        for (unsigned int s = 0; s < n; s++)
 
988
        {
 
989
          // Updating dmats_old with new values and resetting dmats.
 
990
          for (unsigned int t = 0; t < 3; t++)
 
991
          {
 
992
            for (unsigned int u = 0; u < 3; u++)
 
993
            {
 
994
              dmats_old[t][u] = dmats[t][u];
 
995
              dmats[t][u] = 0.00000000;
 
996
            }// end loop over 'u'
 
997
          }// end loop over 't'
 
998
          
 
999
          // Update dmats using an inner product.
 
1000
          if (combinations[r][s] == 0)
 
1001
          {
 
1002
          for (unsigned int t = 0; t < 3; t++)
 
1003
          {
 
1004
            for (unsigned int u = 0; u < 3; u++)
 
1005
            {
 
1006
              for (unsigned int tu = 0; tu < 3; tu++)
 
1007
              {
 
1008
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
1009
              }// end loop over 'tu'
 
1010
            }// end loop over 'u'
 
1011
          }// end loop over 't'
 
1012
          }
 
1013
          
 
1014
          if (combinations[r][s] == 1)
 
1015
          {
 
1016
          for (unsigned int t = 0; t < 3; t++)
 
1017
          {
 
1018
            for (unsigned int u = 0; u < 3; u++)
 
1019
            {
 
1020
              for (unsigned int tu = 0; tu < 3; tu++)
 
1021
              {
 
1022
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
1023
              }// end loop over 'tu'
 
1024
            }// end loop over 'u'
 
1025
          }// end loop over 't'
 
1026
          }
 
1027
          
 
1028
        }// end loop over 's'
 
1029
        for (unsigned int s = 0; s < 3; s++)
 
1030
        {
 
1031
          for (unsigned int t = 0; t < 3; t++)
 
1032
          {
 
1033
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
 
1034
          }// end loop over 't'
 
1035
        }// end loop over 's'
 
1036
      }// end loop over 'r'
 
1037
      
 
1038
      // Transform derivatives back to physical element
 
1039
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1040
      {
 
1041
        for (unsigned int s = 0; s < num_derivatives; s++)
 
1042
        {
 
1043
          values[r] += transform[r][s]*derivatives[s];
 
1044
        }// end loop over 's'
 
1045
      }// end loop over 'r'
 
1046
      
 
1047
      // Delete pointer to array of derivatives on FIAT element
 
1048
      delete [] derivatives;
 
1049
      
 
1050
      // Delete pointer to array of combinations of derivatives and transform
 
1051
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1052
      {
 
1053
        delete [] combinations[r];
 
1054
      }// end loop over 'r'
 
1055
      delete [] combinations;
 
1056
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1057
      {
 
1058
        delete [] transform[r];
 
1059
      }// end loop over 'r'
 
1060
      delete [] transform;
 
1061
        break;
 
1062
      }
 
1063
    case 2:
 
1064
      {
 
1065
        
 
1066
      // Array of basisvalues.
 
1067
      double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
 
1068
      
 
1069
      // Declare helper variables.
 
1070
      unsigned int rr = 0;
 
1071
      unsigned int ss = 0;
 
1072
      double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
 
1073
      
 
1074
      // Compute basisvalues.
 
1075
      basisvalues[0] = 1.00000000;
 
1076
      basisvalues[1] = tmp0;
 
1077
      for (unsigned int r = 0; r < 1; r++)
 
1078
      {
 
1079
        rr = (r + 1)*(r + 1 + 1)/2 + 1;
 
1080
        ss = r*(r + 1)/2;
 
1081
        basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
 
1082
      }// end loop over 'r'
 
1083
      for (unsigned int r = 0; r < 2; r++)
 
1084
      {
 
1085
        for (unsigned int s = 0; s < 2 - r; s++)
 
1086
        {
 
1087
          rr = (r + s)*(r + s + 1)/2 + s;
 
1088
          basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
 
1089
        }// end loop over 's'
 
1090
      }// end loop over 'r'
 
1091
      
 
1092
      // Table(s) of coefficients.
 
1093
      static const double coefficients0[3] = \
 
1094
      {0.47140452, 0.00000000, 0.33333333};
 
1095
      
 
1096
      // Tables of derivatives of the polynomial base (transpose).
 
1097
      static const double dmats0[3][3] = \
 
1098
      {{0.00000000, 0.00000000, 0.00000000},
 
1099
      {4.89897949, 0.00000000, 0.00000000},
 
1100
      {0.00000000, 0.00000000, 0.00000000}};
 
1101
      
 
1102
      static const double dmats1[3][3] = \
 
1103
      {{0.00000000, 0.00000000, 0.00000000},
 
1104
      {2.44948974, 0.00000000, 0.00000000},
 
1105
      {4.24264069, 0.00000000, 0.00000000}};
 
1106
      
 
1107
      // Compute reference derivatives.
 
1108
      // Declare pointer to array of derivatives on FIAT element.
 
1109
      double *derivatives = new double[num_derivatives];
 
1110
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1111
      {
 
1112
        derivatives[r] = 0.00000000;
 
1113
      }// end loop over 'r'
 
1114
      
 
1115
      // Declare derivative matrix (of polynomial basis).
 
1116
      double dmats[3][3] = \
 
1117
      {{1.00000000, 0.00000000, 0.00000000},
 
1118
      {0.00000000, 1.00000000, 0.00000000},
 
1119
      {0.00000000, 0.00000000, 1.00000000}};
 
1120
      
 
1121
      // Declare (auxiliary) derivative matrix (of polynomial basis).
 
1122
      double dmats_old[3][3] = \
 
1123
      {{1.00000000, 0.00000000, 0.00000000},
 
1124
      {0.00000000, 1.00000000, 0.00000000},
 
1125
      {0.00000000, 0.00000000, 1.00000000}};
 
1126
      
 
1127
      // Loop possible derivatives.
 
1128
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1129
      {
 
1130
        // Resetting dmats values to compute next derivative.
 
1131
        for (unsigned int t = 0; t < 3; t++)
 
1132
        {
 
1133
          for (unsigned int u = 0; u < 3; u++)
 
1134
          {
 
1135
            dmats[t][u] = 0.00000000;
 
1136
            if (t == u)
 
1137
            {
 
1138
            dmats[t][u] = 1.00000000;
 
1139
            }
 
1140
            
 
1141
          }// end loop over 'u'
 
1142
        }// end loop over 't'
 
1143
        
 
1144
        // Looping derivative order to generate dmats.
 
1145
        for (unsigned int s = 0; s < n; s++)
 
1146
        {
 
1147
          // Updating dmats_old with new values and resetting dmats.
 
1148
          for (unsigned int t = 0; t < 3; t++)
 
1149
          {
 
1150
            for (unsigned int u = 0; u < 3; u++)
 
1151
            {
 
1152
              dmats_old[t][u] = dmats[t][u];
 
1153
              dmats[t][u] = 0.00000000;
 
1154
            }// end loop over 'u'
 
1155
          }// end loop over 't'
 
1156
          
 
1157
          // Update dmats using an inner product.
 
1158
          if (combinations[r][s] == 0)
 
1159
          {
 
1160
          for (unsigned int t = 0; t < 3; t++)
 
1161
          {
 
1162
            for (unsigned int u = 0; u < 3; u++)
 
1163
            {
 
1164
              for (unsigned int tu = 0; tu < 3; tu++)
 
1165
              {
 
1166
                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
 
1167
              }// end loop over 'tu'
 
1168
            }// end loop over 'u'
 
1169
          }// end loop over 't'
 
1170
          }
 
1171
          
 
1172
          if (combinations[r][s] == 1)
 
1173
          {
 
1174
          for (unsigned int t = 0; t < 3; t++)
 
1175
          {
 
1176
            for (unsigned int u = 0; u < 3; u++)
 
1177
            {
 
1178
              for (unsigned int tu = 0; tu < 3; tu++)
 
1179
              {
 
1180
                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
 
1181
              }// end loop over 'tu'
 
1182
            }// end loop over 'u'
 
1183
          }// end loop over 't'
 
1184
          }
 
1185
          
 
1186
        }// end loop over 's'
 
1187
        for (unsigned int s = 0; s < 3; s++)
 
1188
        {
 
1189
          for (unsigned int t = 0; t < 3; t++)
 
1190
          {
 
1191
            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
 
1192
          }// end loop over 't'
 
1193
        }// end loop over 's'
 
1194
      }// end loop over 'r'
 
1195
      
 
1196
      // Transform derivatives back to physical element
 
1197
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1198
      {
 
1199
        for (unsigned int s = 0; s < num_derivatives; s++)
 
1200
        {
 
1201
          values[r] += transform[r][s]*derivatives[s];
 
1202
        }// end loop over 's'
 
1203
      }// end loop over 'r'
 
1204
      
 
1205
      // Delete pointer to array of derivatives on FIAT element
 
1206
      delete [] derivatives;
 
1207
      
 
1208
      // Delete pointer to array of combinations of derivatives and transform
 
1209
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1210
      {
 
1211
        delete [] combinations[r];
 
1212
      }// end loop over 'r'
 
1213
      delete [] combinations;
 
1214
      for (unsigned int r = 0; r < num_derivatives; r++)
 
1215
      {
 
1216
        delete [] transform[r];
 
1217
      }// end loop over 'r'
 
1218
      delete [] transform;
 
1219
        break;
 
1220
      }
 
1221
    }
 
1222
    
825
1223
  }
826
1224
 
827
1225
  /// Evaluate order n derivatives of all basis functions at given point in cell
1559
1957
    }// end loop over 'r'
1560
1958
    
1561
1959
    // Compute element tensor using UFL quadrature representation
1562
 
    // Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
 
1960
    // Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
1563
1961
    switch (facet)
1564
1962
    {
1565
1963
    case 0:
1743
2141
    }// end loop over 'r'
1744
2142
    
1745
2143
    // Compute element tensor using UFL quadrature representation
1746
 
    // Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
 
2144
    // Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
1747
2145
    switch (facet0)
1748
2146
    {
1749
2147
    case 0: