← Back to branch summary

~ubuntu-branches/ubuntu/trusty/ffc/trusty

« back to all changes in this revision

Viewing changes to test/regression/reference/tensor/ElementRestriction.h

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2010-02-03 20:22:35 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20100203202235-fe8d0kajuvgy2sqn
Tags: 0.9.0-1
* New upstream release.
* debian/control: Bump Standards-Version (no changes needed).
* Update debian/copyright and debian/copyright_hints.

Show diffs side-by-side

added added

removed removed

Lines of Context:
test/regression/reference/tensor/ElementRestriction.h
1
 
// This code conforms with the UFC specification version 1.0
2
 
// and was automatically generated by FFC version 0.7.0.
3
 
 
4
 
#ifndef __ELEMENTRESTRICTION_H
5
 
#define __ELEMENTRESTRICTION_H
6
 
 
7
 
#include <cmath>
8
 
#include <stdexcept>
9
 
#include <ufc.h>
10
 
 
11
 
/// This class defines the interface for a finite element.
12
 
 
13
 
class elementrestriction_0_finite_element_0_0: public ufc::finite_element
14
 
{
15
 
public:
16
 
 
17
 
  /// Constructor
18
 
  elementrestriction_0_finite_element_0_0() : ufc::finite_element()
19
 
  {
20
 
    // Do nothing
21
 
  }
22
 
 
23
 
  /// Destructor
24
 
  virtual ~elementrestriction_0_finite_element_0_0()
25
 
  {
26
 
    // Do nothing
27
 
  }
28
 
 
29
 
  /// Return a string identifying the finite element
30
 
  virtual const char* signature() const
31
 
  {
32
 
    return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 4)";
33
 
  }
34
 
 
35
 
  /// Return the cell shape
36
 
  virtual ufc::shape cell_shape() const
37
 
  {
38
 
    return ufc::triangle;
39
 
  }
40
 
 
41
 
  /// Return the dimension of the finite element function space
42
 
  virtual unsigned int space_dimension() const
43
 
  {
44
 
    return 12;
45
 
  }
46
 
 
47
 
  /// Return the rank of the value space
48
 
  virtual unsigned int value_rank() const
49
 
  {
50
 
    return 0;
51
 
  }
52
 
 
53
 
  /// Return the dimension of the value space for axis i
54
 
  virtual unsigned int value_dimension(unsigned int i) const
55
 
  {
56
 
    return 1;
57
 
  }
58
 
 
59
 
  /// Evaluate basis function i at given point in cell
60
 
  virtual void evaluate_basis(unsigned int i,
61
 
                              double* values,
62
 
                              const double* coordinates,
63
 
                              const ufc::cell& c) const
64
 
  {
65
 
    // Extract vertex coordinates
66
 
    const double * const * element_coordinates = c.coordinates;
67
 
    
68
 
    // Compute Jacobian of affine map from reference cell
69
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
70
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
71
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
72
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
73
 
    
74
 
    // Compute determinant of Jacobian
75
 
    const double detJ = J_00*J_11 - J_01*J_10;
76
 
    
77
 
    // Compute inverse of Jacobian
78
 
    
79
 
    // Get coordinates and map to the reference (UFC) element
80
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
81
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
82
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
83
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
84
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
85
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
86
 
    
87
 
    // Map coordinates to the reference square
88
 
    if (std::abs(y - 1.0) < 1e-08)
89
 
      x = -1.0;
90
 
    else
91
 
      x = 2.0 *x/(1.0 - y) - 1.0;
92
 
    y = 2.0*y - 1.0;
93
 
    
94
 
    // Reset values
95
 
    *values = 0;
96
 
    
97
 
    // Map degree of freedom to element degree of freedom
98
 
    const unsigned int dof = i;
99
 
    
100
 
    // Generate scalings
101
 
    const double scalings_y_0 = 1;
102
 
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
103
 
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
104
 
    const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
105
 
    const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
106
 
    
107
 
    // Compute psitilde_a
108
 
    const double psitilde_a_0 = 1;
109
 
    const double psitilde_a_1 = x;
110
 
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
111
 
    const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;
112
 
    const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
113
 
    
114
 
    // Compute psitilde_bs
115
 
    const double psitilde_bs_0_0 = 1;
116
 
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
117
 
    const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
118
 
    const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
119
 
    const double psitilde_bs_0_4 = 0.0285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571*psitilde_bs_0_2;
120
 
    const double psitilde_bs_1_0 = 1;
121
 
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
122
 
    const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
123
 
    const double psitilde_bs_1_3 = 0.285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714*psitilde_bs_1_1;
124
 
    const double psitilde_bs_2_0 = 1;
125
 
    const double psitilde_bs_2_1 = 3.5*y + 2.5;
126
 
    const double psitilde_bs_2_2 = 1.02040816*psitilde_bs_2_1 + 2.57142857*y*psitilde_bs_2_1 - 0.551020408*psitilde_bs_2_0;
127
 
    const double psitilde_bs_3_0 = 1;
128
 
    const double psitilde_bs_3_1 = 4.5*y + 3.5;
129
 
    const double psitilde_bs_4_0 = 1;
130
 
    
131
 
    // Compute basisvalues
132
 
    const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
133
 
    const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
134
 
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
135
 
    const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
136
 
    const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
137
 
    const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
138
 
    const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
139
 
    const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
140
 
    const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
141
 
    const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
142
 
    const double basisvalue10 = 4.74341649*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
143
 
    const double basisvalue11 = 4.18330013*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
144
 
    const double basisvalue12 = 3.53553391*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
145
 
    const double basisvalue13 = 2.73861279*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
146
 
    const double basisvalue14 = 1.58113883*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
147
 
    
148
 
    // Table(s) of coefficients
149
 
    static const double coefficients0[12][15] = \
150
 
    {{0, -0.0412393049, -0.0238095238, 0.0289800295, 0.0224478343, 0.0129602632, -0.0395942581, -0.0334632557, -0.0259205264, -0.0149652229, 0.0321247254, 0.0283313448, 0.0239443566, 0.0185472189, 0.0107082418},
151
 
    {0, 0.0412393049, -0.0238095238, 0.0289800295, -0.0224478343, 0.0129602632, 0.0395942581, -0.0334632557, 0.0259205264, -0.0149652229, 0.0321247254, -0.0283313448, 0.0239443566, -0.0185472189, 0.0107082418},
152
 
    {0, 0, 0.0476190476, 0, 0, 0.0388807896, 0, 0, 0, 0.0598608915, 0, 0, 0, 0, 0.0535412091},
153
 
    {0.125707872, 0.131965776, -0.0253968254, 0.139104142, -0.0718330698, 0.0311046317, 0.0633508129, 0.0267706045, -0.0622092633, 0.0478887132, 0, 0.0566626896, -0.0838052481, 0.083462485, -0.0535412091},
154
 
    {-0.0314269681, 0.010997148, 0.00634920635, 0, 0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, -0.139104142, 0.107082418},
155
 
    {0.125707872, 0.0439885919, 0.126984127, 0, 0.0359165349, 0.155523158, 0, 0, 0.103682106, -0.0119721783, 0, 0, 0, 0.0927360944, -0.107082418},
156
 
    {0.125707872, -0.131965776, -0.0253968254, 0.139104142, 0.0718330698, 0.0311046317, -0.0633508129, 0.0267706045, 0.0622092633, 0.0478887132, 0, -0.0566626896, -0.0838052481, -0.083462485, -0.0535412091},
157
 
    {-0.0314269681, -0.010997148, 0.00634920635, 0, -0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, 0.139104142, 0.107082418},
158
 
    {0.125707872, -0.0439885919, 0.126984127, 0, -0.0359165349, 0.155523158, 0, 0, -0.103682106, -0.0119721783, 0, 0, 0, -0.0927360944, -0.107082418},
159
 
    {0.125707872, -0.0879771839, -0.101587302, 0.0927360944, 0.107749605, 0.0725774739, 0.0791885161, -0.0133853023, -0.0518410528, -0.0419026241, -0.128498902, -0.0566626896, -0.0119721783, 0.00927360944, 0.0107082418},
160
 
    {-0.0314269681, 0, -0.0126984127, -0.243432248, 0, 0.0544331054, 0, 0.0936971159, 0, -0.0419026241, 0.192748353, 0, -0.0239443566, 0, 0.0107082418},
161
 
    {0.125707872, 0.0879771839, -0.101587302, 0.0927360944, -0.107749605, 0.0725774739, -0.0791885161, -0.0133853023, 0.0518410528, -0.0419026241, -0.128498902, 0.0566626896, -0.0119721783, -0.00927360944, 0.0107082418}};
162
 
    
163
 
    // Extract relevant coefficients
164
 
    const double coeff0_0 = coefficients0[dof][0];
165
 
    const double coeff0_1 = coefficients0[dof][1];
166
 
    const double coeff0_2 = coefficients0[dof][2];
167
 
    const double coeff0_3 = coefficients0[dof][3];
168
 
    const double coeff0_4 = coefficients0[dof][4];
169
 
    const double coeff0_5 = coefficients0[dof][5];
170
 
    const double coeff0_6 = coefficients0[dof][6];
171
 
    const double coeff0_7 = coefficients0[dof][7];
172
 
    const double coeff0_8 = coefficients0[dof][8];
173
 
    const double coeff0_9 = coefficients0[dof][9];
174
 
    const double coeff0_10 = coefficients0[dof][10];
175
 
    const double coeff0_11 = coefficients0[dof][11];
176
 
    const double coeff0_12 = coefficients0[dof][12];
177
 
    const double coeff0_13 = coefficients0[dof][13];
178
 
    const double coeff0_14 = coefficients0[dof][14];
179
 
    
180
 
    // Compute value(s)
181
 
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5 + coeff0_6*basisvalue6 + coeff0_7*basisvalue7 + coeff0_8*basisvalue8 + coeff0_9*basisvalue9 + coeff0_10*basisvalue10 + coeff0_11*basisvalue11 + coeff0_12*basisvalue12 + coeff0_13*basisvalue13 + coeff0_14*basisvalue14;
182
 
  }
183
 
 
184
 
  /// Evaluate all basis functions at given point in cell
185
 
  virtual void evaluate_basis_all(double* values,
186
 
                                  const double* coordinates,
187
 
                                  const ufc::cell& c) const
188
 
  {
189
 
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
190
 
  }
191
 
 
192
 
  /// Evaluate order n derivatives of basis function i at given point in cell
193
 
  virtual void evaluate_basis_derivatives(unsigned int i,
194
 
                                          unsigned int n,
195
 
                                          double* values,
196
 
                                          const double* coordinates,
197
 
                                          const ufc::cell& c) const
198
 
  {
199
 
    // Extract vertex coordinates
200
 
    const double * const * element_coordinates = c.coordinates;
201
 
    
202
 
    // Compute Jacobian of affine map from reference cell
203
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
204
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
205
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
206
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
207
 
    
208
 
    // Compute determinant of Jacobian
209
 
    const double detJ = J_00*J_11 - J_01*J_10;
210
 
    
211
 
    // Compute inverse of Jacobian
212
 
    
213
 
    // Get coordinates and map to the reference (UFC) element
214
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
215
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
216
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
217
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
218
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
219
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
220
 
    
221
 
    // Map coordinates to the reference square
222
 
    if (std::abs(y - 1.0) < 1e-08)
223
 
      x = -1.0;
224
 
    else
225
 
      x = 2.0 *x/(1.0 - y) - 1.0;
226
 
    y = 2.0*y - 1.0;
227
 
    
228
 
    // Compute number of derivatives
229
 
    unsigned int num_derivatives = 1;
230
 
    
231
 
    for (unsigned int j = 0; j < n; j++)
232
 
      num_derivatives *= 2;
233
 
    
234
 
    
235
 
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
236
 
    unsigned int **combinations = new unsigned int *[num_derivatives];
237
 
    
238
 
    for (unsigned int j = 0; j < num_derivatives; j++)
239
 
    {
240
 
      combinations[j] = new unsigned int [n];
241
 
      for (unsigned int k = 0; k < n; k++)
242
 
        combinations[j][k] = 0;
243
 
    }
244
 
    
245
 
    // Generate combinations of derivatives
246
 
    for (unsigned int row = 1; row < num_derivatives; row++)
247
 
    {
248
 
      for (unsigned int num = 0; num < row; num++)
249
 
      {
250
 
        for (unsigned int col = n-1; col+1 > 0; col--)
251
 
        {
252
 
          if (combinations[row][col] + 1 > 1)
253
 
            combinations[row][col] = 0;
254
 
          else
255
 
          {
256
 
            combinations[row][col] += 1;
257
 
            break;
258
 
          }
259
 
        }
260
 
      }
261
 
    }
262
 
    
263
 
    // Compute inverse of Jacobian
264
 
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
265
 
    
266
 
    // Declare transformation matrix
267
 
    // Declare pointer to two dimensional array and initialise
268
 
    double **transform = new double *[num_derivatives];
269
 
    
270
 
    for (unsigned int j = 0; j < num_derivatives; j++)
271
 
    {
272
 
      transform[j] = new double [num_derivatives];
273
 
      for (unsigned int k = 0; k < num_derivatives; k++)
274
 
        transform[j][k] = 1;
275
 
    }
276
 
    
277
 
    // Construct transformation matrix
278
 
    for (unsigned int row = 0; row < num_derivatives; row++)
279
 
    {
280
 
      for (unsigned int col = 0; col < num_derivatives; col++)
281
 
      {
282
 
        for (unsigned int k = 0; k < n; k++)
283
 
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
284
 
      }
285
 
    }
286
 
    
287
 
    // Reset values
288
 
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
289
 
      values[j] = 0;
290
 
    
291
 
    // Map degree of freedom to element degree of freedom
292
 
    const unsigned int dof = i;
293
 
    
294
 
    // Generate scalings
295
 
    const double scalings_y_0 = 1;
296
 
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
297
 
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
298
 
    const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
299
 
    const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
300
 
    
301
 
    // Compute psitilde_a
302
 
    const double psitilde_a_0 = 1;
303
 
    const double psitilde_a_1 = x;
304
 
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
305
 
    const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;
306
 
    const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
307
 
    
308
 
    // Compute psitilde_bs
309
 
    const double psitilde_bs_0_0 = 1;
310
 
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
311
 
    const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
312
 
    const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
313
 
    const double psitilde_bs_0_4 = 0.0285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571*psitilde_bs_0_2;
314
 
    const double psitilde_bs_1_0 = 1;
315
 
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
316
 
    const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
317
 
    const double psitilde_bs_1_3 = 0.285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714*psitilde_bs_1_1;
318
 
    const double psitilde_bs_2_0 = 1;
319
 
    const double psitilde_bs_2_1 = 3.5*y + 2.5;
320
 
    const double psitilde_bs_2_2 = 1.02040816*psitilde_bs_2_1 + 2.57142857*y*psitilde_bs_2_1 - 0.551020408*psitilde_bs_2_0;
321
 
    const double psitilde_bs_3_0 = 1;
322
 
    const double psitilde_bs_3_1 = 4.5*y + 3.5;
323
 
    const double psitilde_bs_4_0 = 1;
324
 
    
325
 
    // Compute basisvalues
326
 
    const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
327
 
    const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
328
 
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
329
 
    const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
330
 
    const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
331
 
    const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
332
 
    const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
333
 
    const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
334
 
    const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
335
 
    const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
336
 
    const double basisvalue10 = 4.74341649*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
337
 
    const double basisvalue11 = 4.18330013*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
338
 
    const double basisvalue12 = 3.53553391*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
339
 
    const double basisvalue13 = 2.73861279*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
340
 
    const double basisvalue14 = 1.58113883*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
341
 
    
342
 
    // Table(s) of coefficients
343
 
    static const double coefficients0[12][15] = \
344
 
    {{0, -0.0412393049, -0.0238095238, 0.0289800295, 0.0224478343, 0.0129602632, -0.0395942581, -0.0334632557, -0.0259205264, -0.0149652229, 0.0321247254, 0.0283313448, 0.0239443566, 0.0185472189, 0.0107082418},
345
 
    {0, 0.0412393049, -0.0238095238, 0.0289800295, -0.0224478343, 0.0129602632, 0.0395942581, -0.0334632557, 0.0259205264, -0.0149652229, 0.0321247254, -0.0283313448, 0.0239443566, -0.0185472189, 0.0107082418},
346
 
    {0, 0, 0.0476190476, 0, 0, 0.0388807896, 0, 0, 0, 0.0598608915, 0, 0, 0, 0, 0.0535412091},
347
 
    {0.125707872, 0.131965776, -0.0253968254, 0.139104142, -0.0718330698, 0.0311046317, 0.0633508129, 0.0267706045, -0.0622092633, 0.0478887132, 0, 0.0566626896, -0.0838052481, 0.083462485, -0.0535412091},
348
 
    {-0.0314269681, 0.010997148, 0.00634920635, 0, 0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, -0.139104142, 0.107082418},
349
 
    {0.125707872, 0.0439885919, 0.126984127, 0, 0.0359165349, 0.155523158, 0, 0, 0.103682106, -0.0119721783, 0, 0, 0, 0.0927360944, -0.107082418},
350
 
    {0.125707872, -0.131965776, -0.0253968254, 0.139104142, 0.0718330698, 0.0311046317, -0.0633508129, 0.0267706045, 0.0622092633, 0.0478887132, 0, -0.0566626896, -0.0838052481, -0.083462485, -0.0535412091},
351
 
    {-0.0314269681, -0.010997148, 0.00634920635, 0, -0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, 0.139104142, 0.107082418},
352
 
    {0.125707872, -0.0439885919, 0.126984127, 0, -0.0359165349, 0.155523158, 0, 0, -0.103682106, -0.0119721783, 0, 0, 0, -0.0927360944, -0.107082418},
353
 
    {0.125707872, -0.0879771839, -0.101587302, 0.0927360944, 0.107749605, 0.0725774739, 0.0791885161, -0.0133853023, -0.0518410528, -0.0419026241, -0.128498902, -0.0566626896, -0.0119721783, 0.00927360944, 0.0107082418},
354
 
    {-0.0314269681, 0, -0.0126984127, -0.243432248, 0, 0.0544331054, 0, 0.0936971159, 0, -0.0419026241, 0.192748353, 0, -0.0239443566, 0, 0.0107082418},
355
 
    {0.125707872, 0.0879771839, -0.101587302, 0.0927360944, -0.107749605, 0.0725774739, -0.0791885161, -0.0133853023, 0.0518410528, -0.0419026241, -0.128498902, 0.0566626896, -0.0119721783, -0.00927360944, 0.0107082418}};
356
 
    
357
 
    // Interesting (new) part
358
 
    // Tables of derivatives of the polynomial base (transpose)
359
 
    static const double dmats0[15][15] = \
360
 
    {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
361
 
    {4.89897949, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
362
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
363
 
    {0, 9.48683298, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
364
 
    {4, 0, 7.07106781, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
365
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
366
 
    {5.29150262, 0, -2.99332591, 13.662601, 0, 0.611010093, 0, 0, 0, 0, 0, 0, 0, 0, 0},
367
 
    {0, 4.38178046, 0, 0, 12.5219807, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
368
 
    {3.46410162, 0, 7.83836718, 0, 0, 8.4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
369
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
370
 
    {0, 10.9544512, 0, 0, -3.83325939, 0, 17.7482393, 0, 0.553283335, 0, 0, 0, 0, 0, 0},
371
 
    {4.73286383, 0, 3.34664011, 4.3643578, 0, -5.07468038, 0, 17.0084013, 0, 1.52127766, 0, 0, 0, 0, 0},
372
 
    {0, 2.44948974, 0, 0, 9.14285714, 0, 0, 0, 14.8461498, 0, 0, 0, 0, 0, 0},
373
 
    {3.09838668, 0, 7.66811581, 0, 0, 10.7331263, 0, 0, 0, 9.29516003, 0, 0, 0, 0, 0},
374
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
375
 
    
376
 
    static const double dmats1[15][15] = \
377
 
    {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
378
 
    {2.44948974, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
379
 
    {4.24264069, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
380
 
    {2.5819889, 4.74341649, -0.912870929, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
381
 
    {2, 6.12372436, 3.53553391, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
382
 
    {-2.30940108, 0, 8.16496581, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
383
 
    {2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.305505046, 0, 0, 0, 0, 0, 0, 0, 0, 0},
384
 
    {2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0, 0, 0, 0, 0, 0, 0, 0, 0},
385
 
    {1.73205081, -5.09116882, 3.91918359, 0, 9.69948452, 4.2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
386
 
    {5, 0, -2.82842712, 0, 0, 12.1243557, 0, 0, 0, 0, 0, 0, 0, 0, 0},
387
 
    {2.68328157, 5.47722558, -1.8973666, 7.42307489, -1.91662969, 0.663940002, 8.87411967, -1.07142857, 0.276641668, -0.0958314847, 0, 0, 0, 0, 0},
388
 
    {2.36643191, 2.89827535, 1.67332005, 2.1821789, 5.74704893, -2.53734019, 10.0623059, 8.50420064, -2.19577516, 0.760638829, 0, 0, 0, 0, 0},
389
 
    {2, 1.22474487, 3.53553391, -7.37711114, 4.57142857, 1.6495722, 0, 11.4997782, 7.42307489, -2.57142857, 0, 0, 0, 0, 0},
390
 
    {1.54919334, 6.64078309, 3.8340579, 0, -6.19677335, 5.36656315, 0, 0, 13.4164079, 4.64758002, 0, 0, 0, 0, 0},
391
 
    {-3.57770876, 0, 8.85437745, 0, 0, -3.09838668, 0, 0, 0, 16.0996894, 0, 0, 0, 0, 0}};
392
 
    
393
 
    // Compute reference derivatives
394
 
    // Declare pointer to array of derivatives on FIAT element
395
 
    double *derivatives = new double [num_derivatives];
396
 
    
397
 
    // Declare coefficients
398
 
    double coeff0_0 = 0;
399
 
    double coeff0_1 = 0;
400
 
    double coeff0_2 = 0;
401
 
    double coeff0_3 = 0;
402
 
    double coeff0_4 = 0;
403
 
    double coeff0_5 = 0;
404
 
    double coeff0_6 = 0;
405
 
    double coeff0_7 = 0;
406
 
    double coeff0_8 = 0;
407
 
    double coeff0_9 = 0;
408
 
    double coeff0_10 = 0;
409
 
    double coeff0_11 = 0;
410
 
    double coeff0_12 = 0;
411
 
    double coeff0_13 = 0;
412
 
    double coeff0_14 = 0;
413
 
    
414
 
    // Declare new coefficients
415
 
    double new_coeff0_0 = 0;
416
 
    double new_coeff0_1 = 0;
417
 
    double new_coeff0_2 = 0;
418
 
    double new_coeff0_3 = 0;
419
 
    double new_coeff0_4 = 0;
420
 
    double new_coeff0_5 = 0;
421
 
    double new_coeff0_6 = 0;
422
 
    double new_coeff0_7 = 0;
423
 
    double new_coeff0_8 = 0;
424
 
    double new_coeff0_9 = 0;
425
 
    double new_coeff0_10 = 0;
426
 
    double new_coeff0_11 = 0;
427
 
    double new_coeff0_12 = 0;
428
 
    double new_coeff0_13 = 0;
429
 
    double new_coeff0_14 = 0;
430
 
    
431
 
    // Loop possible derivatives
432
 
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
433
 
    {
434
 
      // Get values from coefficients array
435
 
      new_coeff0_0 = coefficients0[dof][0];
436
 
      new_coeff0_1 = coefficients0[dof][1];
437
 
      new_coeff0_2 = coefficients0[dof][2];
438
 
      new_coeff0_3 = coefficients0[dof][3];
439
 
      new_coeff0_4 = coefficients0[dof][4];
440
 
      new_coeff0_5 = coefficients0[dof][5];
441
 
      new_coeff0_6 = coefficients0[dof][6];
442
 
      new_coeff0_7 = coefficients0[dof][7];
443
 
      new_coeff0_8 = coefficients0[dof][8];
444
 
      new_coeff0_9 = coefficients0[dof][9];
445
 
      new_coeff0_10 = coefficients0[dof][10];
446
 
      new_coeff0_11 = coefficients0[dof][11];
447
 
      new_coeff0_12 = coefficients0[dof][12];
448
 
      new_coeff0_13 = coefficients0[dof][13];
449
 
      new_coeff0_14 = coefficients0[dof][14];
450
 
    
451
 
      // Loop derivative order
452
 
      for (unsigned int j = 0; j < n; j++)
453
 
      {
454
 
        // Update old coefficients
455
 
        coeff0_0 = new_coeff0_0;
456
 
        coeff0_1 = new_coeff0_1;
457
 
        coeff0_2 = new_coeff0_2;
458
 
        coeff0_3 = new_coeff0_3;
459
 
        coeff0_4 = new_coeff0_4;
460
 
        coeff0_5 = new_coeff0_5;
461
 
        coeff0_6 = new_coeff0_6;
462
 
        coeff0_7 = new_coeff0_7;
463
 
        coeff0_8 = new_coeff0_8;
464
 
        coeff0_9 = new_coeff0_9;
465
 
        coeff0_10 = new_coeff0_10;
466
 
        coeff0_11 = new_coeff0_11;
467
 
        coeff0_12 = new_coeff0_12;
468
 
        coeff0_13 = new_coeff0_13;
469
 
        coeff0_14 = new_coeff0_14;
470
 
    
471
 
        if(combinations[deriv_num][j] == 0)
472
 
        {
473
 
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0] + coeff0_6*dmats0[6][0] + coeff0_7*dmats0[7][0] + coeff0_8*dmats0[8][0] + coeff0_9*dmats0[9][0] + coeff0_10*dmats0[10][0] + coeff0_11*dmats0[11][0] + coeff0_12*dmats0[12][0] + coeff0_13*dmats0[13][0] + coeff0_14*dmats0[14][0];
474
 
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1] + coeff0_6*dmats0[6][1] + coeff0_7*dmats0[7][1] + coeff0_8*dmats0[8][1] + coeff0_9*dmats0[9][1] + coeff0_10*dmats0[10][1] + coeff0_11*dmats0[11][1] + coeff0_12*dmats0[12][1] + coeff0_13*dmats0[13][1] + coeff0_14*dmats0[14][1];
475
 
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2] + coeff0_6*dmats0[6][2] + coeff0_7*dmats0[7][2] + coeff0_8*dmats0[8][2] + coeff0_9*dmats0[9][2] + coeff0_10*dmats0[10][2] + coeff0_11*dmats0[11][2] + coeff0_12*dmats0[12][2] + coeff0_13*dmats0[13][2] + coeff0_14*dmats0[14][2];
476
 
          new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3] + coeff0_6*dmats0[6][3] + coeff0_7*dmats0[7][3] + coeff0_8*dmats0[8][3] + coeff0_9*dmats0[9][3] + coeff0_10*dmats0[10][3] + coeff0_11*dmats0[11][3] + coeff0_12*dmats0[12][3] + coeff0_13*dmats0[13][3] + coeff0_14*dmats0[14][3];
477
 
          new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4] + coeff0_6*dmats0[6][4] + coeff0_7*dmats0[7][4] + coeff0_8*dmats0[8][4] + coeff0_9*dmats0[9][4] + coeff0_10*dmats0[10][4] + coeff0_11*dmats0[11][4] + coeff0_12*dmats0[12][4] + coeff0_13*dmats0[13][4] + coeff0_14*dmats0[14][4];
478
 
          new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5] + coeff0_6*dmats0[6][5] + coeff0_7*dmats0[7][5] + coeff0_8*dmats0[8][5] + coeff0_9*dmats0[9][5] + coeff0_10*dmats0[10][5] + coeff0_11*dmats0[11][5] + coeff0_12*dmats0[12][5] + coeff0_13*dmats0[13][5] + coeff0_14*dmats0[14][5];
479
 
          new_coeff0_6 = coeff0_0*dmats0[0][6] + coeff0_1*dmats0[1][6] + coeff0_2*dmats0[2][6] + coeff0_3*dmats0[3][6] + coeff0_4*dmats0[4][6] + coeff0_5*dmats0[5][6] + coeff0_6*dmats0[6][6] + coeff0_7*dmats0[7][6] + coeff0_8*dmats0[8][6] + coeff0_9*dmats0[9][6] + coeff0_10*dmats0[10][6] + coeff0_11*dmats0[11][6] + coeff0_12*dmats0[12][6] + coeff0_13*dmats0[13][6] + coeff0_14*dmats0[14][6];
480
 
          new_coeff0_7 = coeff0_0*dmats0[0][7] + coeff0_1*dmats0[1][7] + coeff0_2*dmats0[2][7] + coeff0_3*dmats0[3][7] + coeff0_4*dmats0[4][7] + coeff0_5*dmats0[5][7] + coeff0_6*dmats0[6][7] + coeff0_7*dmats0[7][7] + coeff0_8*dmats0[8][7] + coeff0_9*dmats0[9][7] + coeff0_10*dmats0[10][7] + coeff0_11*dmats0[11][7] + coeff0_12*dmats0[12][7] + coeff0_13*dmats0[13][7] + coeff0_14*dmats0[14][7];
481
 
          new_coeff0_8 = coeff0_0*dmats0[0][8] + coeff0_1*dmats0[1][8] + coeff0_2*dmats0[2][8] + coeff0_3*dmats0[3][8] + coeff0_4*dmats0[4][8] + coeff0_5*dmats0[5][8] + coeff0_6*dmats0[6][8] + coeff0_7*dmats0[7][8] + coeff0_8*dmats0[8][8] + coeff0_9*dmats0[9][8] + coeff0_10*dmats0[10][8] + coeff0_11*dmats0[11][8] + coeff0_12*dmats0[12][8] + coeff0_13*dmats0[13][8] + coeff0_14*dmats0[14][8];
482
 
          new_coeff0_9 = coeff0_0*dmats0[0][9] + coeff0_1*dmats0[1][9] + coeff0_2*dmats0[2][9] + coeff0_3*dmats0[3][9] + coeff0_4*dmats0[4][9] + coeff0_5*dmats0[5][9] + coeff0_6*dmats0[6][9] + coeff0_7*dmats0[7][9] + coeff0_8*dmats0[8][9] + coeff0_9*dmats0[9][9] + coeff0_10*dmats0[10][9] + coeff0_11*dmats0[11][9] + coeff0_12*dmats0[12][9] + coeff0_13*dmats0[13][9] + coeff0_14*dmats0[14][9];
483
 
          new_coeff0_10 = coeff0_0*dmats0[0][10] + coeff0_1*dmats0[1][10] + coeff0_2*dmats0[2][10] + coeff0_3*dmats0[3][10] + coeff0_4*dmats0[4][10] + coeff0_5*dmats0[5][10] + coeff0_6*dmats0[6][10] + coeff0_7*dmats0[7][10] + coeff0_8*dmats0[8][10] + coeff0_9*dmats0[9][10] + coeff0_10*dmats0[10][10] + coeff0_11*dmats0[11][10] + coeff0_12*dmats0[12][10] + coeff0_13*dmats0[13][10] + coeff0_14*dmats0[14][10];
484
 
          new_coeff0_11 = coeff0_0*dmats0[0][11] + coeff0_1*dmats0[1][11] + coeff0_2*dmats0[2][11] + coeff0_3*dmats0[3][11] + coeff0_4*dmats0[4][11] + coeff0_5*dmats0[5][11] + coeff0_6*dmats0[6][11] + coeff0_7*dmats0[7][11] + coeff0_8*dmats0[8][11] + coeff0_9*dmats0[9][11] + coeff0_10*dmats0[10][11] + coeff0_11*dmats0[11][11] + coeff0_12*dmats0[12][11] + coeff0_13*dmats0[13][11] + coeff0_14*dmats0[14][11];
485
 
          new_coeff0_12 = coeff0_0*dmats0[0][12] + coeff0_1*dmats0[1][12] + coeff0_2*dmats0[2][12] + coeff0_3*dmats0[3][12] + coeff0_4*dmats0[4][12] + coeff0_5*dmats0[5][12] + coeff0_6*dmats0[6][12] + coeff0_7*dmats0[7][12] + coeff0_8*dmats0[8][12] + coeff0_9*dmats0[9][12] + coeff0_10*dmats0[10][12] + coeff0_11*dmats0[11][12] + coeff0_12*dmats0[12][12] + coeff0_13*dmats0[13][12] + coeff0_14*dmats0[14][12];
486
 
          new_coeff0_13 = coeff0_0*dmats0[0][13] + coeff0_1*dmats0[1][13] + coeff0_2*dmats0[2][13] + coeff0_3*dmats0[3][13] + coeff0_4*dmats0[4][13] + coeff0_5*dmats0[5][13] + coeff0_6*dmats0[6][13] + coeff0_7*dmats0[7][13] + coeff0_8*dmats0[8][13] + coeff0_9*dmats0[9][13] + coeff0_10*dmats0[10][13] + coeff0_11*dmats0[11][13] + coeff0_12*dmats0[12][13] + coeff0_13*dmats0[13][13] + coeff0_14*dmats0[14][13];
487
 
          new_coeff0_14 = coeff0_0*dmats0[0][14] + coeff0_1*dmats0[1][14] + coeff0_2*dmats0[2][14] + coeff0_3*dmats0[3][14] + coeff0_4*dmats0[4][14] + coeff0_5*dmats0[5][14] + coeff0_6*dmats0[6][14] + coeff0_7*dmats0[7][14] + coeff0_8*dmats0[8][14] + coeff0_9*dmats0[9][14] + coeff0_10*dmats0[10][14] + coeff0_11*dmats0[11][14] + coeff0_12*dmats0[12][14] + coeff0_13*dmats0[13][14] + coeff0_14*dmats0[14][14];
488
 
        }
489
 
        if(combinations[deriv_num][j] == 1)
490
 
        {
491
 
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0] + coeff0_6*dmats1[6][0] + coeff0_7*dmats1[7][0] + coeff0_8*dmats1[8][0] + coeff0_9*dmats1[9][0] + coeff0_10*dmats1[10][0] + coeff0_11*dmats1[11][0] + coeff0_12*dmats1[12][0] + coeff0_13*dmats1[13][0] + coeff0_14*dmats1[14][0];
492
 
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1] + coeff0_6*dmats1[6][1] + coeff0_7*dmats1[7][1] + coeff0_8*dmats1[8][1] + coeff0_9*dmats1[9][1] + coeff0_10*dmats1[10][1] + coeff0_11*dmats1[11][1] + coeff0_12*dmats1[12][1] + coeff0_13*dmats1[13][1] + coeff0_14*dmats1[14][1];
493
 
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2] + coeff0_6*dmats1[6][2] + coeff0_7*dmats1[7][2] + coeff0_8*dmats1[8][2] + coeff0_9*dmats1[9][2] + coeff0_10*dmats1[10][2] + coeff0_11*dmats1[11][2] + coeff0_12*dmats1[12][2] + coeff0_13*dmats1[13][2] + coeff0_14*dmats1[14][2];
494
 
          new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3] + coeff0_6*dmats1[6][3] + coeff0_7*dmats1[7][3] + coeff0_8*dmats1[8][3] + coeff0_9*dmats1[9][3] + coeff0_10*dmats1[10][3] + coeff0_11*dmats1[11][3] + coeff0_12*dmats1[12][3] + coeff0_13*dmats1[13][3] + coeff0_14*dmats1[14][3];
495
 
          new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4] + coeff0_6*dmats1[6][4] + coeff0_7*dmats1[7][4] + coeff0_8*dmats1[8][4] + coeff0_9*dmats1[9][4] + coeff0_10*dmats1[10][4] + coeff0_11*dmats1[11][4] + coeff0_12*dmats1[12][4] + coeff0_13*dmats1[13][4] + coeff0_14*dmats1[14][4];
496
 
          new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5] + coeff0_6*dmats1[6][5] + coeff0_7*dmats1[7][5] + coeff0_8*dmats1[8][5] + coeff0_9*dmats1[9][5] + coeff0_10*dmats1[10][5] + coeff0_11*dmats1[11][5] + coeff0_12*dmats1[12][5] + coeff0_13*dmats1[13][5] + coeff0_14*dmats1[14][5];
497
 
          new_coeff0_6 = coeff0_0*dmats1[0][6] + coeff0_1*dmats1[1][6] + coeff0_2*dmats1[2][6] + coeff0_3*dmats1[3][6] + coeff0_4*dmats1[4][6] + coeff0_5*dmats1[5][6] + coeff0_6*dmats1[6][6] + coeff0_7*dmats1[7][6] + coeff0_8*dmats1[8][6] + coeff0_9*dmats1[9][6] + coeff0_10*dmats1[10][6] + coeff0_11*dmats1[11][6] + coeff0_12*dmats1[12][6] + coeff0_13*dmats1[13][6] + coeff0_14*dmats1[14][6];
498
 
          new_coeff0_7 = coeff0_0*dmats1[0][7] + coeff0_1*dmats1[1][7] + coeff0_2*dmats1[2][7] + coeff0_3*dmats1[3][7] + coeff0_4*dmats1[4][7] + coeff0_5*dmats1[5][7] + coeff0_6*dmats1[6][7] + coeff0_7*dmats1[7][7] + coeff0_8*dmats1[8][7] + coeff0_9*dmats1[9][7] + coeff0_10*dmats1[10][7] + coeff0_11*dmats1[11][7] + coeff0_12*dmats1[12][7] + coeff0_13*dmats1[13][7] + coeff0_14*dmats1[14][7];
499
 
          new_coeff0_8 = coeff0_0*dmats1[0][8] + coeff0_1*dmats1[1][8] + coeff0_2*dmats1[2][8] + coeff0_3*dmats1[3][8] + coeff0_4*dmats1[4][8] + coeff0_5*dmats1[5][8] + coeff0_6*dmats1[6][8] + coeff0_7*dmats1[7][8] + coeff0_8*dmats1[8][8] + coeff0_9*dmats1[9][8] + coeff0_10*dmats1[10][8] + coeff0_11*dmats1[11][8] + coeff0_12*dmats1[12][8] + coeff0_13*dmats1[13][8] + coeff0_14*dmats1[14][8];
500
 
          new_coeff0_9 = coeff0_0*dmats1[0][9] + coeff0_1*dmats1[1][9] + coeff0_2*dmats1[2][9] + coeff0_3*dmats1[3][9] + coeff0_4*dmats1[4][9] + coeff0_5*dmats1[5][9] + coeff0_6*dmats1[6][9] + coeff0_7*dmats1[7][9] + coeff0_8*dmats1[8][9] + coeff0_9*dmats1[9][9] + coeff0_10*dmats1[10][9] + coeff0_11*dmats1[11][9] + coeff0_12*dmats1[12][9] + coeff0_13*dmats1[13][9] + coeff0_14*dmats1[14][9];
501
 
          new_coeff0_10 = coeff0_0*dmats1[0][10] + coeff0_1*dmats1[1][10] + coeff0_2*dmats1[2][10] + coeff0_3*dmats1[3][10] + coeff0_4*dmats1[4][10] + coeff0_5*dmats1[5][10] + coeff0_6*dmats1[6][10] + coeff0_7*dmats1[7][10] + coeff0_8*dmats1[8][10] + coeff0_9*dmats1[9][10] + coeff0_10*dmats1[10][10] + coeff0_11*dmats1[11][10] + coeff0_12*dmats1[12][10] + coeff0_13*dmats1[13][10] + coeff0_14*dmats1[14][10];
502
 
          new_coeff0_11 = coeff0_0*dmats1[0][11] + coeff0_1*dmats1[1][11] + coeff0_2*dmats1[2][11] + coeff0_3*dmats1[3][11] + coeff0_4*dmats1[4][11] + coeff0_5*dmats1[5][11] + coeff0_6*dmats1[6][11] + coeff0_7*dmats1[7][11] + coeff0_8*dmats1[8][11] + coeff0_9*dmats1[9][11] + coeff0_10*dmats1[10][11] + coeff0_11*dmats1[11][11] + coeff0_12*dmats1[12][11] + coeff0_13*dmats1[13][11] + coeff0_14*dmats1[14][11];
503
 
          new_coeff0_12 = coeff0_0*dmats1[0][12] + coeff0_1*dmats1[1][12] + coeff0_2*dmats1[2][12] + coeff0_3*dmats1[3][12] + coeff0_4*dmats1[4][12] + coeff0_5*dmats1[5][12] + coeff0_6*dmats1[6][12] + coeff0_7*dmats1[7][12] + coeff0_8*dmats1[8][12] + coeff0_9*dmats1[9][12] + coeff0_10*dmats1[10][12] + coeff0_11*dmats1[11][12] + coeff0_12*dmats1[12][12] + coeff0_13*dmats1[13][12] + coeff0_14*dmats1[14][12];
504
 
          new_coeff0_13 = coeff0_0*dmats1[0][13] + coeff0_1*dmats1[1][13] + coeff0_2*dmats1[2][13] + coeff0_3*dmats1[3][13] + coeff0_4*dmats1[4][13] + coeff0_5*dmats1[5][13] + coeff0_6*dmats1[6][13] + coeff0_7*dmats1[7][13] + coeff0_8*dmats1[8][13] + coeff0_9*dmats1[9][13] + coeff0_10*dmats1[10][13] + coeff0_11*dmats1[11][13] + coeff0_12*dmats1[12][13] + coeff0_13*dmats1[13][13] + coeff0_14*dmats1[14][13];
505
 
          new_coeff0_14 = coeff0_0*dmats1[0][14] + coeff0_1*dmats1[1][14] + coeff0_2*dmats1[2][14] + coeff0_3*dmats1[3][14] + coeff0_4*dmats1[4][14] + coeff0_5*dmats1[5][14] + coeff0_6*dmats1[6][14] + coeff0_7*dmats1[7][14] + coeff0_8*dmats1[8][14] + coeff0_9*dmats1[9][14] + coeff0_10*dmats1[10][14] + coeff0_11*dmats1[11][14] + coeff0_12*dmats1[12][14] + coeff0_13*dmats1[13][14] + coeff0_14*dmats1[14][14];
506
 
        }
507
 
    
508
 
      }
509
 
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
510
 
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5 + new_coeff0_6*basisvalue6 + new_coeff0_7*basisvalue7 + new_coeff0_8*basisvalue8 + new_coeff0_9*basisvalue9 + new_coeff0_10*basisvalue10 + new_coeff0_11*basisvalue11 + new_coeff0_12*basisvalue12 + new_coeff0_13*basisvalue13 + new_coeff0_14*basisvalue14;
511
 
    }
512
 
    
513
 
    // Transform derivatives back to physical element
514
 
    for (unsigned int row = 0; row < num_derivatives; row++)
515
 
    {
516
 
      for (unsigned int col = 0; col < num_derivatives; col++)
517
 
      {
518
 
        values[row] += transform[row][col]*derivatives[col];
519
 
      }
520
 
    }
521
 
    // Delete pointer to array of derivatives on FIAT element
522
 
    delete [] derivatives;
523
 
    
524
 
    // Delete pointer to array of combinations of derivatives and transform
525
 
    for (unsigned int row = 0; row < num_derivatives; row++)
526
 
    {
527
 
      delete [] combinations[row];
528
 
      delete [] transform[row];
529
 
    }
530
 
    
531
 
    delete [] combinations;
532
 
    delete [] transform;
533
 
  }
534
 
 
535
 
  /// Evaluate order n derivatives of all basis functions at given point in cell
536
 
  virtual void evaluate_basis_derivatives_all(unsigned int n,
537
 
                                              double* values,
538
 
                                              const double* coordinates,
539
 
                                              const ufc::cell& c) const
540
 
  {
541
 
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
542
 
  }
543
 
 
544
 
  /// Evaluate linear functional for dof i on the function f
545
 
  virtual double evaluate_dof(unsigned int i,
546
 
                              const ufc::function& f,
547
 
                              const ufc::cell& c) const
548
 
  {
549
 
    // The reference points, direction and weights:
550
 
    static const double X[12][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.75, 0.25}}, {{0.5, 0.5}}, {{0.25, 0.75}}, {{0, 0.25}}, {{0, 0.5}}, {{0, 0.75}}, {{0.25, 0}}, {{0.5, 0}}, {{0.75, 0}}};
551
 
    static const double W[12][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
552
 
    static const double D[12][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
553
 
    
554
 
    const double * const * x = c.coordinates;
555
 
    double result = 0.0;
556
 
    // Iterate over the points:
557
 
    // Evaluate basis functions for affine mapping
558
 
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
559
 
    const double w1 = X[i][0][0];
560
 
    const double w2 = X[i][0][1];
561
 
    
562
 
    // Compute affine mapping y = F(X)
563
 
    double y[2];
564
 
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
565
 
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
566
 
    
567
 
    // Evaluate function at physical points
568
 
    double values[1];
569
 
    f.evaluate(values, y, c);
570
 
    
571
 
    // Map function values using appropriate mapping
572
 
    // Affine map: Do nothing
573
 
    
574
 
    // Note that we do not map the weights (yet).
575
 
    
576
 
    // Take directional components
577
 
    for(int k = 0; k < 1; k++)
578
 
      result += values[k]*D[i][0][k];
579
 
    // Multiply by weights
580
 
    result *= W[i][0];
581
 
    
582
 
    return result;
583
 
  }
584
 
 
585
 
  /// Evaluate linear functionals for all dofs on the function f
586
 
  virtual void evaluate_dofs(double* values,
587
 
                             const ufc::function& f,
588
 
                             const ufc::cell& c) const
589
 
  {
590
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
591
 
  }
592
 
 
593
 
  /// Interpolate vertex values from dof values
594
 
  virtual void interpolate_vertex_values(double* vertex_values,
595
 
                                         const double* dof_values,
596
 
                                         const ufc::cell& c) const
597
 
  {
598
 
    // Evaluate at vertices and use affine mapping
599
 
    vertex_values[0] = dof_values[0];
600
 
    vertex_values[1] = dof_values[1];
601
 
    vertex_values[2] = dof_values[2];
602
 
  }
603
 
 
604
 
  /// Return the number of sub elements (for a mixed element)
605
 
  virtual unsigned int num_sub_elements() const
606
 
  {
607
 
    return 1;
608
 
  }
609
 
 
610
 
  /// Create a new finite element for sub element i (for a mixed element)
611
 
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
612
 
  {
613
 
    return new elementrestriction_0_finite_element_0_0();
614
 
  }
615
 
 
616
 
};
617
 
 
618
 
/// This class defines the interface for a finite element.
619
 
 
620
 
class elementrestriction_0_finite_element_0_1: public ufc::finite_element
621
 
{
622
 
public:
623
 
 
624
 
  /// Constructor
625
 
  elementrestriction_0_finite_element_0_1() : ufc::finite_element()
626
 
  {
627
 
    // Do nothing
628
 
  }
629
 
 
630
 
  /// Destructor
631
 
  virtual ~elementrestriction_0_finite_element_0_1()
632
 
  {
633
 
    // Do nothing
634
 
  }
635
 
 
636
 
  /// Return a string identifying the finite element
637
 
  virtual const char* signature() const
638
 
  {
639
 
    return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 4)";
640
 
  }
641
 
 
642
 
  /// Return the cell shape
643
 
  virtual ufc::shape cell_shape() const
644
 
  {
645
 
    return ufc::triangle;
646
 
  }
647
 
 
648
 
  /// Return the dimension of the finite element function space
649
 
  virtual unsigned int space_dimension() const
650
 
  {
651
 
    return 12;
652
 
  }
653
 
 
654
 
  /// Return the rank of the value space
655
 
  virtual unsigned int value_rank() const
656
 
  {
657
 
    return 0;
658
 
  }
659
 
 
660
 
  /// Return the dimension of the value space for axis i
661
 
  virtual unsigned int value_dimension(unsigned int i) const
662
 
  {
663
 
    return 1;
664
 
  }
665
 
 
666
 
  /// Evaluate basis function i at given point in cell
667
 
  virtual void evaluate_basis(unsigned int i,
668
 
                              double* values,
669
 
                              const double* coordinates,
670
 
                              const ufc::cell& c) const
671
 
  {
672
 
    // Extract vertex coordinates
673
 
    const double * const * element_coordinates = c.coordinates;
674
 
    
675
 
    // Compute Jacobian of affine map from reference cell
676
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
677
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
678
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
679
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
680
 
    
681
 
    // Compute determinant of Jacobian
682
 
    const double detJ = J_00*J_11 - J_01*J_10;
683
 
    
684
 
    // Compute inverse of Jacobian
685
 
    
686
 
    // Get coordinates and map to the reference (UFC) element
687
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
688
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
689
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
690
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
691
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
692
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
693
 
    
694
 
    // Map coordinates to the reference square
695
 
    if (std::abs(y - 1.0) < 1e-08)
696
 
      x = -1.0;
697
 
    else
698
 
      x = 2.0 *x/(1.0 - y) - 1.0;
699
 
    y = 2.0*y - 1.0;
700
 
    
701
 
    // Reset values
702
 
    *values = 0;
703
 
    
704
 
    // Map degree of freedom to element degree of freedom
705
 
    const unsigned int dof = i;
706
 
    
707
 
    // Generate scalings
708
 
    const double scalings_y_0 = 1;
709
 
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
710
 
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
711
 
    const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
712
 
    const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
713
 
    
714
 
    // Compute psitilde_a
715
 
    const double psitilde_a_0 = 1;
716
 
    const double psitilde_a_1 = x;
717
 
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
718
 
    const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;
719
 
    const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
720
 
    
721
 
    // Compute psitilde_bs
722
 
    const double psitilde_bs_0_0 = 1;
723
 
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
724
 
    const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
725
 
    const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
726
 
    const double psitilde_bs_0_4 = 0.0285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571*psitilde_bs_0_2;
727
 
    const double psitilde_bs_1_0 = 1;
728
 
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
729
 
    const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
730
 
    const double psitilde_bs_1_3 = 0.285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714*psitilde_bs_1_1;
731
 
    const double psitilde_bs_2_0 = 1;
732
 
    const double psitilde_bs_2_1 = 3.5*y + 2.5;
733
 
    const double psitilde_bs_2_2 = 1.02040816*psitilde_bs_2_1 + 2.57142857*y*psitilde_bs_2_1 - 0.551020408*psitilde_bs_2_0;
734
 
    const double psitilde_bs_3_0 = 1;
735
 
    const double psitilde_bs_3_1 = 4.5*y + 3.5;
736
 
    const double psitilde_bs_4_0 = 1;
737
 
    
738
 
    // Compute basisvalues
739
 
    const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
740
 
    const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
741
 
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
742
 
    const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
743
 
    const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
744
 
    const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
745
 
    const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
746
 
    const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
747
 
    const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
748
 
    const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
749
 
    const double basisvalue10 = 4.74341649*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
750
 
    const double basisvalue11 = 4.18330013*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
751
 
    const double basisvalue12 = 3.53553391*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
752
 
    const double basisvalue13 = 2.73861279*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
753
 
    const double basisvalue14 = 1.58113883*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
754
 
    
755
 
    // Table(s) of coefficients
756
 
    static const double coefficients0[12][15] = \
757
 
    {{0, -0.0412393049, -0.0238095238, 0.0289800295, 0.0224478343, 0.0129602632, -0.0395942581, -0.0334632557, -0.0259205264, -0.0149652229, 0.0321247254, 0.0283313448, 0.0239443566, 0.0185472189, 0.0107082418},
758
 
    {0, 0.0412393049, -0.0238095238, 0.0289800295, -0.0224478343, 0.0129602632, 0.0395942581, -0.0334632557, 0.0259205264, -0.0149652229, 0.0321247254, -0.0283313448, 0.0239443566, -0.0185472189, 0.0107082418},
759
 
    {0, 0, 0.0476190476, 0, 0, 0.0388807896, 0, 0, 0, 0.0598608915, 0, 0, 0, 0, 0.0535412091},
760
 
    {0.125707872, 0.131965776, -0.0253968254, 0.139104142, -0.0718330698, 0.0311046317, 0.0633508129, 0.0267706045, -0.0622092633, 0.0478887132, 0, 0.0566626896, -0.0838052481, 0.083462485, -0.0535412091},
761
 
    {-0.0314269681, 0.010997148, 0.00634920635, 0, 0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, -0.139104142, 0.107082418},
762
 
    {0.125707872, 0.0439885919, 0.126984127, 0, 0.0359165349, 0.155523158, 0, 0, 0.103682106, -0.0119721783, 0, 0, 0, 0.0927360944, -0.107082418},
763
 
    {0.125707872, -0.131965776, -0.0253968254, 0.139104142, 0.0718330698, 0.0311046317, -0.0633508129, 0.0267706045, 0.0622092633, 0.0478887132, 0, -0.0566626896, -0.0838052481, -0.083462485, -0.0535412091},
764
 
    {-0.0314269681, -0.010997148, 0.00634920635, 0, -0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, 0.139104142, 0.107082418},
765
 
    {0.125707872, -0.0439885919, 0.126984127, 0, -0.0359165349, 0.155523158, 0, 0, -0.103682106, -0.0119721783, 0, 0, 0, -0.0927360944, -0.107082418},
766
 
    {0.125707872, -0.0879771839, -0.101587302, 0.0927360944, 0.107749605, 0.0725774739, 0.0791885161, -0.0133853023, -0.0518410528, -0.0419026241, -0.128498902, -0.0566626896, -0.0119721783, 0.00927360944, 0.0107082418},
767
 
    {-0.0314269681, 0, -0.0126984127, -0.243432248, 0, 0.0544331054, 0, 0.0936971159, 0, -0.0419026241, 0.192748353, 0, -0.0239443566, 0, 0.0107082418},
768
 
    {0.125707872, 0.0879771839, -0.101587302, 0.0927360944, -0.107749605, 0.0725774739, -0.0791885161, -0.0133853023, 0.0518410528, -0.0419026241, -0.128498902, 0.0566626896, -0.0119721783, -0.00927360944, 0.0107082418}};
769
 
    
770
 
    // Extract relevant coefficients
771
 
    const double coeff0_0 = coefficients0[dof][0];
772
 
    const double coeff0_1 = coefficients0[dof][1];
773
 
    const double coeff0_2 = coefficients0[dof][2];
774
 
    const double coeff0_3 = coefficients0[dof][3];
775
 
    const double coeff0_4 = coefficients0[dof][4];
776
 
    const double coeff0_5 = coefficients0[dof][5];
777
 
    const double coeff0_6 = coefficients0[dof][6];
778
 
    const double coeff0_7 = coefficients0[dof][7];
779
 
    const double coeff0_8 = coefficients0[dof][8];
780
 
    const double coeff0_9 = coefficients0[dof][9];
781
 
    const double coeff0_10 = coefficients0[dof][10];
782
 
    const double coeff0_11 = coefficients0[dof][11];
783
 
    const double coeff0_12 = coefficients0[dof][12];
784
 
    const double coeff0_13 = coefficients0[dof][13];
785
 
    const double coeff0_14 = coefficients0[dof][14];
786
 
    
787
 
    // Compute value(s)
788
 
    *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5 + coeff0_6*basisvalue6 + coeff0_7*basisvalue7 + coeff0_8*basisvalue8 + coeff0_9*basisvalue9 + coeff0_10*basisvalue10 + coeff0_11*basisvalue11 + coeff0_12*basisvalue12 + coeff0_13*basisvalue13 + coeff0_14*basisvalue14;
789
 
  }
790
 
 
791
 
  /// Evaluate all basis functions at given point in cell
792
 
  virtual void evaluate_basis_all(double* values,
793
 
                                  const double* coordinates,
794
 
                                  const ufc::cell& c) const
795
 
  {
796
 
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
797
 
  }
798
 
 
799
 
  /// Evaluate order n derivatives of basis function i at given point in cell
800
 
  virtual void evaluate_basis_derivatives(unsigned int i,
801
 
                                          unsigned int n,
802
 
                                          double* values,
803
 
                                          const double* coordinates,
804
 
                                          const ufc::cell& c) const
805
 
  {
806
 
    // Extract vertex coordinates
807
 
    const double * const * element_coordinates = c.coordinates;
808
 
    
809
 
    // Compute Jacobian of affine map from reference cell
810
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
811
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
812
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
813
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
814
 
    
815
 
    // Compute determinant of Jacobian
816
 
    const double detJ = J_00*J_11 - J_01*J_10;
817
 
    
818
 
    // Compute inverse of Jacobian
819
 
    
820
 
    // Get coordinates and map to the reference (UFC) element
821
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
822
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
823
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
824
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
825
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
826
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
827
 
    
828
 
    // Map coordinates to the reference square
829
 
    if (std::abs(y - 1.0) < 1e-08)
830
 
      x = -1.0;
831
 
    else
832
 
      x = 2.0 *x/(1.0 - y) - 1.0;
833
 
    y = 2.0*y - 1.0;
834
 
    
835
 
    // Compute number of derivatives
836
 
    unsigned int num_derivatives = 1;
837
 
    
838
 
    for (unsigned int j = 0; j < n; j++)
839
 
      num_derivatives *= 2;
840
 
    
841
 
    
842
 
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
843
 
    unsigned int **combinations = new unsigned int *[num_derivatives];
844
 
    
845
 
    for (unsigned int j = 0; j < num_derivatives; j++)
846
 
    {
847
 
      combinations[j] = new unsigned int [n];
848
 
      for (unsigned int k = 0; k < n; k++)
849
 
        combinations[j][k] = 0;
850
 
    }
851
 
    
852
 
    // Generate combinations of derivatives
853
 
    for (unsigned int row = 1; row < num_derivatives; row++)
854
 
    {
855
 
      for (unsigned int num = 0; num < row; num++)
856
 
      {
857
 
        for (unsigned int col = n-1; col+1 > 0; col--)
858
 
        {
859
 
          if (combinations[row][col] + 1 > 1)
860
 
            combinations[row][col] = 0;
861
 
          else
862
 
          {
863
 
            combinations[row][col] += 1;
864
 
            break;
865
 
          }
866
 
        }
867
 
      }
868
 
    }
869
 
    
870
 
    // Compute inverse of Jacobian
871
 
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
872
 
    
873
 
    // Declare transformation matrix
874
 
    // Declare pointer to two dimensional array and initialise
875
 
    double **transform = new double *[num_derivatives];
876
 
    
877
 
    for (unsigned int j = 0; j < num_derivatives; j++)
878
 
    {
879
 
      transform[j] = new double [num_derivatives];
880
 
      for (unsigned int k = 0; k < num_derivatives; k++)
881
 
        transform[j][k] = 1;
882
 
    }
883
 
    
884
 
    // Construct transformation matrix
885
 
    for (unsigned int row = 0; row < num_derivatives; row++)
886
 
    {
887
 
      for (unsigned int col = 0; col < num_derivatives; col++)
888
 
      {
889
 
        for (unsigned int k = 0; k < n; k++)
890
 
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
891
 
      }
892
 
    }
893
 
    
894
 
    // Reset values
895
 
    for (unsigned int j = 0; j < 1*num_derivatives; j++)
896
 
      values[j] = 0;
897
 
    
898
 
    // Map degree of freedom to element degree of freedom
899
 
    const unsigned int dof = i;
900
 
    
901
 
    // Generate scalings
902
 
    const double scalings_y_0 = 1;
903
 
    const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
904
 
    const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
905
 
    const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
906
 
    const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
907
 
    
908
 
    // Compute psitilde_a
909
 
    const double psitilde_a_0 = 1;
910
 
    const double psitilde_a_1 = x;
911
 
    const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
912
 
    const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;
913
 
    const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
914
 
    
915
 
    // Compute psitilde_bs
916
 
    const double psitilde_bs_0_0 = 1;
917
 
    const double psitilde_bs_0_1 = 1.5*y + 0.5;
918
 
    const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
919
 
    const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
920
 
    const double psitilde_bs_0_4 = 0.0285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571*psitilde_bs_0_2;
921
 
    const double psitilde_bs_1_0 = 1;
922
 
    const double psitilde_bs_1_1 = 2.5*y + 1.5;
923
 
    const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
924
 
    const double psitilde_bs_1_3 = 0.285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714*psitilde_bs_1_1;
925
 
    const double psitilde_bs_2_0 = 1;
926
 
    const double psitilde_bs_2_1 = 3.5*y + 2.5;
927
 
    const double psitilde_bs_2_2 = 1.02040816*psitilde_bs_2_1 + 2.57142857*y*psitilde_bs_2_1 - 0.551020408*psitilde_bs_2_0;
928
 
    const double psitilde_bs_3_0 = 1;
929
 
    const double psitilde_bs_3_1 = 4.5*y + 3.5;
930
 
    const double psitilde_bs_4_0 = 1;
931
 
    
932
 
    // Compute basisvalues
933
 
    const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
934
 
    const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
935
 
    const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
936
 
    const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
937
 
    const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
938
 
    const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
939
 
    const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
940
 
    const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
941
 
    const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
942
 
    const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
943
 
    const double basisvalue10 = 4.74341649*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
944
 
    const double basisvalue11 = 4.18330013*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
945
 
    const double basisvalue12 = 3.53553391*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
946
 
    const double basisvalue13 = 2.73861279*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
947
 
    const double basisvalue14 = 1.58113883*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
948
 
    
949
 
    // Table(s) of coefficients
950
 
    static const double coefficients0[12][15] = \
951
 
    {{0, -0.0412393049, -0.0238095238, 0.0289800295, 0.0224478343, 0.0129602632, -0.0395942581, -0.0334632557, -0.0259205264, -0.0149652229, 0.0321247254, 0.0283313448, 0.0239443566, 0.0185472189, 0.0107082418},
952
 
    {0, 0.0412393049, -0.0238095238, 0.0289800295, -0.0224478343, 0.0129602632, 0.0395942581, -0.0334632557, 0.0259205264, -0.0149652229, 0.0321247254, -0.0283313448, 0.0239443566, -0.0185472189, 0.0107082418},
953
 
    {0, 0, 0.0476190476, 0, 0, 0.0388807896, 0, 0, 0, 0.0598608915, 0, 0, 0, 0, 0.0535412091},
954
 
    {0.125707872, 0.131965776, -0.0253968254, 0.139104142, -0.0718330698, 0.0311046317, 0.0633508129, 0.0267706045, -0.0622092633, 0.0478887132, 0, 0.0566626896, -0.0838052481, 0.083462485, -0.0535412091},
955
 
    {-0.0314269681, 0.010997148, 0.00634920635, 0, 0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, -0.139104142, 0.107082418},
956
 
    {0.125707872, 0.0439885919, 0.126984127, 0, 0.0359165349, 0.155523158, 0, 0, 0.103682106, -0.0119721783, 0, 0, 0, 0.0927360944, -0.107082418},
957
 
    {0.125707872, -0.131965776, -0.0253968254, 0.139104142, 0.0718330698, 0.0311046317, -0.0633508129, 0.0267706045, 0.0622092633, 0.0478887132, 0, -0.0566626896, -0.0838052481, -0.083462485, -0.0535412091},
958
 
    {-0.0314269681, -0.010997148, 0.00634920635, 0, -0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, 0.139104142, 0.107082418},
959
 
    {0.125707872, -0.0439885919, 0.126984127, 0, -0.0359165349, 0.155523158, 0, 0, -0.103682106, -0.0119721783, 0, 0, 0, -0.0927360944, -0.107082418},
960
 
    {0.125707872, -0.0879771839, -0.101587302, 0.0927360944, 0.107749605, 0.0725774739, 0.0791885161, -0.0133853023, -0.0518410528, -0.0419026241, -0.128498902, -0.0566626896, -0.0119721783, 0.00927360944, 0.0107082418},
961
 
    {-0.0314269681, 0, -0.0126984127, -0.243432248, 0, 0.0544331054, 0, 0.0936971159, 0, -0.0419026241, 0.192748353, 0, -0.0239443566, 0, 0.0107082418},
962
 
    {0.125707872, 0.0879771839, -0.101587302, 0.0927360944, -0.107749605, 0.0725774739, -0.0791885161, -0.0133853023, 0.0518410528, -0.0419026241, -0.128498902, 0.0566626896, -0.0119721783, -0.00927360944, 0.0107082418}};
963
 
    
964
 
    // Interesting (new) part
965
 
    // Tables of derivatives of the polynomial base (transpose)
966
 
    static const double dmats0[15][15] = \
967
 
    {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
968
 
    {4.89897949, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
969
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
970
 
    {0, 9.48683298, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
971
 
    {4, 0, 7.07106781, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
972
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
973
 
    {5.29150262, 0, -2.99332591, 13.662601, 0, 0.611010093, 0, 0, 0, 0, 0, 0, 0, 0, 0},
974
 
    {0, 4.38178046, 0, 0, 12.5219807, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
975
 
    {3.46410162, 0, 7.83836718, 0, 0, 8.4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
976
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
977
 
    {0, 10.9544512, 0, 0, -3.83325939, 0, 17.7482393, 0, 0.553283335, 0, 0, 0, 0, 0, 0},
978
 
    {4.73286383, 0, 3.34664011, 4.3643578, 0, -5.07468038, 0, 17.0084013, 0, 1.52127766, 0, 0, 0, 0, 0},
979
 
    {0, 2.44948974, 0, 0, 9.14285714, 0, 0, 0, 14.8461498, 0, 0, 0, 0, 0, 0},
980
 
    {3.09838668, 0, 7.66811581, 0, 0, 10.7331263, 0, 0, 0, 9.29516003, 0, 0, 0, 0, 0},
981
 
    {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
982
 
    
983
 
    static const double dmats1[15][15] = \
984
 
    {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
985
 
    {2.44948974, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
986
 
    {4.24264069, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
987
 
    {2.5819889, 4.74341649, -0.912870929, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
988
 
    {2, 6.12372436, 3.53553391, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
989
 
    {-2.30940108, 0, 8.16496581, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
990
 
    {2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.305505046, 0, 0, 0, 0, 0, 0, 0, 0, 0},
991
 
    {2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0, 0, 0, 0, 0, 0, 0, 0, 0},
992
 
    {1.73205081, -5.09116882, 3.91918359, 0, 9.69948452, 4.2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
993
 
    {5, 0, -2.82842712, 0, 0, 12.1243557, 0, 0, 0, 0, 0, 0, 0, 0, 0},
994
 
    {2.68328157, 5.47722558, -1.8973666, 7.42307489, -1.91662969, 0.663940002, 8.87411967, -1.07142857, 0.276641668, -0.0958314847, 0, 0, 0, 0, 0},
995
 
    {2.36643191, 2.89827535, 1.67332005, 2.1821789, 5.74704893, -2.53734019, 10.0623059, 8.50420064, -2.19577516, 0.760638829, 0, 0, 0, 0, 0},
996
 
    {2, 1.22474487, 3.53553391, -7.37711114, 4.57142857, 1.6495722, 0, 11.4997782, 7.42307489, -2.57142857, 0, 0, 0, 0, 0},
997
 
    {1.54919334, 6.64078309, 3.8340579, 0, -6.19677335, 5.36656315, 0, 0, 13.4164079, 4.64758002, 0, 0, 0, 0, 0},
998
 
    {-3.57770876, 0, 8.85437745, 0, 0, -3.09838668, 0, 0, 0, 16.0996894, 0, 0, 0, 0, 0}};
999
 
    
1000
 
    // Compute reference derivatives
1001
 
    // Declare pointer to array of derivatives on FIAT element
1002
 
    double *derivatives = new double [num_derivatives];
1003
 
    
1004
 
    // Declare coefficients
1005
 
    double coeff0_0 = 0;
1006
 
    double coeff0_1 = 0;
1007
 
    double coeff0_2 = 0;
1008
 
    double coeff0_3 = 0;
1009
 
    double coeff0_4 = 0;
1010
 
    double coeff0_5 = 0;
1011
 
    double coeff0_6 = 0;
1012
 
    double coeff0_7 = 0;
1013
 
    double coeff0_8 = 0;
1014
 
    double coeff0_9 = 0;
1015
 
    double coeff0_10 = 0;
1016
 
    double coeff0_11 = 0;
1017
 
    double coeff0_12 = 0;
1018
 
    double coeff0_13 = 0;
1019
 
    double coeff0_14 = 0;
1020
 
    
1021
 
    // Declare new coefficients
1022
 
    double new_coeff0_0 = 0;
1023
 
    double new_coeff0_1 = 0;
1024
 
    double new_coeff0_2 = 0;
1025
 
    double new_coeff0_3 = 0;
1026
 
    double new_coeff0_4 = 0;
1027
 
    double new_coeff0_5 = 0;
1028
 
    double new_coeff0_6 = 0;
1029
 
    double new_coeff0_7 = 0;
1030
 
    double new_coeff0_8 = 0;
1031
 
    double new_coeff0_9 = 0;
1032
 
    double new_coeff0_10 = 0;
1033
 
    double new_coeff0_11 = 0;
1034
 
    double new_coeff0_12 = 0;
1035
 
    double new_coeff0_13 = 0;
1036
 
    double new_coeff0_14 = 0;
1037
 
    
1038
 
    // Loop possible derivatives
1039
 
    for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1040
 
    {
1041
 
      // Get values from coefficients array
1042
 
      new_coeff0_0 = coefficients0[dof][0];
1043
 
      new_coeff0_1 = coefficients0[dof][1];
1044
 
      new_coeff0_2 = coefficients0[dof][2];
1045
 
      new_coeff0_3 = coefficients0[dof][3];
1046
 
      new_coeff0_4 = coefficients0[dof][4];
1047
 
      new_coeff0_5 = coefficients0[dof][5];
1048
 
      new_coeff0_6 = coefficients0[dof][6];
1049
 
      new_coeff0_7 = coefficients0[dof][7];
1050
 
      new_coeff0_8 = coefficients0[dof][8];
1051
 
      new_coeff0_9 = coefficients0[dof][9];
1052
 
      new_coeff0_10 = coefficients0[dof][10];
1053
 
      new_coeff0_11 = coefficients0[dof][11];
1054
 
      new_coeff0_12 = coefficients0[dof][12];
1055
 
      new_coeff0_13 = coefficients0[dof][13];
1056
 
      new_coeff0_14 = coefficients0[dof][14];
1057
 
    
1058
 
      // Loop derivative order
1059
 
      for (unsigned int j = 0; j < n; j++)
1060
 
      {
1061
 
        // Update old coefficients
1062
 
        coeff0_0 = new_coeff0_0;
1063
 
        coeff0_1 = new_coeff0_1;
1064
 
        coeff0_2 = new_coeff0_2;
1065
 
        coeff0_3 = new_coeff0_3;
1066
 
        coeff0_4 = new_coeff0_4;
1067
 
        coeff0_5 = new_coeff0_5;
1068
 
        coeff0_6 = new_coeff0_6;
1069
 
        coeff0_7 = new_coeff0_7;
1070
 
        coeff0_8 = new_coeff0_8;
1071
 
        coeff0_9 = new_coeff0_9;
1072
 
        coeff0_10 = new_coeff0_10;
1073
 
        coeff0_11 = new_coeff0_11;
1074
 
        coeff0_12 = new_coeff0_12;
1075
 
        coeff0_13 = new_coeff0_13;
1076
 
        coeff0_14 = new_coeff0_14;
1077
 
    
1078
 
        if(combinations[deriv_num][j] == 0)
1079
 
        {
1080
 
          new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0] + coeff0_6*dmats0[6][0] + coeff0_7*dmats0[7][0] + coeff0_8*dmats0[8][0] + coeff0_9*dmats0[9][0] + coeff0_10*dmats0[10][0] + coeff0_11*dmats0[11][0] + coeff0_12*dmats0[12][0] + coeff0_13*dmats0[13][0] + coeff0_14*dmats0[14][0];
1081
 
          new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1] + coeff0_6*dmats0[6][1] + coeff0_7*dmats0[7][1] + coeff0_8*dmats0[8][1] + coeff0_9*dmats0[9][1] + coeff0_10*dmats0[10][1] + coeff0_11*dmats0[11][1] + coeff0_12*dmats0[12][1] + coeff0_13*dmats0[13][1] + coeff0_14*dmats0[14][1];
1082
 
          new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2] + coeff0_6*dmats0[6][2] + coeff0_7*dmats0[7][2] + coeff0_8*dmats0[8][2] + coeff0_9*dmats0[9][2] + coeff0_10*dmats0[10][2] + coeff0_11*dmats0[11][2] + coeff0_12*dmats0[12][2] + coeff0_13*dmats0[13][2] + coeff0_14*dmats0[14][2];
1083
 
          new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3] + coeff0_6*dmats0[6][3] + coeff0_7*dmats0[7][3] + coeff0_8*dmats0[8][3] + coeff0_9*dmats0[9][3] + coeff0_10*dmats0[10][3] + coeff0_11*dmats0[11][3] + coeff0_12*dmats0[12][3] + coeff0_13*dmats0[13][3] + coeff0_14*dmats0[14][3];
1084
 
          new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4] + coeff0_6*dmats0[6][4] + coeff0_7*dmats0[7][4] + coeff0_8*dmats0[8][4] + coeff0_9*dmats0[9][4] + coeff0_10*dmats0[10][4] + coeff0_11*dmats0[11][4] + coeff0_12*dmats0[12][4] + coeff0_13*dmats0[13][4] + coeff0_14*dmats0[14][4];
1085
 
          new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5] + coeff0_6*dmats0[6][5] + coeff0_7*dmats0[7][5] + coeff0_8*dmats0[8][5] + coeff0_9*dmats0[9][5] + coeff0_10*dmats0[10][5] + coeff0_11*dmats0[11][5] + coeff0_12*dmats0[12][5] + coeff0_13*dmats0[13][5] + coeff0_14*dmats0[14][5];
1086
 
          new_coeff0_6 = coeff0_0*dmats0[0][6] + coeff0_1*dmats0[1][6] + coeff0_2*dmats0[2][6] + coeff0_3*dmats0[3][6] + coeff0_4*dmats0[4][6] + coeff0_5*dmats0[5][6] + coeff0_6*dmats0[6][6] + coeff0_7*dmats0[7][6] + coeff0_8*dmats0[8][6] + coeff0_9*dmats0[9][6] + coeff0_10*dmats0[10][6] + coeff0_11*dmats0[11][6] + coeff0_12*dmats0[12][6] + coeff0_13*dmats0[13][6] + coeff0_14*dmats0[14][6];
1087
 
          new_coeff0_7 = coeff0_0*dmats0[0][7] + coeff0_1*dmats0[1][7] + coeff0_2*dmats0[2][7] + coeff0_3*dmats0[3][7] + coeff0_4*dmats0[4][7] + coeff0_5*dmats0[5][7] + coeff0_6*dmats0[6][7] + coeff0_7*dmats0[7][7] + coeff0_8*dmats0[8][7] + coeff0_9*dmats0[9][7] + coeff0_10*dmats0[10][7] + coeff0_11*dmats0[11][7] + coeff0_12*dmats0[12][7] + coeff0_13*dmats0[13][7] + coeff0_14*dmats0[14][7];
1088
 
          new_coeff0_8 = coeff0_0*dmats0[0][8] + coeff0_1*dmats0[1][8] + coeff0_2*dmats0[2][8] + coeff0_3*dmats0[3][8] + coeff0_4*dmats0[4][8] + coeff0_5*dmats0[5][8] + coeff0_6*dmats0[6][8] + coeff0_7*dmats0[7][8] + coeff0_8*dmats0[8][8] + coeff0_9*dmats0[9][8] + coeff0_10*dmats0[10][8] + coeff0_11*dmats0[11][8] + coeff0_12*dmats0[12][8] + coeff0_13*dmats0[13][8] + coeff0_14*dmats0[14][8];
1089
 
          new_coeff0_9 = coeff0_0*dmats0[0][9] + coeff0_1*dmats0[1][9] + coeff0_2*dmats0[2][9] + coeff0_3*dmats0[3][9] + coeff0_4*dmats0[4][9] + coeff0_5*dmats0[5][9] + coeff0_6*dmats0[6][9] + coeff0_7*dmats0[7][9] + coeff0_8*dmats0[8][9] + coeff0_9*dmats0[9][9] + coeff0_10*dmats0[10][9] + coeff0_11*dmats0[11][9] + coeff0_12*dmats0[12][9] + coeff0_13*dmats0[13][9] + coeff0_14*dmats0[14][9];
1090
 
          new_coeff0_10 = coeff0_0*dmats0[0][10] + coeff0_1*dmats0[1][10] + coeff0_2*dmats0[2][10] + coeff0_3*dmats0[3][10] + coeff0_4*dmats0[4][10] + coeff0_5*dmats0[5][10] + coeff0_6*dmats0[6][10] + coeff0_7*dmats0[7][10] + coeff0_8*dmats0[8][10] + coeff0_9*dmats0[9][10] + coeff0_10*dmats0[10][10] + coeff0_11*dmats0[11][10] + coeff0_12*dmats0[12][10] + coeff0_13*dmats0[13][10] + coeff0_14*dmats0[14][10];
1091
 
          new_coeff0_11 = coeff0_0*dmats0[0][11] + coeff0_1*dmats0[1][11] + coeff0_2*dmats0[2][11] + coeff0_3*dmats0[3][11] + coeff0_4*dmats0[4][11] + coeff0_5*dmats0[5][11] + coeff0_6*dmats0[6][11] + coeff0_7*dmats0[7][11] + coeff0_8*dmats0[8][11] + coeff0_9*dmats0[9][11] + coeff0_10*dmats0[10][11] + coeff0_11*dmats0[11][11] + coeff0_12*dmats0[12][11] + coeff0_13*dmats0[13][11] + coeff0_14*dmats0[14][11];
1092
 
          new_coeff0_12 = coeff0_0*dmats0[0][12] + coeff0_1*dmats0[1][12] + coeff0_2*dmats0[2][12] + coeff0_3*dmats0[3][12] + coeff0_4*dmats0[4][12] + coeff0_5*dmats0[5][12] + coeff0_6*dmats0[6][12] + coeff0_7*dmats0[7][12] + coeff0_8*dmats0[8][12] + coeff0_9*dmats0[9][12] + coeff0_10*dmats0[10][12] + coeff0_11*dmats0[11][12] + coeff0_12*dmats0[12][12] + coeff0_13*dmats0[13][12] + coeff0_14*dmats0[14][12];
1093
 
          new_coeff0_13 = coeff0_0*dmats0[0][13] + coeff0_1*dmats0[1][13] + coeff0_2*dmats0[2][13] + coeff0_3*dmats0[3][13] + coeff0_4*dmats0[4][13] + coeff0_5*dmats0[5][13] + coeff0_6*dmats0[6][13] + coeff0_7*dmats0[7][13] + coeff0_8*dmats0[8][13] + coeff0_9*dmats0[9][13] + coeff0_10*dmats0[10][13] + coeff0_11*dmats0[11][13] + coeff0_12*dmats0[12][13] + coeff0_13*dmats0[13][13] + coeff0_14*dmats0[14][13];
1094
 
          new_coeff0_14 = coeff0_0*dmats0[0][14] + coeff0_1*dmats0[1][14] + coeff0_2*dmats0[2][14] + coeff0_3*dmats0[3][14] + coeff0_4*dmats0[4][14] + coeff0_5*dmats0[5][14] + coeff0_6*dmats0[6][14] + coeff0_7*dmats0[7][14] + coeff0_8*dmats0[8][14] + coeff0_9*dmats0[9][14] + coeff0_10*dmats0[10][14] + coeff0_11*dmats0[11][14] + coeff0_12*dmats0[12][14] + coeff0_13*dmats0[13][14] + coeff0_14*dmats0[14][14];
1095
 
        }
1096
 
        if(combinations[deriv_num][j] == 1)
1097
 
        {
1098
 
          new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0] + coeff0_6*dmats1[6][0] + coeff0_7*dmats1[7][0] + coeff0_8*dmats1[8][0] + coeff0_9*dmats1[9][0] + coeff0_10*dmats1[10][0] + coeff0_11*dmats1[11][0] + coeff0_12*dmats1[12][0] + coeff0_13*dmats1[13][0] + coeff0_14*dmats1[14][0];
1099
 
          new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1] + coeff0_6*dmats1[6][1] + coeff0_7*dmats1[7][1] + coeff0_8*dmats1[8][1] + coeff0_9*dmats1[9][1] + coeff0_10*dmats1[10][1] + coeff0_11*dmats1[11][1] + coeff0_12*dmats1[12][1] + coeff0_13*dmats1[13][1] + coeff0_14*dmats1[14][1];
1100
 
          new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2] + coeff0_6*dmats1[6][2] + coeff0_7*dmats1[7][2] + coeff0_8*dmats1[8][2] + coeff0_9*dmats1[9][2] + coeff0_10*dmats1[10][2] + coeff0_11*dmats1[11][2] + coeff0_12*dmats1[12][2] + coeff0_13*dmats1[13][2] + coeff0_14*dmats1[14][2];
1101
 
          new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3] + coeff0_6*dmats1[6][3] + coeff0_7*dmats1[7][3] + coeff0_8*dmats1[8][3] + coeff0_9*dmats1[9][3] + coeff0_10*dmats1[10][3] + coeff0_11*dmats1[11][3] + coeff0_12*dmats1[12][3] + coeff0_13*dmats1[13][3] + coeff0_14*dmats1[14][3];
1102
 
          new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4] + coeff0_6*dmats1[6][4] + coeff0_7*dmats1[7][4] + coeff0_8*dmats1[8][4] + coeff0_9*dmats1[9][4] + coeff0_10*dmats1[10][4] + coeff0_11*dmats1[11][4] + coeff0_12*dmats1[12][4] + coeff0_13*dmats1[13][4] + coeff0_14*dmats1[14][4];
1103
 
          new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5] + coeff0_6*dmats1[6][5] + coeff0_7*dmats1[7][5] + coeff0_8*dmats1[8][5] + coeff0_9*dmats1[9][5] + coeff0_10*dmats1[10][5] + coeff0_11*dmats1[11][5] + coeff0_12*dmats1[12][5] + coeff0_13*dmats1[13][5] + coeff0_14*dmats1[14][5];
1104
 
          new_coeff0_6 = coeff0_0*dmats1[0][6] + coeff0_1*dmats1[1][6] + coeff0_2*dmats1[2][6] + coeff0_3*dmats1[3][6] + coeff0_4*dmats1[4][6] + coeff0_5*dmats1[5][6] + coeff0_6*dmats1[6][6] + coeff0_7*dmats1[7][6] + coeff0_8*dmats1[8][6] + coeff0_9*dmats1[9][6] + coeff0_10*dmats1[10][6] + coeff0_11*dmats1[11][6] + coeff0_12*dmats1[12][6] + coeff0_13*dmats1[13][6] + coeff0_14*dmats1[14][6];
1105
 
          new_coeff0_7 = coeff0_0*dmats1[0][7] + coeff0_1*dmats1[1][7] + coeff0_2*dmats1[2][7] + coeff0_3*dmats1[3][7] + coeff0_4*dmats1[4][7] + coeff0_5*dmats1[5][7] + coeff0_6*dmats1[6][7] + coeff0_7*dmats1[7][7] + coeff0_8*dmats1[8][7] + coeff0_9*dmats1[9][7] + coeff0_10*dmats1[10][7] + coeff0_11*dmats1[11][7] + coeff0_12*dmats1[12][7] + coeff0_13*dmats1[13][7] + coeff0_14*dmats1[14][7];
1106
 
          new_coeff0_8 = coeff0_0*dmats1[0][8] + coeff0_1*dmats1[1][8] + coeff0_2*dmats1[2][8] + coeff0_3*dmats1[3][8] + coeff0_4*dmats1[4][8] + coeff0_5*dmats1[5][8] + coeff0_6*dmats1[6][8] + coeff0_7*dmats1[7][8] + coeff0_8*dmats1[8][8] + coeff0_9*dmats1[9][8] + coeff0_10*dmats1[10][8] + coeff0_11*dmats1[11][8] + coeff0_12*dmats1[12][8] + coeff0_13*dmats1[13][8] + coeff0_14*dmats1[14][8];
1107
 
          new_coeff0_9 = coeff0_0*dmats1[0][9] + coeff0_1*dmats1[1][9] + coeff0_2*dmats1[2][9] + coeff0_3*dmats1[3][9] + coeff0_4*dmats1[4][9] + coeff0_5*dmats1[5][9] + coeff0_6*dmats1[6][9] + coeff0_7*dmats1[7][9] + coeff0_8*dmats1[8][9] + coeff0_9*dmats1[9][9] + coeff0_10*dmats1[10][9] + coeff0_11*dmats1[11][9] + coeff0_12*dmats1[12][9] + coeff0_13*dmats1[13][9] + coeff0_14*dmats1[14][9];
1108
 
          new_coeff0_10 = coeff0_0*dmats1[0][10] + coeff0_1*dmats1[1][10] + coeff0_2*dmats1[2][10] + coeff0_3*dmats1[3][10] + coeff0_4*dmats1[4][10] + coeff0_5*dmats1[5][10] + coeff0_6*dmats1[6][10] + coeff0_7*dmats1[7][10] + coeff0_8*dmats1[8][10] + coeff0_9*dmats1[9][10] + coeff0_10*dmats1[10][10] + coeff0_11*dmats1[11][10] + coeff0_12*dmats1[12][10] + coeff0_13*dmats1[13][10] + coeff0_14*dmats1[14][10];
1109
 
          new_coeff0_11 = coeff0_0*dmats1[0][11] + coeff0_1*dmats1[1][11] + coeff0_2*dmats1[2][11] + coeff0_3*dmats1[3][11] + coeff0_4*dmats1[4][11] + coeff0_5*dmats1[5][11] + coeff0_6*dmats1[6][11] + coeff0_7*dmats1[7][11] + coeff0_8*dmats1[8][11] + coeff0_9*dmats1[9][11] + coeff0_10*dmats1[10][11] + coeff0_11*dmats1[11][11] + coeff0_12*dmats1[12][11] + coeff0_13*dmats1[13][11] + coeff0_14*dmats1[14][11];
1110
 
          new_coeff0_12 = coeff0_0*dmats1[0][12] + coeff0_1*dmats1[1][12] + coeff0_2*dmats1[2][12] + coeff0_3*dmats1[3][12] + coeff0_4*dmats1[4][12] + coeff0_5*dmats1[5][12] + coeff0_6*dmats1[6][12] + coeff0_7*dmats1[7][12] + coeff0_8*dmats1[8][12] + coeff0_9*dmats1[9][12] + coeff0_10*dmats1[10][12] + coeff0_11*dmats1[11][12] + coeff0_12*dmats1[12][12] + coeff0_13*dmats1[13][12] + coeff0_14*dmats1[14][12];
1111
 
          new_coeff0_13 = coeff0_0*dmats1[0][13] + coeff0_1*dmats1[1][13] + coeff0_2*dmats1[2][13] + coeff0_3*dmats1[3][13] + coeff0_4*dmats1[4][13] + coeff0_5*dmats1[5][13] + coeff0_6*dmats1[6][13] + coeff0_7*dmats1[7][13] + coeff0_8*dmats1[8][13] + coeff0_9*dmats1[9][13] + coeff0_10*dmats1[10][13] + coeff0_11*dmats1[11][13] + coeff0_12*dmats1[12][13] + coeff0_13*dmats1[13][13] + coeff0_14*dmats1[14][13];
1112
 
          new_coeff0_14 = coeff0_0*dmats1[0][14] + coeff0_1*dmats1[1][14] + coeff0_2*dmats1[2][14] + coeff0_3*dmats1[3][14] + coeff0_4*dmats1[4][14] + coeff0_5*dmats1[5][14] + coeff0_6*dmats1[6][14] + coeff0_7*dmats1[7][14] + coeff0_8*dmats1[8][14] + coeff0_9*dmats1[9][14] + coeff0_10*dmats1[10][14] + coeff0_11*dmats1[11][14] + coeff0_12*dmats1[12][14] + coeff0_13*dmats1[13][14] + coeff0_14*dmats1[14][14];
1113
 
        }
1114
 
    
1115
 
      }
1116
 
      // Compute derivatives on reference element as dot product of coefficients and basisvalues
1117
 
      derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5 + new_coeff0_6*basisvalue6 + new_coeff0_7*basisvalue7 + new_coeff0_8*basisvalue8 + new_coeff0_9*basisvalue9 + new_coeff0_10*basisvalue10 + new_coeff0_11*basisvalue11 + new_coeff0_12*basisvalue12 + new_coeff0_13*basisvalue13 + new_coeff0_14*basisvalue14;
1118
 
    }
1119
 
    
1120
 
    // Transform derivatives back to physical element
1121
 
    for (unsigned int row = 0; row < num_derivatives; row++)
1122
 
    {
1123
 
      for (unsigned int col = 0; col < num_derivatives; col++)
1124
 
      {
1125
 
        values[row] += transform[row][col]*derivatives[col];
1126
 
      }
1127
 
    }
1128
 
    // Delete pointer to array of derivatives on FIAT element
1129
 
    delete [] derivatives;
1130
 
    
1131
 
    // Delete pointer to array of combinations of derivatives and transform
1132
 
    for (unsigned int row = 0; row < num_derivatives; row++)
1133
 
    {
1134
 
      delete [] combinations[row];
1135
 
      delete [] transform[row];
1136
 
    }
1137
 
    
1138
 
    delete [] combinations;
1139
 
    delete [] transform;
1140
 
  }
1141
 
 
1142
 
  /// Evaluate order n derivatives of all basis functions at given point in cell
1143
 
  virtual void evaluate_basis_derivatives_all(unsigned int n,
1144
 
                                              double* values,
1145
 
                                              const double* coordinates,
1146
 
                                              const ufc::cell& c) const
1147
 
  {
1148
 
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
1149
 
  }
1150
 
 
1151
 
  /// Evaluate linear functional for dof i on the function f
1152
 
  virtual double evaluate_dof(unsigned int i,
1153
 
                              const ufc::function& f,
1154
 
                              const ufc::cell& c) const
1155
 
  {
1156
 
    // The reference points, direction and weights:
1157
 
    static const double X[12][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.75, 0.25}}, {{0.5, 0.5}}, {{0.25, 0.75}}, {{0, 0.25}}, {{0, 0.5}}, {{0, 0.75}}, {{0.25, 0}}, {{0.5, 0}}, {{0.75, 0}}};
1158
 
    static const double W[12][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
1159
 
    static const double D[12][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
1160
 
    
1161
 
    const double * const * x = c.coordinates;
1162
 
    double result = 0.0;
1163
 
    // Iterate over the points:
1164
 
    // Evaluate basis functions for affine mapping
1165
 
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
1166
 
    const double w1 = X[i][0][0];
1167
 
    const double w2 = X[i][0][1];
1168
 
    
1169
 
    // Compute affine mapping y = F(X)
1170
 
    double y[2];
1171
 
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
1172
 
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
1173
 
    
1174
 
    // Evaluate function at physical points
1175
 
    double values[1];
1176
 
    f.evaluate(values, y, c);
1177
 
    
1178
 
    // Map function values using appropriate mapping
1179
 
    // Affine map: Do nothing
1180
 
    
1181
 
    // Note that we do not map the weights (yet).
1182
 
    
1183
 
    // Take directional components
1184
 
    for(int k = 0; k < 1; k++)
1185
 
      result += values[k]*D[i][0][k];
1186
 
    // Multiply by weights
1187
 
    result *= W[i][0];
1188
 
    
1189
 
    return result;
1190
 
  }
1191
 
 
1192
 
  /// Evaluate linear functionals for all dofs on the function f
1193
 
  virtual void evaluate_dofs(double* values,
1194
 
                             const ufc::function& f,
1195
 
                             const ufc::cell& c) const
1196
 
  {
1197
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
1198
 
  }
1199
 
 
1200
 
  /// Interpolate vertex values from dof values
1201
 
  virtual void interpolate_vertex_values(double* vertex_values,
1202
 
                                         const double* dof_values,
1203
 
                                         const ufc::cell& c) const
1204
 
  {
1205
 
    // Evaluate at vertices and use affine mapping
1206
 
    vertex_values[0] = dof_values[0];
1207
 
    vertex_values[1] = dof_values[1];
1208
 
    vertex_values[2] = dof_values[2];
1209
 
  }
1210
 
 
1211
 
  /// Return the number of sub elements (for a mixed element)
1212
 
  virtual unsigned int num_sub_elements() const
1213
 
  {
1214
 
    return 1;
1215
 
  }
1216
 
 
1217
 
  /// Create a new finite element for sub element i (for a mixed element)
1218
 
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
1219
 
  {
1220
 
    return new elementrestriction_0_finite_element_0_1();
1221
 
  }
1222
 
 
1223
 
};
1224
 
 
1225
 
/// This class defines the interface for a finite element.
1226
 
 
1227
 
class elementrestriction_0_finite_element_0: public ufc::finite_element
1228
 
{
1229
 
public:
1230
 
 
1231
 
  /// Constructor
1232
 
  elementrestriction_0_finite_element_0() : ufc::finite_element()
1233
 
  {
1234
 
    // Do nothing
1235
 
  }
1236
 
 
1237
 
  /// Destructor
1238
 
  virtual ~elementrestriction_0_finite_element_0()
1239
 
  {
1240
 
    // Do nothing
1241
 
  }
1242
 
 
1243
 
  /// Return a string identifying the finite element
1244
 
  virtual const char* signature() const
1245
 
  {
1246
 
    return "MixedElement(*[FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 4), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 4), Cell('interval', 1, Space(1)))], **{'value_shape': (2,) })";
1247
 
  }
1248
 
 
1249
 
  /// Return the cell shape
1250
 
  virtual ufc::shape cell_shape() const
1251
 
  {
1252
 
    return ufc::triangle;
1253
 
  }
1254
 
 
1255
 
  /// Return the dimension of the finite element function space
1256
 
  virtual unsigned int space_dimension() const
1257
 
  {
1258
 
    return 24;
1259
 
  }
1260
 
 
1261
 
  /// Return the rank of the value space
1262
 
  virtual unsigned int value_rank() const
1263
 
  {
1264
 
    return 1;
1265
 
  }
1266
 
 
1267
 
  /// Return the dimension of the value space for axis i
1268
 
  virtual unsigned int value_dimension(unsigned int i) const
1269
 
  {
1270
 
    return 2;
1271
 
  }
1272
 
 
1273
 
  /// Evaluate basis function i at given point in cell
1274
 
  virtual void evaluate_basis(unsigned int i,
1275
 
                              double* values,
1276
 
                              const double* coordinates,
1277
 
                              const ufc::cell& c) const
1278
 
  {
1279
 
    // Extract vertex coordinates
1280
 
    const double * const * element_coordinates = c.coordinates;
1281
 
    
1282
 
    // Compute Jacobian of affine map from reference cell
1283
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1284
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1285
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1286
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1287
 
    
1288
 
    // Compute determinant of Jacobian
1289
 
    const double detJ = J_00*J_11 - J_01*J_10;
1290
 
    
1291
 
    // Compute inverse of Jacobian
1292
 
    
1293
 
    // Get coordinates and map to the reference (UFC) element
1294
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1295
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
1296
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1297
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1298
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
1299
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1300
 
    
1301
 
    // Map coordinates to the reference square
1302
 
    if (std::abs(y - 1.0) < 1e-08)
1303
 
      x = -1.0;
1304
 
    else
1305
 
      x = 2.0 *x/(1.0 - y) - 1.0;
1306
 
    y = 2.0*y - 1.0;
1307
 
    
1308
 
    // Reset values
1309
 
    values[0] = 0;
1310
 
    values[1] = 0;
1311
 
    
1312
 
    if (0 <= i && i <= 11)
1313
 
    {
1314
 
      // Map degree of freedom to element degree of freedom
1315
 
      const unsigned int dof = i;
1316
 
    
1317
 
      // Generate scalings
1318
 
      const double scalings_y_0 = 1;
1319
 
      const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1320
 
      const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
1321
 
      const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
1322
 
      const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
1323
 
    
1324
 
      // Compute psitilde_a
1325
 
      const double psitilde_a_0 = 1;
1326
 
      const double psitilde_a_1 = x;
1327
 
      const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
1328
 
      const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;
1329
 
      const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
1330
 
    
1331
 
      // Compute psitilde_bs
1332
 
      const double psitilde_bs_0_0 = 1;
1333
 
      const double psitilde_bs_0_1 = 1.5*y + 0.5;
1334
 
      const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
1335
 
      const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
1336
 
      const double psitilde_bs_0_4 = 0.0285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571*psitilde_bs_0_2;
1337
 
      const double psitilde_bs_1_0 = 1;
1338
 
      const double psitilde_bs_1_1 = 2.5*y + 1.5;
1339
 
      const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
1340
 
      const double psitilde_bs_1_3 = 0.285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714*psitilde_bs_1_1;
1341
 
      const double psitilde_bs_2_0 = 1;
1342
 
      const double psitilde_bs_2_1 = 3.5*y + 2.5;
1343
 
      const double psitilde_bs_2_2 = 1.02040816*psitilde_bs_2_1 + 2.57142857*y*psitilde_bs_2_1 - 0.551020408*psitilde_bs_2_0;
1344
 
      const double psitilde_bs_3_0 = 1;
1345
 
      const double psitilde_bs_3_1 = 4.5*y + 3.5;
1346
 
      const double psitilde_bs_4_0 = 1;
1347
 
    
1348
 
      // Compute basisvalues
1349
 
      const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1350
 
      const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1351
 
      const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1352
 
      const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
1353
 
      const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
1354
 
      const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
1355
 
      const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
1356
 
      const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
1357
 
      const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
1358
 
      const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
1359
 
      const double basisvalue10 = 4.74341649*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
1360
 
      const double basisvalue11 = 4.18330013*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
1361
 
      const double basisvalue12 = 3.53553391*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
1362
 
      const double basisvalue13 = 2.73861279*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
1363
 
      const double basisvalue14 = 1.58113883*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
1364
 
    
1365
 
      // Table(s) of coefficients
1366
 
      static const double coefficients0[12][15] =   \
1367
 
      {{0, -0.0412393049, -0.0238095238, 0.0289800295, 0.0224478343, 0.0129602632, -0.0395942581, -0.0334632557, -0.0259205264, -0.0149652229, 0.0321247254, 0.0283313448, 0.0239443566, 0.0185472189, 0.0107082418},
1368
 
      {0, 0.0412393049, -0.0238095238, 0.0289800295, -0.0224478343, 0.0129602632, 0.0395942581, -0.0334632557, 0.0259205264, -0.0149652229, 0.0321247254, -0.0283313448, 0.0239443566, -0.0185472189, 0.0107082418},
1369
 
      {0, 0, 0.0476190476, 0, 0, 0.0388807896, 0, 0, 0, 0.0598608915, 0, 0, 0, 0, 0.0535412091},
1370
 
      {0.125707872, 0.131965776, -0.0253968254, 0.139104142, -0.0718330698, 0.0311046317, 0.0633508129, 0.0267706045, -0.0622092633, 0.0478887132, 0, 0.0566626896, -0.0838052481, 0.083462485, -0.0535412091},
1371
 
      {-0.0314269681, 0.010997148, 0.00634920635, 0, 0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, -0.139104142, 0.107082418},
1372
 
      {0.125707872, 0.0439885919, 0.126984127, 0, 0.0359165349, 0.155523158, 0, 0, 0.103682106, -0.0119721783, 0, 0, 0, 0.0927360944, -0.107082418},
1373
 
      {0.125707872, -0.131965776, -0.0253968254, 0.139104142, 0.0718330698, 0.0311046317, -0.0633508129, 0.0267706045, 0.0622092633, 0.0478887132, 0, -0.0566626896, -0.0838052481, -0.083462485, -0.0535412091},
1374
 
      {-0.0314269681, -0.010997148, 0.00634920635, 0, -0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, 0.139104142, 0.107082418},
1375
 
      {0.125707872, -0.0439885919, 0.126984127, 0, -0.0359165349, 0.155523158, 0, 0, -0.103682106, -0.0119721783, 0, 0, 0, -0.0927360944, -0.107082418},
1376
 
      {0.125707872, -0.0879771839, -0.101587302, 0.0927360944, 0.107749605, 0.0725774739, 0.0791885161, -0.0133853023, -0.0518410528, -0.0419026241, -0.128498902, -0.0566626896, -0.0119721783, 0.00927360944, 0.0107082418},
1377
 
      {-0.0314269681, 0, -0.0126984127, -0.243432248, 0, 0.0544331054, 0, 0.0936971159, 0, -0.0419026241, 0.192748353, 0, -0.0239443566, 0, 0.0107082418},
1378
 
      {0.125707872, 0.0879771839, -0.101587302, 0.0927360944, -0.107749605, 0.0725774739, -0.0791885161, -0.0133853023, 0.0518410528, -0.0419026241, -0.128498902, 0.0566626896, -0.0119721783, -0.00927360944, 0.0107082418}};
1379
 
    
1380
 
      // Extract relevant coefficients
1381
 
      const double coeff0_0 =   coefficients0[dof][0];
1382
 
      const double coeff0_1 =   coefficients0[dof][1];
1383
 
      const double coeff0_2 =   coefficients0[dof][2];
1384
 
      const double coeff0_3 =   coefficients0[dof][3];
1385
 
      const double coeff0_4 =   coefficients0[dof][4];
1386
 
      const double coeff0_5 =   coefficients0[dof][5];
1387
 
      const double coeff0_6 =   coefficients0[dof][6];
1388
 
      const double coeff0_7 =   coefficients0[dof][7];
1389
 
      const double coeff0_8 =   coefficients0[dof][8];
1390
 
      const double coeff0_9 =   coefficients0[dof][9];
1391
 
      const double coeff0_10 =   coefficients0[dof][10];
1392
 
      const double coeff0_11 =   coefficients0[dof][11];
1393
 
      const double coeff0_12 =   coefficients0[dof][12];
1394
 
      const double coeff0_13 =   coefficients0[dof][13];
1395
 
      const double coeff0_14 =   coefficients0[dof][14];
1396
 
    
1397
 
      // Compute value(s)
1398
 
      values[0] = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5 + coeff0_6*basisvalue6 + coeff0_7*basisvalue7 + coeff0_8*basisvalue8 + coeff0_9*basisvalue9 + coeff0_10*basisvalue10 + coeff0_11*basisvalue11 + coeff0_12*basisvalue12 + coeff0_13*basisvalue13 + coeff0_14*basisvalue14;
1399
 
    }
1400
 
    
1401
 
    if (12 <= i && i <= 23)
1402
 
    {
1403
 
      // Map degree of freedom to element degree of freedom
1404
 
      const unsigned int dof = i - 12;
1405
 
    
1406
 
      // Generate scalings
1407
 
      const double scalings_y_0 = 1;
1408
 
      const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1409
 
      const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
1410
 
      const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
1411
 
      const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
1412
 
    
1413
 
      // Compute psitilde_a
1414
 
      const double psitilde_a_0 = 1;
1415
 
      const double psitilde_a_1 = x;
1416
 
      const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
1417
 
      const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;
1418
 
      const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
1419
 
    
1420
 
      // Compute psitilde_bs
1421
 
      const double psitilde_bs_0_0 = 1;
1422
 
      const double psitilde_bs_0_1 = 1.5*y + 0.5;
1423
 
      const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
1424
 
      const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
1425
 
      const double psitilde_bs_0_4 = 0.0285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571*psitilde_bs_0_2;
1426
 
      const double psitilde_bs_1_0 = 1;
1427
 
      const double psitilde_bs_1_1 = 2.5*y + 1.5;
1428
 
      const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
1429
 
      const double psitilde_bs_1_3 = 0.285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714*psitilde_bs_1_1;
1430
 
      const double psitilde_bs_2_0 = 1;
1431
 
      const double psitilde_bs_2_1 = 3.5*y + 2.5;
1432
 
      const double psitilde_bs_2_2 = 1.02040816*psitilde_bs_2_1 + 2.57142857*y*psitilde_bs_2_1 - 0.551020408*psitilde_bs_2_0;
1433
 
      const double psitilde_bs_3_0 = 1;
1434
 
      const double psitilde_bs_3_1 = 4.5*y + 3.5;
1435
 
      const double psitilde_bs_4_0 = 1;
1436
 
    
1437
 
      // Compute basisvalues
1438
 
      const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1439
 
      const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1440
 
      const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1441
 
      const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
1442
 
      const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
1443
 
      const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
1444
 
      const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
1445
 
      const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
1446
 
      const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
1447
 
      const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
1448
 
      const double basisvalue10 = 4.74341649*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
1449
 
      const double basisvalue11 = 4.18330013*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
1450
 
      const double basisvalue12 = 3.53553391*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
1451
 
      const double basisvalue13 = 2.73861279*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
1452
 
      const double basisvalue14 = 1.58113883*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
1453
 
    
1454
 
      // Table(s) of coefficients
1455
 
      static const double coefficients0[12][15] =   \
1456
 
      {{0, -0.0412393049, -0.0238095238, 0.0289800295, 0.0224478343, 0.0129602632, -0.0395942581, -0.0334632557, -0.0259205264, -0.0149652229, 0.0321247254, 0.0283313448, 0.0239443566, 0.0185472189, 0.0107082418},
1457
 
      {0, 0.0412393049, -0.0238095238, 0.0289800295, -0.0224478343, 0.0129602632, 0.0395942581, -0.0334632557, 0.0259205264, -0.0149652229, 0.0321247254, -0.0283313448, 0.0239443566, -0.0185472189, 0.0107082418},
1458
 
      {0, 0, 0.0476190476, 0, 0, 0.0388807896, 0, 0, 0, 0.0598608915, 0, 0, 0, 0, 0.0535412091},
1459
 
      {0.125707872, 0.131965776, -0.0253968254, 0.139104142, -0.0718330698, 0.0311046317, 0.0633508129, 0.0267706045, -0.0622092633, 0.0478887132, 0, 0.0566626896, -0.0838052481, 0.083462485, -0.0535412091},
1460
 
      {-0.0314269681, 0.010997148, 0.00634920635, 0, 0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, -0.139104142, 0.107082418},
1461
 
      {0.125707872, 0.0439885919, 0.126984127, 0, 0.0359165349, 0.155523158, 0, 0, 0.103682106, -0.0119721783, 0, 0, 0, 0.0927360944, -0.107082418},
1462
 
      {0.125707872, -0.131965776, -0.0253968254, 0.139104142, 0.0718330698, 0.0311046317, -0.0633508129, 0.0267706045, 0.0622092633, 0.0478887132, 0, -0.0566626896, -0.0838052481, -0.083462485, -0.0535412091},
1463
 
      {-0.0314269681, -0.010997148, 0.00634920635, 0, -0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, 0.139104142, 0.107082418},
1464
 
      {0.125707872, -0.0439885919, 0.126984127, 0, -0.0359165349, 0.155523158, 0, 0, -0.103682106, -0.0119721783, 0, 0, 0, -0.0927360944, -0.107082418},
1465
 
      {0.125707872, -0.0879771839, -0.101587302, 0.0927360944, 0.107749605, 0.0725774739, 0.0791885161, -0.0133853023, -0.0518410528, -0.0419026241, -0.128498902, -0.0566626896, -0.0119721783, 0.00927360944, 0.0107082418},
1466
 
      {-0.0314269681, 0, -0.0126984127, -0.243432248, 0, 0.0544331054, 0, 0.0936971159, 0, -0.0419026241, 0.192748353, 0, -0.0239443566, 0, 0.0107082418},
1467
 
      {0.125707872, 0.0879771839, -0.101587302, 0.0927360944, -0.107749605, 0.0725774739, -0.0791885161, -0.0133853023, 0.0518410528, -0.0419026241, -0.128498902, 0.0566626896, -0.0119721783, -0.00927360944, 0.0107082418}};
1468
 
    
1469
 
      // Extract relevant coefficients
1470
 
      const double coeff0_0 =   coefficients0[dof][0];
1471
 
      const double coeff0_1 =   coefficients0[dof][1];
1472
 
      const double coeff0_2 =   coefficients0[dof][2];
1473
 
      const double coeff0_3 =   coefficients0[dof][3];
1474
 
      const double coeff0_4 =   coefficients0[dof][4];
1475
 
      const double coeff0_5 =   coefficients0[dof][5];
1476
 
      const double coeff0_6 =   coefficients0[dof][6];
1477
 
      const double coeff0_7 =   coefficients0[dof][7];
1478
 
      const double coeff0_8 =   coefficients0[dof][8];
1479
 
      const double coeff0_9 =   coefficients0[dof][9];
1480
 
      const double coeff0_10 =   coefficients0[dof][10];
1481
 
      const double coeff0_11 =   coefficients0[dof][11];
1482
 
      const double coeff0_12 =   coefficients0[dof][12];
1483
 
      const double coeff0_13 =   coefficients0[dof][13];
1484
 
      const double coeff0_14 =   coefficients0[dof][14];
1485
 
    
1486
 
      // Compute value(s)
1487
 
      values[1] = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5 + coeff0_6*basisvalue6 + coeff0_7*basisvalue7 + coeff0_8*basisvalue8 + coeff0_9*basisvalue9 + coeff0_10*basisvalue10 + coeff0_11*basisvalue11 + coeff0_12*basisvalue12 + coeff0_13*basisvalue13 + coeff0_14*basisvalue14;
1488
 
    }
1489
 
    
1490
 
  }
1491
 
 
1492
 
  /// Evaluate all basis functions at given point in cell
1493
 
  virtual void evaluate_basis_all(double* values,
1494
 
                                  const double* coordinates,
1495
 
                                  const ufc::cell& c) const
1496
 
  {
1497
 
    throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
1498
 
  }
1499
 
 
1500
 
  /// Evaluate order n derivatives of basis function i at given point in cell
1501
 
  virtual void evaluate_basis_derivatives(unsigned int i,
1502
 
                                          unsigned int n,
1503
 
                                          double* values,
1504
 
                                          const double* coordinates,
1505
 
                                          const ufc::cell& c) const
1506
 
  {
1507
 
    // Extract vertex coordinates
1508
 
    const double * const * element_coordinates = c.coordinates;
1509
 
    
1510
 
    // Compute Jacobian of affine map from reference cell
1511
 
    const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
1512
 
    const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
1513
 
    const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
1514
 
    const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
1515
 
    
1516
 
    // Compute determinant of Jacobian
1517
 
    const double detJ = J_00*J_11 - J_01*J_10;
1518
 
    
1519
 
    // Compute inverse of Jacobian
1520
 
    
1521
 
    // Get coordinates and map to the reference (UFC) element
1522
 
    double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
1523
 
                element_coordinates[0][0]*element_coordinates[2][1] +\
1524
 
                J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
1525
 
    double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
1526
 
                element_coordinates[1][0]*element_coordinates[0][1] -\
1527
 
                J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
1528
 
    
1529
 
    // Map coordinates to the reference square
1530
 
    if (std::abs(y - 1.0) < 1e-08)
1531
 
      x = -1.0;
1532
 
    else
1533
 
      x = 2.0 *x/(1.0 - y) - 1.0;
1534
 
    y = 2.0*y - 1.0;
1535
 
    
1536
 
    // Compute number of derivatives
1537
 
    unsigned int num_derivatives = 1;
1538
 
    
1539
 
    for (unsigned int j = 0; j < n; j++)
1540
 
      num_derivatives *= 2;
1541
 
    
1542
 
    
1543
 
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1544
 
    unsigned int **combinations = new unsigned int *[num_derivatives];
1545
 
    
1546
 
    for (unsigned int j = 0; j < num_derivatives; j++)
1547
 
    {
1548
 
      combinations[j] = new unsigned int [n];
1549
 
      for (unsigned int k = 0; k < n; k++)
1550
 
        combinations[j][k] = 0;
1551
 
    }
1552
 
    
1553
 
    // Generate combinations of derivatives
1554
 
    for (unsigned int row = 1; row < num_derivatives; row++)
1555
 
    {
1556
 
      for (unsigned int num = 0; num < row; num++)
1557
 
      {
1558
 
        for (unsigned int col = n-1; col+1 > 0; col--)
1559
 
        {
1560
 
          if (combinations[row][col] + 1 > 1)
1561
 
            combinations[row][col] = 0;
1562
 
          else
1563
 
          {
1564
 
            combinations[row][col] += 1;
1565
 
            break;
1566
 
          }
1567
 
        }
1568
 
      }
1569
 
    }
1570
 
    
1571
 
    // Compute inverse of Jacobian
1572
 
    const double Jinv[2][2] =  {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
1573
 
    
1574
 
    // Declare transformation matrix
1575
 
    // Declare pointer to two dimensional array and initialise
1576
 
    double **transform = new double *[num_derivatives];
1577
 
    
1578
 
    for (unsigned int j = 0; j < num_derivatives; j++)
1579
 
    {
1580
 
      transform[j] = new double [num_derivatives];
1581
 
      for (unsigned int k = 0; k < num_derivatives; k++)
1582
 
        transform[j][k] = 1;
1583
 
    }
1584
 
    
1585
 
    // Construct transformation matrix
1586
 
    for (unsigned int row = 0; row < num_derivatives; row++)
1587
 
    {
1588
 
      for (unsigned int col = 0; col < num_derivatives; col++)
1589
 
      {
1590
 
        for (unsigned int k = 0; k < n; k++)
1591
 
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1592
 
      }
1593
 
    }
1594
 
    
1595
 
    // Reset values
1596
 
    for (unsigned int j = 0; j < 2*num_derivatives; j++)
1597
 
      values[j] = 0;
1598
 
    
1599
 
    if (0 <= i && i <= 11)
1600
 
    {
1601
 
      // Map degree of freedom to element degree of freedom
1602
 
      const unsigned int dof = i;
1603
 
    
1604
 
      // Generate scalings
1605
 
      const double scalings_y_0 = 1;
1606
 
      const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1607
 
      const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
1608
 
      const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
1609
 
      const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
1610
 
    
1611
 
      // Compute psitilde_a
1612
 
      const double psitilde_a_0 = 1;
1613
 
      const double psitilde_a_1 = x;
1614
 
      const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
1615
 
      const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;
1616
 
      const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
1617
 
    
1618
 
      // Compute psitilde_bs
1619
 
      const double psitilde_bs_0_0 = 1;
1620
 
      const double psitilde_bs_0_1 = 1.5*y + 0.5;
1621
 
      const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
1622
 
      const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
1623
 
      const double psitilde_bs_0_4 = 0.0285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571*psitilde_bs_0_2;
1624
 
      const double psitilde_bs_1_0 = 1;
1625
 
      const double psitilde_bs_1_1 = 2.5*y + 1.5;
1626
 
      const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
1627
 
      const double psitilde_bs_1_3 = 0.285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714*psitilde_bs_1_1;
1628
 
      const double psitilde_bs_2_0 = 1;
1629
 
      const double psitilde_bs_2_1 = 3.5*y + 2.5;
1630
 
      const double psitilde_bs_2_2 = 1.02040816*psitilde_bs_2_1 + 2.57142857*y*psitilde_bs_2_1 - 0.551020408*psitilde_bs_2_0;
1631
 
      const double psitilde_bs_3_0 = 1;
1632
 
      const double psitilde_bs_3_1 = 4.5*y + 3.5;
1633
 
      const double psitilde_bs_4_0 = 1;
1634
 
    
1635
 
      // Compute basisvalues
1636
 
      const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1637
 
      const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1638
 
      const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1639
 
      const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
1640
 
      const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
1641
 
      const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
1642
 
      const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
1643
 
      const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
1644
 
      const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
1645
 
      const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
1646
 
      const double basisvalue10 = 4.74341649*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
1647
 
      const double basisvalue11 = 4.18330013*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
1648
 
      const double basisvalue12 = 3.53553391*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
1649
 
      const double basisvalue13 = 2.73861279*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
1650
 
      const double basisvalue14 = 1.58113883*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
1651
 
    
1652
 
      // Table(s) of coefficients
1653
 
      static const double coefficients0[12][15] =   \
1654
 
      {{0, -0.0412393049, -0.0238095238, 0.0289800295, 0.0224478343, 0.0129602632, -0.0395942581, -0.0334632557, -0.0259205264, -0.0149652229, 0.0321247254, 0.0283313448, 0.0239443566, 0.0185472189, 0.0107082418},
1655
 
      {0, 0.0412393049, -0.0238095238, 0.0289800295, -0.0224478343, 0.0129602632, 0.0395942581, -0.0334632557, 0.0259205264, -0.0149652229, 0.0321247254, -0.0283313448, 0.0239443566, -0.0185472189, 0.0107082418},
1656
 
      {0, 0, 0.0476190476, 0, 0, 0.0388807896, 0, 0, 0, 0.0598608915, 0, 0, 0, 0, 0.0535412091},
1657
 
      {0.125707872, 0.131965776, -0.0253968254, 0.139104142, -0.0718330698, 0.0311046317, 0.0633508129, 0.0267706045, -0.0622092633, 0.0478887132, 0, 0.0566626896, -0.0838052481, 0.083462485, -0.0535412091},
1658
 
      {-0.0314269681, 0.010997148, 0.00634920635, 0, 0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, -0.139104142, 0.107082418},
1659
 
      {0.125707872, 0.0439885919, 0.126984127, 0, 0.0359165349, 0.155523158, 0, 0, 0.103682106, -0.0119721783, 0, 0, 0, 0.0927360944, -0.107082418},
1660
 
      {0.125707872, -0.131965776, -0.0253968254, 0.139104142, 0.0718330698, 0.0311046317, -0.0633508129, 0.0267706045, 0.0622092633, 0.0478887132, 0, -0.0566626896, -0.0838052481, -0.083462485, -0.0535412091},
1661
 
      {-0.0314269681, -0.010997148, 0.00634920635, 0, -0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, 0.139104142, 0.107082418},
1662
 
      {0.125707872, -0.0439885919, 0.126984127, 0, -0.0359165349, 0.155523158, 0, 0, -0.103682106, -0.0119721783, 0, 0, 0, -0.0927360944, -0.107082418},
1663
 
      {0.125707872, -0.0879771839, -0.101587302, 0.0927360944, 0.107749605, 0.0725774739, 0.0791885161, -0.0133853023, -0.0518410528, -0.0419026241, -0.128498902, -0.0566626896, -0.0119721783, 0.00927360944, 0.0107082418},
1664
 
      {-0.0314269681, 0, -0.0126984127, -0.243432248, 0, 0.0544331054, 0, 0.0936971159, 0, -0.0419026241, 0.192748353, 0, -0.0239443566, 0, 0.0107082418},
1665
 
      {0.125707872, 0.0879771839, -0.101587302, 0.0927360944, -0.107749605, 0.0725774739, -0.0791885161, -0.0133853023, 0.0518410528, -0.0419026241, -0.128498902, 0.0566626896, -0.0119721783, -0.00927360944, 0.0107082418}};
1666
 
    
1667
 
      // Interesting (new) part
1668
 
      // Tables of derivatives of the polynomial base (transpose)
1669
 
      static const double dmats0[15][15] =   \
1670
 
      {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1671
 
      {4.89897949, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1672
 
      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1673
 
      {0, 9.48683298, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1674
 
      {4, 0, 7.07106781, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1675
 
      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1676
 
      {5.29150262, 0, -2.99332591, 13.662601, 0, 0.611010093, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1677
 
      {0, 4.38178046, 0, 0, 12.5219807, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1678
 
      {3.46410162, 0, 7.83836718, 0, 0, 8.4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1679
 
      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1680
 
      {0, 10.9544512, 0, 0, -3.83325939, 0, 17.7482393, 0, 0.553283335, 0, 0, 0, 0, 0, 0},
1681
 
      {4.73286383, 0, 3.34664011, 4.3643578, 0, -5.07468038, 0, 17.0084013, 0, 1.52127766, 0, 0, 0, 0, 0},
1682
 
      {0, 2.44948974, 0, 0, 9.14285714, 0, 0, 0, 14.8461498, 0, 0, 0, 0, 0, 0},
1683
 
      {3.09838668, 0, 7.66811581, 0, 0, 10.7331263, 0, 0, 0, 9.29516003, 0, 0, 0, 0, 0},
1684
 
      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
1685
 
    
1686
 
      static const double dmats1[15][15] =   \
1687
 
      {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1688
 
      {2.44948974, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1689
 
      {4.24264069, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1690
 
      {2.5819889, 4.74341649, -0.912870929, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1691
 
      {2, 6.12372436, 3.53553391, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1692
 
      {-2.30940108, 0, 8.16496581, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1693
 
      {2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.305505046, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1694
 
      {2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1695
 
      {1.73205081, -5.09116882, 3.91918359, 0, 9.69948452, 4.2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1696
 
      {5, 0, -2.82842712, 0, 0, 12.1243557, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1697
 
      {2.68328157, 5.47722558, -1.8973666, 7.42307489, -1.91662969, 0.663940002, 8.87411967, -1.07142857, 0.276641668, -0.0958314847, 0, 0, 0, 0, 0},
1698
 
      {2.36643191, 2.89827535, 1.67332005, 2.1821789, 5.74704893, -2.53734019, 10.0623059, 8.50420064, -2.19577516, 0.760638829, 0, 0, 0, 0, 0},
1699
 
      {2, 1.22474487, 3.53553391, -7.37711114, 4.57142857, 1.6495722, 0, 11.4997782, 7.42307489, -2.57142857, 0, 0, 0, 0, 0},
1700
 
      {1.54919334, 6.64078309, 3.8340579, 0, -6.19677335, 5.36656315, 0, 0, 13.4164079, 4.64758002, 0, 0, 0, 0, 0},
1701
 
      {-3.57770876, 0, 8.85437745, 0, 0, -3.09838668, 0, 0, 0, 16.0996894, 0, 0, 0, 0, 0}};
1702
 
    
1703
 
      // Compute reference derivatives
1704
 
      // Declare pointer to array of derivatives on FIAT element
1705
 
      double *derivatives = new double [num_derivatives];
1706
 
    
1707
 
      // Declare coefficients
1708
 
      double coeff0_0 = 0;
1709
 
      double coeff0_1 = 0;
1710
 
      double coeff0_2 = 0;
1711
 
      double coeff0_3 = 0;
1712
 
      double coeff0_4 = 0;
1713
 
      double coeff0_5 = 0;
1714
 
      double coeff0_6 = 0;
1715
 
      double coeff0_7 = 0;
1716
 
      double coeff0_8 = 0;
1717
 
      double coeff0_9 = 0;
1718
 
      double coeff0_10 = 0;
1719
 
      double coeff0_11 = 0;
1720
 
      double coeff0_12 = 0;
1721
 
      double coeff0_13 = 0;
1722
 
      double coeff0_14 = 0;
1723
 
    
1724
 
      // Declare new coefficients
1725
 
      double new_coeff0_0 = 0;
1726
 
      double new_coeff0_1 = 0;
1727
 
      double new_coeff0_2 = 0;
1728
 
      double new_coeff0_3 = 0;
1729
 
      double new_coeff0_4 = 0;
1730
 
      double new_coeff0_5 = 0;
1731
 
      double new_coeff0_6 = 0;
1732
 
      double new_coeff0_7 = 0;
1733
 
      double new_coeff0_8 = 0;
1734
 
      double new_coeff0_9 = 0;
1735
 
      double new_coeff0_10 = 0;
1736
 
      double new_coeff0_11 = 0;
1737
 
      double new_coeff0_12 = 0;
1738
 
      double new_coeff0_13 = 0;
1739
 
      double new_coeff0_14 = 0;
1740
 
    
1741
 
      // Loop possible derivatives
1742
 
      for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1743
 
      {
1744
 
        // Get values from coefficients array
1745
 
        new_coeff0_0 = coefficients0[dof][0];
1746
 
        new_coeff0_1 = coefficients0[dof][1];
1747
 
        new_coeff0_2 = coefficients0[dof][2];
1748
 
        new_coeff0_3 = coefficients0[dof][3];
1749
 
        new_coeff0_4 = coefficients0[dof][4];
1750
 
        new_coeff0_5 = coefficients0[dof][5];
1751
 
        new_coeff0_6 = coefficients0[dof][6];
1752
 
        new_coeff0_7 = coefficients0[dof][7];
1753
 
        new_coeff0_8 = coefficients0[dof][8];
1754
 
        new_coeff0_9 = coefficients0[dof][9];
1755
 
        new_coeff0_10 = coefficients0[dof][10];
1756
 
        new_coeff0_11 = coefficients0[dof][11];
1757
 
        new_coeff0_12 = coefficients0[dof][12];
1758
 
        new_coeff0_13 = coefficients0[dof][13];
1759
 
        new_coeff0_14 = coefficients0[dof][14];
1760
 
    
1761
 
        // Loop derivative order
1762
 
        for (unsigned int j = 0; j < n; j++)
1763
 
        {
1764
 
          // Update old coefficients
1765
 
          coeff0_0 = new_coeff0_0;
1766
 
          coeff0_1 = new_coeff0_1;
1767
 
          coeff0_2 = new_coeff0_2;
1768
 
          coeff0_3 = new_coeff0_3;
1769
 
          coeff0_4 = new_coeff0_4;
1770
 
          coeff0_5 = new_coeff0_5;
1771
 
          coeff0_6 = new_coeff0_6;
1772
 
          coeff0_7 = new_coeff0_7;
1773
 
          coeff0_8 = new_coeff0_8;
1774
 
          coeff0_9 = new_coeff0_9;
1775
 
          coeff0_10 = new_coeff0_10;
1776
 
          coeff0_11 = new_coeff0_11;
1777
 
          coeff0_12 = new_coeff0_12;
1778
 
          coeff0_13 = new_coeff0_13;
1779
 
          coeff0_14 = new_coeff0_14;
1780
 
    
1781
 
          if(combinations[deriv_num][j] == 0)
1782
 
          {
1783
 
            new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0] + coeff0_6*dmats0[6][0] + coeff0_7*dmats0[7][0] + coeff0_8*dmats0[8][0] + coeff0_9*dmats0[9][0] + coeff0_10*dmats0[10][0] + coeff0_11*dmats0[11][0] + coeff0_12*dmats0[12][0] + coeff0_13*dmats0[13][0] + coeff0_14*dmats0[14][0];
1784
 
            new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1] + coeff0_6*dmats0[6][1] + coeff0_7*dmats0[7][1] + coeff0_8*dmats0[8][1] + coeff0_9*dmats0[9][1] + coeff0_10*dmats0[10][1] + coeff0_11*dmats0[11][1] + coeff0_12*dmats0[12][1] + coeff0_13*dmats0[13][1] + coeff0_14*dmats0[14][1];
1785
 
            new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2] + coeff0_6*dmats0[6][2] + coeff0_7*dmats0[7][2] + coeff0_8*dmats0[8][2] + coeff0_9*dmats0[9][2] + coeff0_10*dmats0[10][2] + coeff0_11*dmats0[11][2] + coeff0_12*dmats0[12][2] + coeff0_13*dmats0[13][2] + coeff0_14*dmats0[14][2];
1786
 
            new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3] + coeff0_6*dmats0[6][3] + coeff0_7*dmats0[7][3] + coeff0_8*dmats0[8][3] + coeff0_9*dmats0[9][3] + coeff0_10*dmats0[10][3] + coeff0_11*dmats0[11][3] + coeff0_12*dmats0[12][3] + coeff0_13*dmats0[13][3] + coeff0_14*dmats0[14][3];
1787
 
            new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4] + coeff0_6*dmats0[6][4] + coeff0_7*dmats0[7][4] + coeff0_8*dmats0[8][4] + coeff0_9*dmats0[9][4] + coeff0_10*dmats0[10][4] + coeff0_11*dmats0[11][4] + coeff0_12*dmats0[12][4] + coeff0_13*dmats0[13][4] + coeff0_14*dmats0[14][4];
1788
 
            new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5] + coeff0_6*dmats0[6][5] + coeff0_7*dmats0[7][5] + coeff0_8*dmats0[8][5] + coeff0_9*dmats0[9][5] + coeff0_10*dmats0[10][5] + coeff0_11*dmats0[11][5] + coeff0_12*dmats0[12][5] + coeff0_13*dmats0[13][5] + coeff0_14*dmats0[14][5];
1789
 
            new_coeff0_6 = coeff0_0*dmats0[0][6] + coeff0_1*dmats0[1][6] + coeff0_2*dmats0[2][6] + coeff0_3*dmats0[3][6] + coeff0_4*dmats0[4][6] + coeff0_5*dmats0[5][6] + coeff0_6*dmats0[6][6] + coeff0_7*dmats0[7][6] + coeff0_8*dmats0[8][6] + coeff0_9*dmats0[9][6] + coeff0_10*dmats0[10][6] + coeff0_11*dmats0[11][6] + coeff0_12*dmats0[12][6] + coeff0_13*dmats0[13][6] + coeff0_14*dmats0[14][6];
1790
 
            new_coeff0_7 = coeff0_0*dmats0[0][7] + coeff0_1*dmats0[1][7] + coeff0_2*dmats0[2][7] + coeff0_3*dmats0[3][7] + coeff0_4*dmats0[4][7] + coeff0_5*dmats0[5][7] + coeff0_6*dmats0[6][7] + coeff0_7*dmats0[7][7] + coeff0_8*dmats0[8][7] + coeff0_9*dmats0[9][7] + coeff0_10*dmats0[10][7] + coeff0_11*dmats0[11][7] + coeff0_12*dmats0[12][7] + coeff0_13*dmats0[13][7] + coeff0_14*dmats0[14][7];
1791
 
            new_coeff0_8 = coeff0_0*dmats0[0][8] + coeff0_1*dmats0[1][8] + coeff0_2*dmats0[2][8] + coeff0_3*dmats0[3][8] + coeff0_4*dmats0[4][8] + coeff0_5*dmats0[5][8] + coeff0_6*dmats0[6][8] + coeff0_7*dmats0[7][8] + coeff0_8*dmats0[8][8] + coeff0_9*dmats0[9][8] + coeff0_10*dmats0[10][8] + coeff0_11*dmats0[11][8] + coeff0_12*dmats0[12][8] + coeff0_13*dmats0[13][8] + coeff0_14*dmats0[14][8];
1792
 
            new_coeff0_9 = coeff0_0*dmats0[0][9] + coeff0_1*dmats0[1][9] + coeff0_2*dmats0[2][9] + coeff0_3*dmats0[3][9] + coeff0_4*dmats0[4][9] + coeff0_5*dmats0[5][9] + coeff0_6*dmats0[6][9] + coeff0_7*dmats0[7][9] + coeff0_8*dmats0[8][9] + coeff0_9*dmats0[9][9] + coeff0_10*dmats0[10][9] + coeff0_11*dmats0[11][9] + coeff0_12*dmats0[12][9] + coeff0_13*dmats0[13][9] + coeff0_14*dmats0[14][9];
1793
 
            new_coeff0_10 = coeff0_0*dmats0[0][10] + coeff0_1*dmats0[1][10] + coeff0_2*dmats0[2][10] + coeff0_3*dmats0[3][10] + coeff0_4*dmats0[4][10] + coeff0_5*dmats0[5][10] + coeff0_6*dmats0[6][10] + coeff0_7*dmats0[7][10] + coeff0_8*dmats0[8][10] + coeff0_9*dmats0[9][10] + coeff0_10*dmats0[10][10] + coeff0_11*dmats0[11][10] + coeff0_12*dmats0[12][10] + coeff0_13*dmats0[13][10] + coeff0_14*dmats0[14][10];
1794
 
            new_coeff0_11 = coeff0_0*dmats0[0][11] + coeff0_1*dmats0[1][11] + coeff0_2*dmats0[2][11] + coeff0_3*dmats0[3][11] + coeff0_4*dmats0[4][11] + coeff0_5*dmats0[5][11] + coeff0_6*dmats0[6][11] + coeff0_7*dmats0[7][11] + coeff0_8*dmats0[8][11] + coeff0_9*dmats0[9][11] + coeff0_10*dmats0[10][11] + coeff0_11*dmats0[11][11] + coeff0_12*dmats0[12][11] + coeff0_13*dmats0[13][11] + coeff0_14*dmats0[14][11];
1795
 
            new_coeff0_12 = coeff0_0*dmats0[0][12] + coeff0_1*dmats0[1][12] + coeff0_2*dmats0[2][12] + coeff0_3*dmats0[3][12] + coeff0_4*dmats0[4][12] + coeff0_5*dmats0[5][12] + coeff0_6*dmats0[6][12] + coeff0_7*dmats0[7][12] + coeff0_8*dmats0[8][12] + coeff0_9*dmats0[9][12] + coeff0_10*dmats0[10][12] + coeff0_11*dmats0[11][12] + coeff0_12*dmats0[12][12] + coeff0_13*dmats0[13][12] + coeff0_14*dmats0[14][12];
1796
 
            new_coeff0_13 = coeff0_0*dmats0[0][13] + coeff0_1*dmats0[1][13] + coeff0_2*dmats0[2][13] + coeff0_3*dmats0[3][13] + coeff0_4*dmats0[4][13] + coeff0_5*dmats0[5][13] + coeff0_6*dmats0[6][13] + coeff0_7*dmats0[7][13] + coeff0_8*dmats0[8][13] + coeff0_9*dmats0[9][13] + coeff0_10*dmats0[10][13] + coeff0_11*dmats0[11][13] + coeff0_12*dmats0[12][13] + coeff0_13*dmats0[13][13] + coeff0_14*dmats0[14][13];
1797
 
            new_coeff0_14 = coeff0_0*dmats0[0][14] + coeff0_1*dmats0[1][14] + coeff0_2*dmats0[2][14] + coeff0_3*dmats0[3][14] + coeff0_4*dmats0[4][14] + coeff0_5*dmats0[5][14] + coeff0_6*dmats0[6][14] + coeff0_7*dmats0[7][14] + coeff0_8*dmats0[8][14] + coeff0_9*dmats0[9][14] + coeff0_10*dmats0[10][14] + coeff0_11*dmats0[11][14] + coeff0_12*dmats0[12][14] + coeff0_13*dmats0[13][14] + coeff0_14*dmats0[14][14];
1798
 
          }
1799
 
          if(combinations[deriv_num][j] == 1)
1800
 
          {
1801
 
            new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0] + coeff0_6*dmats1[6][0] + coeff0_7*dmats1[7][0] + coeff0_8*dmats1[8][0] + coeff0_9*dmats1[9][0] + coeff0_10*dmats1[10][0] + coeff0_11*dmats1[11][0] + coeff0_12*dmats1[12][0] + coeff0_13*dmats1[13][0] + coeff0_14*dmats1[14][0];
1802
 
            new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1] + coeff0_6*dmats1[6][1] + coeff0_7*dmats1[7][1] + coeff0_8*dmats1[8][1] + coeff0_9*dmats1[9][1] + coeff0_10*dmats1[10][1] + coeff0_11*dmats1[11][1] + coeff0_12*dmats1[12][1] + coeff0_13*dmats1[13][1] + coeff0_14*dmats1[14][1];
1803
 
            new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2] + coeff0_6*dmats1[6][2] + coeff0_7*dmats1[7][2] + coeff0_8*dmats1[8][2] + coeff0_9*dmats1[9][2] + coeff0_10*dmats1[10][2] + coeff0_11*dmats1[11][2] + coeff0_12*dmats1[12][2] + coeff0_13*dmats1[13][2] + coeff0_14*dmats1[14][2];
1804
 
            new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3] + coeff0_6*dmats1[6][3] + coeff0_7*dmats1[7][3] + coeff0_8*dmats1[8][3] + coeff0_9*dmats1[9][3] + coeff0_10*dmats1[10][3] + coeff0_11*dmats1[11][3] + coeff0_12*dmats1[12][3] + coeff0_13*dmats1[13][3] + coeff0_14*dmats1[14][3];
1805
 
            new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4] + coeff0_6*dmats1[6][4] + coeff0_7*dmats1[7][4] + coeff0_8*dmats1[8][4] + coeff0_9*dmats1[9][4] + coeff0_10*dmats1[10][4] + coeff0_11*dmats1[11][4] + coeff0_12*dmats1[12][4] + coeff0_13*dmats1[13][4] + coeff0_14*dmats1[14][4];
1806
 
            new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5] + coeff0_6*dmats1[6][5] + coeff0_7*dmats1[7][5] + coeff0_8*dmats1[8][5] + coeff0_9*dmats1[9][5] + coeff0_10*dmats1[10][5] + coeff0_11*dmats1[11][5] + coeff0_12*dmats1[12][5] + coeff0_13*dmats1[13][5] + coeff0_14*dmats1[14][5];
1807
 
            new_coeff0_6 = coeff0_0*dmats1[0][6] + coeff0_1*dmats1[1][6] + coeff0_2*dmats1[2][6] + coeff0_3*dmats1[3][6] + coeff0_4*dmats1[4][6] + coeff0_5*dmats1[5][6] + coeff0_6*dmats1[6][6] + coeff0_7*dmats1[7][6] + coeff0_8*dmats1[8][6] + coeff0_9*dmats1[9][6] + coeff0_10*dmats1[10][6] + coeff0_11*dmats1[11][6] + coeff0_12*dmats1[12][6] + coeff0_13*dmats1[13][6] + coeff0_14*dmats1[14][6];
1808
 
            new_coeff0_7 = coeff0_0*dmats1[0][7] + coeff0_1*dmats1[1][7] + coeff0_2*dmats1[2][7] + coeff0_3*dmats1[3][7] + coeff0_4*dmats1[4][7] + coeff0_5*dmats1[5][7] + coeff0_6*dmats1[6][7] + coeff0_7*dmats1[7][7] + coeff0_8*dmats1[8][7] + coeff0_9*dmats1[9][7] + coeff0_10*dmats1[10][7] + coeff0_11*dmats1[11][7] + coeff0_12*dmats1[12][7] + coeff0_13*dmats1[13][7] + coeff0_14*dmats1[14][7];
1809
 
            new_coeff0_8 = coeff0_0*dmats1[0][8] + coeff0_1*dmats1[1][8] + coeff0_2*dmats1[2][8] + coeff0_3*dmats1[3][8] + coeff0_4*dmats1[4][8] + coeff0_5*dmats1[5][8] + coeff0_6*dmats1[6][8] + coeff0_7*dmats1[7][8] + coeff0_8*dmats1[8][8] + coeff0_9*dmats1[9][8] + coeff0_10*dmats1[10][8] + coeff0_11*dmats1[11][8] + coeff0_12*dmats1[12][8] + coeff0_13*dmats1[13][8] + coeff0_14*dmats1[14][8];
1810
 
            new_coeff0_9 = coeff0_0*dmats1[0][9] + coeff0_1*dmats1[1][9] + coeff0_2*dmats1[2][9] + coeff0_3*dmats1[3][9] + coeff0_4*dmats1[4][9] + coeff0_5*dmats1[5][9] + coeff0_6*dmats1[6][9] + coeff0_7*dmats1[7][9] + coeff0_8*dmats1[8][9] + coeff0_9*dmats1[9][9] + coeff0_10*dmats1[10][9] + coeff0_11*dmats1[11][9] + coeff0_12*dmats1[12][9] + coeff0_13*dmats1[13][9] + coeff0_14*dmats1[14][9];
1811
 
            new_coeff0_10 = coeff0_0*dmats1[0][10] + coeff0_1*dmats1[1][10] + coeff0_2*dmats1[2][10] + coeff0_3*dmats1[3][10] + coeff0_4*dmats1[4][10] + coeff0_5*dmats1[5][10] + coeff0_6*dmats1[6][10] + coeff0_7*dmats1[7][10] + coeff0_8*dmats1[8][10] + coeff0_9*dmats1[9][10] + coeff0_10*dmats1[10][10] + coeff0_11*dmats1[11][10] + coeff0_12*dmats1[12][10] + coeff0_13*dmats1[13][10] + coeff0_14*dmats1[14][10];
1812
 
            new_coeff0_11 = coeff0_0*dmats1[0][11] + coeff0_1*dmats1[1][11] + coeff0_2*dmats1[2][11] + coeff0_3*dmats1[3][11] + coeff0_4*dmats1[4][11] + coeff0_5*dmats1[5][11] + coeff0_6*dmats1[6][11] + coeff0_7*dmats1[7][11] + coeff0_8*dmats1[8][11] + coeff0_9*dmats1[9][11] + coeff0_10*dmats1[10][11] + coeff0_11*dmats1[11][11] + coeff0_12*dmats1[12][11] + coeff0_13*dmats1[13][11] + coeff0_14*dmats1[14][11];
1813
 
            new_coeff0_12 = coeff0_0*dmats1[0][12] + coeff0_1*dmats1[1][12] + coeff0_2*dmats1[2][12] + coeff0_3*dmats1[3][12] + coeff0_4*dmats1[4][12] + coeff0_5*dmats1[5][12] + coeff0_6*dmats1[6][12] + coeff0_7*dmats1[7][12] + coeff0_8*dmats1[8][12] + coeff0_9*dmats1[9][12] + coeff0_10*dmats1[10][12] + coeff0_11*dmats1[11][12] + coeff0_12*dmats1[12][12] + coeff0_13*dmats1[13][12] + coeff0_14*dmats1[14][12];
1814
 
            new_coeff0_13 = coeff0_0*dmats1[0][13] + coeff0_1*dmats1[1][13] + coeff0_2*dmats1[2][13] + coeff0_3*dmats1[3][13] + coeff0_4*dmats1[4][13] + coeff0_5*dmats1[5][13] + coeff0_6*dmats1[6][13] + coeff0_7*dmats1[7][13] + coeff0_8*dmats1[8][13] + coeff0_9*dmats1[9][13] + coeff0_10*dmats1[10][13] + coeff0_11*dmats1[11][13] + coeff0_12*dmats1[12][13] + coeff0_13*dmats1[13][13] + coeff0_14*dmats1[14][13];
1815
 
            new_coeff0_14 = coeff0_0*dmats1[0][14] + coeff0_1*dmats1[1][14] + coeff0_2*dmats1[2][14] + coeff0_3*dmats1[3][14] + coeff0_4*dmats1[4][14] + coeff0_5*dmats1[5][14] + coeff0_6*dmats1[6][14] + coeff0_7*dmats1[7][14] + coeff0_8*dmats1[8][14] + coeff0_9*dmats1[9][14] + coeff0_10*dmats1[10][14] + coeff0_11*dmats1[11][14] + coeff0_12*dmats1[12][14] + coeff0_13*dmats1[13][14] + coeff0_14*dmats1[14][14];
1816
 
          }
1817
 
    
1818
 
        }
1819
 
        // Compute derivatives on reference element as dot product of coefficients and basisvalues
1820
 
        derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5 + new_coeff0_6*basisvalue6 + new_coeff0_7*basisvalue7 + new_coeff0_8*basisvalue8 + new_coeff0_9*basisvalue9 + new_coeff0_10*basisvalue10 + new_coeff0_11*basisvalue11 + new_coeff0_12*basisvalue12 + new_coeff0_13*basisvalue13 + new_coeff0_14*basisvalue14;
1821
 
      }
1822
 
    
1823
 
      // Transform derivatives back to physical element
1824
 
      for (unsigned int row = 0; row < num_derivatives; row++)
1825
 
      {
1826
 
        for (unsigned int col = 0; col < num_derivatives; col++)
1827
 
        {
1828
 
          values[row] += transform[row][col]*derivatives[col];
1829
 
        }
1830
 
      }
1831
 
      // Delete pointer to array of derivatives on FIAT element
1832
 
      delete [] derivatives;
1833
 
    
1834
 
      // Delete pointer to array of combinations of derivatives and transform
1835
 
      for (unsigned int row = 0; row < num_derivatives; row++)
1836
 
      {
1837
 
        delete [] combinations[row];
1838
 
        delete [] transform[row];
1839
 
      }
1840
 
    
1841
 
      delete [] combinations;
1842
 
      delete [] transform;
1843
 
    }
1844
 
    
1845
 
    if (12 <= i && i <= 23)
1846
 
    {
1847
 
      // Map degree of freedom to element degree of freedom
1848
 
      const unsigned int dof = i - 12;
1849
 
    
1850
 
      // Generate scalings
1851
 
      const double scalings_y_0 = 1;
1852
 
      const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
1853
 
      const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
1854
 
      const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
1855
 
      const double scalings_y_4 = scalings_y_3*(0.5 - 0.5*y);
1856
 
    
1857
 
      // Compute psitilde_a
1858
 
      const double psitilde_a_0 = 1;
1859
 
      const double psitilde_a_1 = x;
1860
 
      const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
1861
 
      const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;
1862
 
      const double psitilde_a_4 = 1.75*x*psitilde_a_3 - 0.75*psitilde_a_2;
1863
 
    
1864
 
      // Compute psitilde_bs
1865
 
      const double psitilde_bs_0_0 = 1;
1866
 
      const double psitilde_bs_0_1 = 1.5*y + 0.5;
1867
 
      const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;
1868
 
      const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
1869
 
      const double psitilde_bs_0_4 = 0.0285714286*psitilde_bs_0_3 + 1.8*y*psitilde_bs_0_3 - 0.771428571*psitilde_bs_0_2;
1870
 
      const double psitilde_bs_1_0 = 1;
1871
 
      const double psitilde_bs_1_1 = 2.5*y + 1.5;
1872
 
      const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
1873
 
      const double psitilde_bs_1_3 = 0.285714286*psitilde_bs_1_2 + 2*y*psitilde_bs_1_2 - 0.714285714*psitilde_bs_1_1;
1874
 
      const double psitilde_bs_2_0 = 1;
1875
 
      const double psitilde_bs_2_1 = 3.5*y + 2.5;
1876
 
      const double psitilde_bs_2_2 = 1.02040816*psitilde_bs_2_1 + 2.57142857*y*psitilde_bs_2_1 - 0.551020408*psitilde_bs_2_0;
1877
 
      const double psitilde_bs_3_0 = 1;
1878
 
      const double psitilde_bs_3_1 = 4.5*y + 3.5;
1879
 
      const double psitilde_bs_4_0 = 1;
1880
 
    
1881
 
      // Compute basisvalues
1882
 
      const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
1883
 
      const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
1884
 
      const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
1885
 
      const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
1886
 
      const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
1887
 
      const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
1888
 
      const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
1889
 
      const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
1890
 
      const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
1891
 
      const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
1892
 
      const double basisvalue10 = 4.74341649*psitilde_a_4*scalings_y_4*psitilde_bs_4_0;
1893
 
      const double basisvalue11 = 4.18330013*psitilde_a_3*scalings_y_3*psitilde_bs_3_1;
1894
 
      const double basisvalue12 = 3.53553391*psitilde_a_2*scalings_y_2*psitilde_bs_2_2;
1895
 
      const double basisvalue13 = 2.73861279*psitilde_a_1*scalings_y_1*psitilde_bs_1_3;
1896
 
      const double basisvalue14 = 1.58113883*psitilde_a_0*scalings_y_0*psitilde_bs_0_4;
1897
 
    
1898
 
      // Table(s) of coefficients
1899
 
      static const double coefficients0[12][15] =   \
1900
 
      {{0, -0.0412393049, -0.0238095238, 0.0289800295, 0.0224478343, 0.0129602632, -0.0395942581, -0.0334632557, -0.0259205264, -0.0149652229, 0.0321247254, 0.0283313448, 0.0239443566, 0.0185472189, 0.0107082418},
1901
 
      {0, 0.0412393049, -0.0238095238, 0.0289800295, -0.0224478343, 0.0129602632, 0.0395942581, -0.0334632557, 0.0259205264, -0.0149652229, 0.0321247254, -0.0283313448, 0.0239443566, -0.0185472189, 0.0107082418},
1902
 
      {0, 0, 0.0476190476, 0, 0, 0.0388807896, 0, 0, 0, 0.0598608915, 0, 0, 0, 0, 0.0535412091},
1903
 
      {0.125707872, 0.131965776, -0.0253968254, 0.139104142, -0.0718330698, 0.0311046317, 0.0633508129, 0.0267706045, -0.0622092633, 0.0478887132, 0, 0.0566626896, -0.0838052481, 0.083462485, -0.0535412091},
1904
 
      {-0.0314269681, 0.010997148, 0.00634920635, 0, 0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, -0.139104142, 0.107082418},
1905
 
      {0.125707872, 0.0439885919, 0.126984127, 0, 0.0359165349, 0.155523158, 0, 0, 0.103682106, -0.0119721783, 0, 0, 0, 0.0927360944, -0.107082418},
1906
 
      {0.125707872, -0.131965776, -0.0253968254, 0.139104142, 0.0718330698, 0.0311046317, -0.0633508129, 0.0267706045, 0.0622092633, 0.0478887132, 0, -0.0566626896, -0.0838052481, -0.083462485, -0.0535412091},
1907
 
      {-0.0314269681, -0.010997148, 0.00634920635, 0, -0.188561808, -0.163299316, 0, 0.0936971159, 0, -0.0419026241, 0, 0, 0.0838052481, 0.139104142, 0.107082418},
1908
 
      {0.125707872, -0.0439885919, 0.126984127, 0, -0.0359165349, 0.155523158, 0, 0, -0.103682106, -0.0119721783, 0, 0, 0, -0.0927360944, -0.107082418},
1909
 
      {0.125707872, -0.0879771839, -0.101587302, 0.0927360944, 0.107749605, 0.0725774739, 0.0791885161, -0.0133853023, -0.0518410528, -0.0419026241, -0.128498902, -0.0566626896, -0.0119721783, 0.00927360944, 0.0107082418},
1910
 
      {-0.0314269681, 0, -0.0126984127, -0.243432248, 0, 0.0544331054, 0, 0.0936971159, 0, -0.0419026241, 0.192748353, 0, -0.0239443566, 0, 0.0107082418},
1911
 
      {0.125707872, 0.0879771839, -0.101587302, 0.0927360944, -0.107749605, 0.0725774739, -0.0791885161, -0.0133853023, 0.0518410528, -0.0419026241, -0.128498902, 0.0566626896, -0.0119721783, -0.00927360944, 0.0107082418}};
1912
 
    
1913
 
      // Interesting (new) part
1914
 
      // Tables of derivatives of the polynomial base (transpose)
1915
 
      static const double dmats0[15][15] =   \
1916
 
      {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1917
 
      {4.89897949, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1918
 
      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1919
 
      {0, 9.48683298, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1920
 
      {4, 0, 7.07106781, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1921
 
      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1922
 
      {5.29150262, 0, -2.99332591, 13.662601, 0, 0.611010093, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1923
 
      {0, 4.38178046, 0, 0, 12.5219807, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1924
 
      {3.46410162, 0, 7.83836718, 0, 0, 8.4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1925
 
      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1926
 
      {0, 10.9544512, 0, 0, -3.83325939, 0, 17.7482393, 0, 0.553283335, 0, 0, 0, 0, 0, 0},
1927
 
      {4.73286383, 0, 3.34664011, 4.3643578, 0, -5.07468038, 0, 17.0084013, 0, 1.52127766, 0, 0, 0, 0, 0},
1928
 
      {0, 2.44948974, 0, 0, 9.14285714, 0, 0, 0, 14.8461498, 0, 0, 0, 0, 0, 0},
1929
 
      {3.09838668, 0, 7.66811581, 0, 0, 10.7331263, 0, 0, 0, 9.29516003, 0, 0, 0, 0, 0},
1930
 
      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
1931
 
    
1932
 
      static const double dmats1[15][15] =   \
1933
 
      {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1934
 
      {2.44948974, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1935
 
      {4.24264069, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1936
 
      {2.5819889, 4.74341649, -0.912870929, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1937
 
      {2, 6.12372436, 3.53553391, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1938
 
      {-2.30940108, 0, 8.16496581, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1939
 
      {2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.305505046, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1940
 
      {2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1941
 
      {1.73205081, -5.09116882, 3.91918359, 0, 9.69948452, 4.2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1942
 
      {5, 0, -2.82842712, 0, 0, 12.1243557, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1943
 
      {2.68328157, 5.47722558, -1.8973666, 7.42307489, -1.91662969, 0.663940002, 8.87411967, -1.07142857, 0.276641668, -0.0958314847, 0, 0, 0, 0, 0},
1944
 
      {2.36643191, 2.89827535, 1.67332005, 2.1821789, 5.74704893, -2.53734019, 10.0623059, 8.50420064, -2.19577516, 0.760638829, 0, 0, 0, 0, 0},
1945
 
      {2, 1.22474487, 3.53553391, -7.37711114, 4.57142857, 1.6495722, 0, 11.4997782, 7.42307489, -2.57142857, 0, 0, 0, 0, 0},
1946
 
      {1.54919334, 6.64078309, 3.8340579, 0, -6.19677335, 5.36656315, 0, 0, 13.4164079, 4.64758002, 0, 0, 0, 0, 0},
1947
 
      {-3.57770876, 0, 8.85437745, 0, 0, -3.09838668, 0, 0, 0, 16.0996894, 0, 0, 0, 0, 0}};
1948
 
    
1949
 
      // Compute reference derivatives
1950
 
      // Declare pointer to array of derivatives on FIAT element
1951
 
      double *derivatives = new double [num_derivatives];
1952
 
    
1953
 
      // Declare coefficients
1954
 
      double coeff0_0 = 0;
1955
 
      double coeff0_1 = 0;
1956
 
      double coeff0_2 = 0;
1957
 
      double coeff0_3 = 0;
1958
 
      double coeff0_4 = 0;
1959
 
      double coeff0_5 = 0;
1960
 
      double coeff0_6 = 0;
1961
 
      double coeff0_7 = 0;
1962
 
      double coeff0_8 = 0;
1963
 
      double coeff0_9 = 0;
1964
 
      double coeff0_10 = 0;
1965
 
      double coeff0_11 = 0;
1966
 
      double coeff0_12 = 0;
1967
 
      double coeff0_13 = 0;
1968
 
      double coeff0_14 = 0;
1969
 
    
1970
 
      // Declare new coefficients
1971
 
      double new_coeff0_0 = 0;
1972
 
      double new_coeff0_1 = 0;
1973
 
      double new_coeff0_2 = 0;
1974
 
      double new_coeff0_3 = 0;
1975
 
      double new_coeff0_4 = 0;
1976
 
      double new_coeff0_5 = 0;
1977
 
      double new_coeff0_6 = 0;
1978
 
      double new_coeff0_7 = 0;
1979
 
      double new_coeff0_8 = 0;
1980
 
      double new_coeff0_9 = 0;
1981
 
      double new_coeff0_10 = 0;
1982
 
      double new_coeff0_11 = 0;
1983
 
      double new_coeff0_12 = 0;
1984
 
      double new_coeff0_13 = 0;
1985
 
      double new_coeff0_14 = 0;
1986
 
    
1987
 
      // Loop possible derivatives
1988
 
      for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
1989
 
      {
1990
 
        // Get values from coefficients array
1991
 
        new_coeff0_0 = coefficients0[dof][0];
1992
 
        new_coeff0_1 = coefficients0[dof][1];
1993
 
        new_coeff0_2 = coefficients0[dof][2];
1994
 
        new_coeff0_3 = coefficients0[dof][3];
1995
 
        new_coeff0_4 = coefficients0[dof][4];
1996
 
        new_coeff0_5 = coefficients0[dof][5];
1997
 
        new_coeff0_6 = coefficients0[dof][6];
1998
 
        new_coeff0_7 = coefficients0[dof][7];
1999
 
        new_coeff0_8 = coefficients0[dof][8];
2000
 
        new_coeff0_9 = coefficients0[dof][9];
2001
 
        new_coeff0_10 = coefficients0[dof][10];
2002
 
        new_coeff0_11 = coefficients0[dof][11];
2003
 
        new_coeff0_12 = coefficients0[dof][12];
2004
 
        new_coeff0_13 = coefficients0[dof][13];
2005
 
        new_coeff0_14 = coefficients0[dof][14];
2006
 
    
2007
 
        // Loop derivative order
2008
 
        for (unsigned int j = 0; j < n; j++)
2009
 
        {
2010
 
          // Update old coefficients
2011
 
          coeff0_0 = new_coeff0_0;
2012
 
          coeff0_1 = new_coeff0_1;
2013
 
          coeff0_2 = new_coeff0_2;
2014
 
          coeff0_3 = new_coeff0_3;
2015
 
          coeff0_4 = new_coeff0_4;
2016
 
          coeff0_5 = new_coeff0_5;
2017
 
          coeff0_6 = new_coeff0_6;
2018
 
          coeff0_7 = new_coeff0_7;
2019
 
          coeff0_8 = new_coeff0_8;
2020
 
          coeff0_9 = new_coeff0_9;
2021
 
          coeff0_10 = new_coeff0_10;
2022
 
          coeff0_11 = new_coeff0_11;
2023
 
          coeff0_12 = new_coeff0_12;
2024
 
          coeff0_13 = new_coeff0_13;
2025
 
          coeff0_14 = new_coeff0_14;
2026
 
    
2027
 
          if(combinations[deriv_num][j] == 0)
2028
 
          {
2029
 
            new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0] + coeff0_6*dmats0[6][0] + coeff0_7*dmats0[7][0] + coeff0_8*dmats0[8][0] + coeff0_9*dmats0[9][0] + coeff0_10*dmats0[10][0] + coeff0_11*dmats0[11][0] + coeff0_12*dmats0[12][0] + coeff0_13*dmats0[13][0] + coeff0_14*dmats0[14][0];
2030
 
            new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1] + coeff0_6*dmats0[6][1] + coeff0_7*dmats0[7][1] + coeff0_8*dmats0[8][1] + coeff0_9*dmats0[9][1] + coeff0_10*dmats0[10][1] + coeff0_11*dmats0[11][1] + coeff0_12*dmats0[12][1] + coeff0_13*dmats0[13][1] + coeff0_14*dmats0[14][1];
2031
 
            new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2] + coeff0_6*dmats0[6][2] + coeff0_7*dmats0[7][2] + coeff0_8*dmats0[8][2] + coeff0_9*dmats0[9][2] + coeff0_10*dmats0[10][2] + coeff0_11*dmats0[11][2] + coeff0_12*dmats0[12][2] + coeff0_13*dmats0[13][2] + coeff0_14*dmats0[14][2];
2032
 
            new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3] + coeff0_6*dmats0[6][3] + coeff0_7*dmats0[7][3] + coeff0_8*dmats0[8][3] + coeff0_9*dmats0[9][3] + coeff0_10*dmats0[10][3] + coeff0_11*dmats0[11][3] + coeff0_12*dmats0[12][3] + coeff0_13*dmats0[13][3] + coeff0_14*dmats0[14][3];
2033
 
            new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4] + coeff0_6*dmats0[6][4] + coeff0_7*dmats0[7][4] + coeff0_8*dmats0[8][4] + coeff0_9*dmats0[9][4] + coeff0_10*dmats0[10][4] + coeff0_11*dmats0[11][4] + coeff0_12*dmats0[12][4] + coeff0_13*dmats0[13][4] + coeff0_14*dmats0[14][4];
2034
 
            new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5] + coeff0_6*dmats0[6][5] + coeff0_7*dmats0[7][5] + coeff0_8*dmats0[8][5] + coeff0_9*dmats0[9][5] + coeff0_10*dmats0[10][5] + coeff0_11*dmats0[11][5] + coeff0_12*dmats0[12][5] + coeff0_13*dmats0[13][5] + coeff0_14*dmats0[14][5];
2035
 
            new_coeff0_6 = coeff0_0*dmats0[0][6] + coeff0_1*dmats0[1][6] + coeff0_2*dmats0[2][6] + coeff0_3*dmats0[3][6] + coeff0_4*dmats0[4][6] + coeff0_5*dmats0[5][6] + coeff0_6*dmats0[6][6] + coeff0_7*dmats0[7][6] + coeff0_8*dmats0[8][6] + coeff0_9*dmats0[9][6] + coeff0_10*dmats0[10][6] + coeff0_11*dmats0[11][6] + coeff0_12*dmats0[12][6] + coeff0_13*dmats0[13][6] + coeff0_14*dmats0[14][6];
2036
 
            new_coeff0_7 = coeff0_0*dmats0[0][7] + coeff0_1*dmats0[1][7] + coeff0_2*dmats0[2][7] + coeff0_3*dmats0[3][7] + coeff0_4*dmats0[4][7] + coeff0_5*dmats0[5][7] + coeff0_6*dmats0[6][7] + coeff0_7*dmats0[7][7] + coeff0_8*dmats0[8][7] + coeff0_9*dmats0[9][7] + coeff0_10*dmats0[10][7] + coeff0_11*dmats0[11][7] + coeff0_12*dmats0[12][7] + coeff0_13*dmats0[13][7] + coeff0_14*dmats0[14][7];
2037
 
            new_coeff0_8 = coeff0_0*dmats0[0][8] + coeff0_1*dmats0[1][8] + coeff0_2*dmats0[2][8] + coeff0_3*dmats0[3][8] + coeff0_4*dmats0[4][8] + coeff0_5*dmats0[5][8] + coeff0_6*dmats0[6][8] + coeff0_7*dmats0[7][8] + coeff0_8*dmats0[8][8] + coeff0_9*dmats0[9][8] + coeff0_10*dmats0[10][8] + coeff0_11*dmats0[11][8] + coeff0_12*dmats0[12][8] + coeff0_13*dmats0[13][8] + coeff0_14*dmats0[14][8];
2038
 
            new_coeff0_9 = coeff0_0*dmats0[0][9] + coeff0_1*dmats0[1][9] + coeff0_2*dmats0[2][9] + coeff0_3*dmats0[3][9] + coeff0_4*dmats0[4][9] + coeff0_5*dmats0[5][9] + coeff0_6*dmats0[6][9] + coeff0_7*dmats0[7][9] + coeff0_8*dmats0[8][9] + coeff0_9*dmats0[9][9] + coeff0_10*dmats0[10][9] + coeff0_11*dmats0[11][9] + coeff0_12*dmats0[12][9] + coeff0_13*dmats0[13][9] + coeff0_14*dmats0[14][9];
2039
 
            new_coeff0_10 = coeff0_0*dmats0[0][10] + coeff0_1*dmats0[1][10] + coeff0_2*dmats0[2][10] + coeff0_3*dmats0[3][10] + coeff0_4*dmats0[4][10] + coeff0_5*dmats0[5][10] + coeff0_6*dmats0[6][10] + coeff0_7*dmats0[7][10] + coeff0_8*dmats0[8][10] + coeff0_9*dmats0[9][10] + coeff0_10*dmats0[10][10] + coeff0_11*dmats0[11][10] + coeff0_12*dmats0[12][10] + coeff0_13*dmats0[13][10] + coeff0_14*dmats0[14][10];
2040
 
            new_coeff0_11 = coeff0_0*dmats0[0][11] + coeff0_1*dmats0[1][11] + coeff0_2*dmats0[2][11] + coeff0_3*dmats0[3][11] + coeff0_4*dmats0[4][11] + coeff0_5*dmats0[5][11] + coeff0_6*dmats0[6][11] + coeff0_7*dmats0[7][11] + coeff0_8*dmats0[8][11] + coeff0_9*dmats0[9][11] + coeff0_10*dmats0[10][11] + coeff0_11*dmats0[11][11] + coeff0_12*dmats0[12][11] + coeff0_13*dmats0[13][11] + coeff0_14*dmats0[14][11];
2041
 
            new_coeff0_12 = coeff0_0*dmats0[0][12] + coeff0_1*dmats0[1][12] + coeff0_2*dmats0[2][12] + coeff0_3*dmats0[3][12] + coeff0_4*dmats0[4][12] + coeff0_5*dmats0[5][12] + coeff0_6*dmats0[6][12] + coeff0_7*dmats0[7][12] + coeff0_8*dmats0[8][12] + coeff0_9*dmats0[9][12] + coeff0_10*dmats0[10][12] + coeff0_11*dmats0[11][12] + coeff0_12*dmats0[12][12] + coeff0_13*dmats0[13][12] + coeff0_14*dmats0[14][12];
2042
 
            new_coeff0_13 = coeff0_0*dmats0[0][13] + coeff0_1*dmats0[1][13] + coeff0_2*dmats0[2][13] + coeff0_3*dmats0[3][13] + coeff0_4*dmats0[4][13] + coeff0_5*dmats0[5][13] + coeff0_6*dmats0[6][13] + coeff0_7*dmats0[7][13] + coeff0_8*dmats0[8][13] + coeff0_9*dmats0[9][13] + coeff0_10*dmats0[10][13] + coeff0_11*dmats0[11][13] + coeff0_12*dmats0[12][13] + coeff0_13*dmats0[13][13] + coeff0_14*dmats0[14][13];
2043
 
            new_coeff0_14 = coeff0_0*dmats0[0][14] + coeff0_1*dmats0[1][14] + coeff0_2*dmats0[2][14] + coeff0_3*dmats0[3][14] + coeff0_4*dmats0[4][14] + coeff0_5*dmats0[5][14] + coeff0_6*dmats0[6][14] + coeff0_7*dmats0[7][14] + coeff0_8*dmats0[8][14] + coeff0_9*dmats0[9][14] + coeff0_10*dmats0[10][14] + coeff0_11*dmats0[11][14] + coeff0_12*dmats0[12][14] + coeff0_13*dmats0[13][14] + coeff0_14*dmats0[14][14];
2044
 
          }
2045
 
          if(combinations[deriv_num][j] == 1)
2046
 
          {
2047
 
            new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0] + coeff0_6*dmats1[6][0] + coeff0_7*dmats1[7][0] + coeff0_8*dmats1[8][0] + coeff0_9*dmats1[9][0] + coeff0_10*dmats1[10][0] + coeff0_11*dmats1[11][0] + coeff0_12*dmats1[12][0] + coeff0_13*dmats1[13][0] + coeff0_14*dmats1[14][0];
2048
 
            new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1] + coeff0_6*dmats1[6][1] + coeff0_7*dmats1[7][1] + coeff0_8*dmats1[8][1] + coeff0_9*dmats1[9][1] + coeff0_10*dmats1[10][1] + coeff0_11*dmats1[11][1] + coeff0_12*dmats1[12][1] + coeff0_13*dmats1[13][1] + coeff0_14*dmats1[14][1];
2049
 
            new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2] + coeff0_6*dmats1[6][2] + coeff0_7*dmats1[7][2] + coeff0_8*dmats1[8][2] + coeff0_9*dmats1[9][2] + coeff0_10*dmats1[10][2] + coeff0_11*dmats1[11][2] + coeff0_12*dmats1[12][2] + coeff0_13*dmats1[13][2] + coeff0_14*dmats1[14][2];
2050
 
            new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3] + coeff0_6*dmats1[6][3] + coeff0_7*dmats1[7][3] + coeff0_8*dmats1[8][3] + coeff0_9*dmats1[9][3] + coeff0_10*dmats1[10][3] + coeff0_11*dmats1[11][3] + coeff0_12*dmats1[12][3] + coeff0_13*dmats1[13][3] + coeff0_14*dmats1[14][3];
2051
 
            new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4] + coeff0_6*dmats1[6][4] + coeff0_7*dmats1[7][4] + coeff0_8*dmats1[8][4] + coeff0_9*dmats1[9][4] + coeff0_10*dmats1[10][4] + coeff0_11*dmats1[11][4] + coeff0_12*dmats1[12][4] + coeff0_13*dmats1[13][4] + coeff0_14*dmats1[14][4];
2052
 
            new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5] + coeff0_6*dmats1[6][5] + coeff0_7*dmats1[7][5] + coeff0_8*dmats1[8][5] + coeff0_9*dmats1[9][5] + coeff0_10*dmats1[10][5] + coeff0_11*dmats1[11][5] + coeff0_12*dmats1[12][5] + coeff0_13*dmats1[13][5] + coeff0_14*dmats1[14][5];
2053
 
            new_coeff0_6 = coeff0_0*dmats1[0][6] + coeff0_1*dmats1[1][6] + coeff0_2*dmats1[2][6] + coeff0_3*dmats1[3][6] + coeff0_4*dmats1[4][6] + coeff0_5*dmats1[5][6] + coeff0_6*dmats1[6][6] + coeff0_7*dmats1[7][6] + coeff0_8*dmats1[8][6] + coeff0_9*dmats1[9][6] + coeff0_10*dmats1[10][6] + coeff0_11*dmats1[11][6] + coeff0_12*dmats1[12][6] + coeff0_13*dmats1[13][6] + coeff0_14*dmats1[14][6];
2054
 
            new_coeff0_7 = coeff0_0*dmats1[0][7] + coeff0_1*dmats1[1][7] + coeff0_2*dmats1[2][7] + coeff0_3*dmats1[3][7] + coeff0_4*dmats1[4][7] + coeff0_5*dmats1[5][7] + coeff0_6*dmats1[6][7] + coeff0_7*dmats1[7][7] + coeff0_8*dmats1[8][7] + coeff0_9*dmats1[9][7] + coeff0_10*dmats1[10][7] + coeff0_11*dmats1[11][7] + coeff0_12*dmats1[12][7] + coeff0_13*dmats1[13][7] + coeff0_14*dmats1[14][7];
2055
 
            new_coeff0_8 = coeff0_0*dmats1[0][8] + coeff0_1*dmats1[1][8] + coeff0_2*dmats1[2][8] + coeff0_3*dmats1[3][8] + coeff0_4*dmats1[4][8] + coeff0_5*dmats1[5][8] + coeff0_6*dmats1[6][8] + coeff0_7*dmats1[7][8] + coeff0_8*dmats1[8][8] + coeff0_9*dmats1[9][8] + coeff0_10*dmats1[10][8] + coeff0_11*dmats1[11][8] + coeff0_12*dmats1[12][8] + coeff0_13*dmats1[13][8] + coeff0_14*dmats1[14][8];
2056
 
            new_coeff0_9 = coeff0_0*dmats1[0][9] + coeff0_1*dmats1[1][9] + coeff0_2*dmats1[2][9] + coeff0_3*dmats1[3][9] + coeff0_4*dmats1[4][9] + coeff0_5*dmats1[5][9] + coeff0_6*dmats1[6][9] + coeff0_7*dmats1[7][9] + coeff0_8*dmats1[8][9] + coeff0_9*dmats1[9][9] + coeff0_10*dmats1[10][9] + coeff0_11*dmats1[11][9] + coeff0_12*dmats1[12][9] + coeff0_13*dmats1[13][9] + coeff0_14*dmats1[14][9];
2057
 
            new_coeff0_10 = coeff0_0*dmats1[0][10] + coeff0_1*dmats1[1][10] + coeff0_2*dmats1[2][10] + coeff0_3*dmats1[3][10] + coeff0_4*dmats1[4][10] + coeff0_5*dmats1[5][10] + coeff0_6*dmats1[6][10] + coeff0_7*dmats1[7][10] + coeff0_8*dmats1[8][10] + coeff0_9*dmats1[9][10] + coeff0_10*dmats1[10][10] + coeff0_11*dmats1[11][10] + coeff0_12*dmats1[12][10] + coeff0_13*dmats1[13][10] + coeff0_14*dmats1[14][10];
2058
 
            new_coeff0_11 = coeff0_0*dmats1[0][11] + coeff0_1*dmats1[1][11] + coeff0_2*dmats1[2][11] + coeff0_3*dmats1[3][11] + coeff0_4*dmats1[4][11] + coeff0_5*dmats1[5][11] + coeff0_6*dmats1[6][11] + coeff0_7*dmats1[7][11] + coeff0_8*dmats1[8][11] + coeff0_9*dmats1[9][11] + coeff0_10*dmats1[10][11] + coeff0_11*dmats1[11][11] + coeff0_12*dmats1[12][11] + coeff0_13*dmats1[13][11] + coeff0_14*dmats1[14][11];
2059
 
            new_coeff0_12 = coeff0_0*dmats1[0][12] + coeff0_1*dmats1[1][12] + coeff0_2*dmats1[2][12] + coeff0_3*dmats1[3][12] + coeff0_4*dmats1[4][12] + coeff0_5*dmats1[5][12] + coeff0_6*dmats1[6][12] + coeff0_7*dmats1[7][12] + coeff0_8*dmats1[8][12] + coeff0_9*dmats1[9][12] + coeff0_10*dmats1[10][12] + coeff0_11*dmats1[11][12] + coeff0_12*dmats1[12][12] + coeff0_13*dmats1[13][12] + coeff0_14*dmats1[14][12];
2060
 
            new_coeff0_13 = coeff0_0*dmats1[0][13] + coeff0_1*dmats1[1][13] + coeff0_2*dmats1[2][13] + coeff0_3*dmats1[3][13] + coeff0_4*dmats1[4][13] + coeff0_5*dmats1[5][13] + coeff0_6*dmats1[6][13] + coeff0_7*dmats1[7][13] + coeff0_8*dmats1[8][13] + coeff0_9*dmats1[9][13] + coeff0_10*dmats1[10][13] + coeff0_11*dmats1[11][13] + coeff0_12*dmats1[12][13] + coeff0_13*dmats1[13][13] + coeff0_14*dmats1[14][13];
2061
 
            new_coeff0_14 = coeff0_0*dmats1[0][14] + coeff0_1*dmats1[1][14] + coeff0_2*dmats1[2][14] + coeff0_3*dmats1[3][14] + coeff0_4*dmats1[4][14] + coeff0_5*dmats1[5][14] + coeff0_6*dmats1[6][14] + coeff0_7*dmats1[7][14] + coeff0_8*dmats1[8][14] + coeff0_9*dmats1[9][14] + coeff0_10*dmats1[10][14] + coeff0_11*dmats1[11][14] + coeff0_12*dmats1[12][14] + coeff0_13*dmats1[13][14] + coeff0_14*dmats1[14][14];
2062
 
          }
2063
 
    
2064
 
        }
2065
 
        // Compute derivatives on reference element as dot product of coefficients and basisvalues
2066
 
        derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5 + new_coeff0_6*basisvalue6 + new_coeff0_7*basisvalue7 + new_coeff0_8*basisvalue8 + new_coeff0_9*basisvalue9 + new_coeff0_10*basisvalue10 + new_coeff0_11*basisvalue11 + new_coeff0_12*basisvalue12 + new_coeff0_13*basisvalue13 + new_coeff0_14*basisvalue14;
2067
 
      }
2068
 
    
2069
 
      // Transform derivatives back to physical element
2070
 
      for (unsigned int row = 0; row < num_derivatives; row++)
2071
 
      {
2072
 
        for (unsigned int col = 0; col < num_derivatives; col++)
2073
 
        {
2074
 
          values[num_derivatives + row] += transform[row][col]*derivatives[col];
2075
 
        }
2076
 
      }
2077
 
      // Delete pointer to array of derivatives on FIAT element
2078
 
      delete [] derivatives;
2079
 
    
2080
 
      // Delete pointer to array of combinations of derivatives and transform
2081
 
      for (unsigned int row = 0; row < num_derivatives; row++)
2082
 
      {
2083
 
        delete [] combinations[row];
2084
 
        delete [] transform[row];
2085
 
      }
2086
 
    
2087
 
      delete [] combinations;
2088
 
      delete [] transform;
2089
 
    }
2090
 
    
2091
 
  }
2092
 
 
2093
 
  /// Evaluate order n derivatives of all basis functions at given point in cell
2094
 
  virtual void evaluate_basis_derivatives_all(unsigned int n,
2095
 
                                              double* values,
2096
 
                                              const double* coordinates,
2097
 
                                              const ufc::cell& c) const
2098
 
  {
2099
 
    throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
2100
 
  }
2101
 
 
2102
 
  /// Evaluate linear functional for dof i on the function f
2103
 
  virtual double evaluate_dof(unsigned int i,
2104
 
                              const ufc::function& f,
2105
 
                              const ufc::cell& c) const
2106
 
  {
2107
 
    // The reference points, direction and weights:
2108
 
    static const double X[24][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.75, 0.25}}, {{0.5, 0.5}}, {{0.25, 0.75}}, {{0, 0.25}}, {{0, 0.5}}, {{0, 0.75}}, {{0.25, 0}}, {{0.5, 0}}, {{0.75, 0}}, {{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.75, 0.25}}, {{0.5, 0.5}}, {{0.25, 0.75}}, {{0, 0.25}}, {{0, 0.5}}, {{0, 0.75}}, {{0.25, 0}}, {{0.5, 0}}, {{0.75, 0}}};
2109
 
    static const double W[24][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
2110
 
    static const double D[24][1][2] = {{{1, 0}}, {{1, 0}}, {{1, 0}}, {{1, 0}}, {{1, 0}}, {{1, 0}}, {{1, 0}}, {{1, 0}}, {{1, 0}}, {{1, 0}}, {{1, 0}}, {{1, 0}}, {{0, 1}}, {{0, 1}}, {{0, 1}}, {{0, 1}}, {{0, 1}}, {{0, 1}}, {{0, 1}}, {{0, 1}}, {{0, 1}}, {{0, 1}}, {{0, 1}}, {{0, 1}}};
2111
 
    
2112
 
    const double * const * x = c.coordinates;
2113
 
    double result = 0.0;
2114
 
    // Iterate over the points:
2115
 
    // Evaluate basis functions for affine mapping
2116
 
    const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
2117
 
    const double w1 = X[i][0][0];
2118
 
    const double w2 = X[i][0][1];
2119
 
    
2120
 
    // Compute affine mapping y = F(X)
2121
 
    double y[2];
2122
 
    y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
2123
 
    y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
2124
 
    
2125
 
    // Evaluate function at physical points
2126
 
    double values[2];
2127
 
    f.evaluate(values, y, c);
2128
 
    
2129
 
    // Map function values using appropriate mapping
2130
 
    // Affine map: Do nothing
2131
 
    
2132
 
    // Note that we do not map the weights (yet).
2133
 
    
2134
 
    // Take directional components
2135
 
    for(int k = 0; k < 2; k++)
2136
 
      result += values[k]*D[i][0][k];
2137
 
    // Multiply by weights
2138
 
    result *= W[i][0];
2139
 
    
2140
 
    return result;
2141
 
  }
2142
 
 
2143
 
  /// Evaluate linear functionals for all dofs on the function f
2144
 
  virtual void evaluate_dofs(double* values,
2145
 
                             const ufc::function& f,
2146
 
                             const ufc::cell& c) const
2147
 
  {
2148
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2149
 
  }
2150
 
 
2151
 
  /// Interpolate vertex values from dof values
2152
 
  virtual void interpolate_vertex_values(double* vertex_values,
2153
 
                                         const double* dof_values,
2154
 
                                         const ufc::cell& c) const
2155
 
  {
2156
 
    // Evaluate at vertices and use affine mapping
2157
 
    vertex_values[0] = dof_values[0];
2158
 
    vertex_values[2] = dof_values[1];
2159
 
    vertex_values[4] = dof_values[2];
2160
 
    // Evaluate at vertices and use affine mapping
2161
 
    vertex_values[1] = dof_values[12];
2162
 
    vertex_values[3] = dof_values[13];
2163
 
    vertex_values[5] = dof_values[14];
2164
 
  }
2165
 
 
2166
 
  /// Return the number of sub elements (for a mixed element)
2167
 
  virtual unsigned int num_sub_elements() const
2168
 
  {
2169
 
    return 2;
2170
 
  }
2171
 
 
2172
 
  /// Create a new finite element for sub element i (for a mixed element)
2173
 
  virtual ufc::finite_element* create_sub_element(unsigned int i) const
2174
 
  {
2175
 
    switch ( i )
2176
 
    {
2177
 
    case 0:
2178
 
      return new elementrestriction_0_finite_element_0_0();
2179
 
      break;
2180
 
    case 1:
2181
 
      return new elementrestriction_0_finite_element_0_1();
2182
 
      break;
2183
 
    }
2184
 
    return 0;
2185
 
  }
2186
 
 
2187
 
};
2188
 
 
2189
 
/// This class defines the interface for a local-to-global mapping of
2190
 
/// degrees of freedom (dofs).
2191
 
 
2192
 
class elementrestriction_0_dof_map_0_0: public ufc::dof_map
2193
 
{
2194
 
private:
2195
 
 
2196
 
  unsigned int __global_dimension;
2197
 
 
2198
 
public:
2199
 
 
2200
 
  /// Constructor
2201
 
  elementrestriction_0_dof_map_0_0() : ufc::dof_map()
2202
 
  {
2203
 
    __global_dimension = 0;
2204
 
  }
2205
 
 
2206
 
  /// Destructor
2207
 
  virtual ~elementrestriction_0_dof_map_0_0()
2208
 
  {
2209
 
    // Do nothing
2210
 
  }
2211
 
 
2212
 
  /// Return a string identifying the dof map
2213
 
  virtual const char* signature() const
2214
 
  {
2215
 
    return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 4)";
2216
 
  }
2217
 
 
2218
 
  /// Return true iff mesh entities of topological dimension d are needed
2219
 
  virtual bool needs_mesh_entities(unsigned int d) const
2220
 
  {
2221
 
    switch ( d )
2222
 
    {
2223
 
    case 0:
2224
 
      return true;
2225
 
      break;
2226
 
    case 1:
2227
 
      return true;
2228
 
      break;
2229
 
    case 2:
2230
 
      return false;
2231
 
      break;
2232
 
    }
2233
 
    return false;
2234
 
  }
2235
 
 
2236
 
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
2237
 
  virtual bool init_mesh(const ufc::mesh& m)
2238
 
  {
2239
 
    __global_dimension = m.num_entities[0] + 3*m.num_entities[1];
2240
 
    return false;
2241
 
  }
2242
 
 
2243
 
  /// Initialize dof map for given cell
2244
 
  virtual void init_cell(const ufc::mesh& m,
2245
 
                         const ufc::cell& c)
2246
 
  {
2247
 
    // Do nothing
2248
 
  }
2249
 
 
2250
 
  /// Finish initialization of dof map for cells
2251
 
  virtual void init_cell_finalize()
2252
 
  {
2253
 
    // Do nothing
2254
 
  }
2255
 
 
2256
 
  /// Return the dimension of the global finite element function space
2257
 
  virtual unsigned int global_dimension() const
2258
 
  {
2259
 
    return __global_dimension;
2260
 
  }
2261
 
 
2262
 
  /// Return the dimension of the local finite element function space for a cell
2263
 
  virtual unsigned int local_dimension(const ufc::cell& c) const
2264
 
  {
2265
 
    return 12;
2266
 
  }
2267
 
 
2268
 
  /// Return the maximum dimension of the local finite element function space
2269
 
  virtual unsigned int max_local_dimension() const
2270
 
  {
2271
 
    return 12;
2272
 
  }
2273
 
 
2274
 
  // Return the geometric dimension of the coordinates this dof map provides
2275
 
  virtual unsigned int geometric_dimension() const
2276
 
  {
2277
 
    return 2;
2278
 
  }
2279
 
 
2280
 
  /// Return the number of dofs on each cell facet
2281
 
  virtual unsigned int num_facet_dofs() const
2282
 
  {
2283
 
    return 5;
2284
 
  }
2285
 
 
2286
 
  /// Return the number of dofs associated with each cell entity of dimension d
2287
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
2288
 
  {
2289
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2290
 
  }
2291
 
 
2292
 
  /// Tabulate the local-to-global mapping of dofs on a cell
2293
 
  virtual void tabulate_dofs(unsigned int* dofs,
2294
 
                             const ufc::mesh& m,
2295
 
                             const ufc::cell& c) const
2296
 
  {
2297
 
    dofs[0] = c.entity_indices[0][0];
2298
 
    dofs[1] = c.entity_indices[0][1];
2299
 
    dofs[2] = c.entity_indices[0][2];
2300
 
    unsigned int offset = m.num_entities[0];
2301
 
    dofs[3] = offset + 3*c.entity_indices[1][0];
2302
 
    dofs[4] = offset + 3*c.entity_indices[1][0] + 1;
2303
 
    dofs[5] = offset + 3*c.entity_indices[1][0] + 2;
2304
 
    dofs[6] = offset + 3*c.entity_indices[1][1];
2305
 
    dofs[7] = offset + 3*c.entity_indices[1][1] + 1;
2306
 
    dofs[8] = offset + 3*c.entity_indices[1][1] + 2;
2307
 
    dofs[9] = offset + 3*c.entity_indices[1][2];
2308
 
    dofs[10] = offset + 3*c.entity_indices[1][2] + 1;
2309
 
    dofs[11] = offset + 3*c.entity_indices[1][2] + 2;
2310
 
  }
2311
 
 
2312
 
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
2313
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
2314
 
                                   unsigned int facet) const
2315
 
  {
2316
 
    switch ( facet )
2317
 
    {
2318
 
    case 0:
2319
 
      dofs[0] = 1;
2320
 
      dofs[1] = 2;
2321
 
      dofs[2] = 3;
2322
 
      dofs[3] = 4;
2323
 
      dofs[4] = 5;
2324
 
      break;
2325
 
    case 1:
2326
 
      dofs[0] = 0;
2327
 
      dofs[1] = 2;
2328
 
      dofs[2] = 6;
2329
 
      dofs[3] = 7;
2330
 
      dofs[4] = 8;
2331
 
      break;
2332
 
    case 2:
2333
 
      dofs[0] = 0;
2334
 
      dofs[1] = 1;
2335
 
      dofs[2] = 9;
2336
 
      dofs[3] = 10;
2337
 
      dofs[4] = 11;
2338
 
      break;
2339
 
    }
2340
 
  }
2341
 
 
2342
 
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
2343
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
2344
 
                                    unsigned int d, unsigned int i) const
2345
 
  {
2346
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2347
 
  }
2348
 
 
2349
 
  /// Tabulate the coordinates of all dofs on a cell
2350
 
  virtual void tabulate_coordinates(double** coordinates,
2351
 
                                    const ufc::cell& c) const
2352
 
  {
2353
 
    const double * const * x = c.coordinates;
2354
 
    coordinates[0][0] = x[0][0];
2355
 
    coordinates[0][1] = x[0][1];
2356
 
    coordinates[1][0] = x[1][0];
2357
 
    coordinates[1][1] = x[1][1];
2358
 
    coordinates[2][0] = x[2][0];
2359
 
    coordinates[2][1] = x[2][1];
2360
 
    coordinates[3][0] = 0.75*x[1][0] + 0.25*x[2][0];
2361
 
    coordinates[3][1] = 0.75*x[1][1] + 0.25*x[2][1];
2362
 
    coordinates[4][0] = 0.5*x[1][0] + 0.5*x[2][0];
2363
 
    coordinates[4][1] = 0.5*x[1][1] + 0.5*x[2][1];
2364
 
    coordinates[5][0] = 0.25*x[1][0] + 0.75*x[2][0];
2365
 
    coordinates[5][1] = 0.25*x[1][1] + 0.75*x[2][1];
2366
 
    coordinates[6][0] = 0.75*x[0][0] + 0.25*x[2][0];
2367
 
    coordinates[6][1] = 0.75*x[0][1] + 0.25*x[2][1];
2368
 
    coordinates[7][0] = 0.5*x[0][0] + 0.5*x[2][0];
2369
 
    coordinates[7][1] = 0.5*x[0][1] + 0.5*x[2][1];
2370
 
    coordinates[8][0] = 0.25*x[0][0] + 0.75*x[2][0];
2371
 
    coordinates[8][1] = 0.25*x[0][1] + 0.75*x[2][1];
2372
 
    coordinates[9][0] = 0.75*x[0][0] + 0.25*x[1][0];
2373
 
    coordinates[9][1] = 0.75*x[0][1] + 0.25*x[1][1];
2374
 
    coordinates[10][0] = 0.5*x[0][0] + 0.5*x[1][0];
2375
 
    coordinates[10][1] = 0.5*x[0][1] + 0.5*x[1][1];
2376
 
    coordinates[11][0] = 0.25*x[0][0] + 0.75*x[1][0];
2377
 
    coordinates[11][1] = 0.25*x[0][1] + 0.75*x[1][1];
2378
 
  }
2379
 
 
2380
 
  /// Return the number of sub dof maps (for a mixed element)
2381
 
  virtual unsigned int num_sub_dof_maps() const
2382
 
  {
2383
 
    return 1;
2384
 
  }
2385
 
 
2386
 
  /// Create a new dof_map for sub dof map i (for a mixed element)
2387
 
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2388
 
  {
2389
 
    return new elementrestriction_0_dof_map_0_0();
2390
 
  }
2391
 
 
2392
 
};
2393
 
 
2394
 
/// This class defines the interface for a local-to-global mapping of
2395
 
/// degrees of freedom (dofs).
2396
 
 
2397
 
class elementrestriction_0_dof_map_0_1: public ufc::dof_map
2398
 
{
2399
 
private:
2400
 
 
2401
 
  unsigned int __global_dimension;
2402
 
 
2403
 
public:
2404
 
 
2405
 
  /// Constructor
2406
 
  elementrestriction_0_dof_map_0_1() : ufc::dof_map()
2407
 
  {
2408
 
    __global_dimension = 0;
2409
 
  }
2410
 
 
2411
 
  /// Destructor
2412
 
  virtual ~elementrestriction_0_dof_map_0_1()
2413
 
  {
2414
 
    // Do nothing
2415
 
  }
2416
 
 
2417
 
  /// Return a string identifying the dof map
2418
 
  virtual const char* signature() const
2419
 
  {
2420
 
    return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 4)";
2421
 
  }
2422
 
 
2423
 
  /// Return true iff mesh entities of topological dimension d are needed
2424
 
  virtual bool needs_mesh_entities(unsigned int d) const
2425
 
  {
2426
 
    switch ( d )
2427
 
    {
2428
 
    case 0:
2429
 
      return true;
2430
 
      break;
2431
 
    case 1:
2432
 
      return true;
2433
 
      break;
2434
 
    case 2:
2435
 
      return false;
2436
 
      break;
2437
 
    }
2438
 
    return false;
2439
 
  }
2440
 
 
2441
 
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
2442
 
  virtual bool init_mesh(const ufc::mesh& m)
2443
 
  {
2444
 
    __global_dimension = m.num_entities[0] + 3*m.num_entities[1];
2445
 
    return false;
2446
 
  }
2447
 
 
2448
 
  /// Initialize dof map for given cell
2449
 
  virtual void init_cell(const ufc::mesh& m,
2450
 
                         const ufc::cell& c)
2451
 
  {
2452
 
    // Do nothing
2453
 
  }
2454
 
 
2455
 
  /// Finish initialization of dof map for cells
2456
 
  virtual void init_cell_finalize()
2457
 
  {
2458
 
    // Do nothing
2459
 
  }
2460
 
 
2461
 
  /// Return the dimension of the global finite element function space
2462
 
  virtual unsigned int global_dimension() const
2463
 
  {
2464
 
    return __global_dimension;
2465
 
  }
2466
 
 
2467
 
  /// Return the dimension of the local finite element function space for a cell
2468
 
  virtual unsigned int local_dimension(const ufc::cell& c) const
2469
 
  {
2470
 
    return 12;
2471
 
  }
2472
 
 
2473
 
  /// Return the maximum dimension of the local finite element function space
2474
 
  virtual unsigned int max_local_dimension() const
2475
 
  {
2476
 
    return 12;
2477
 
  }
2478
 
 
2479
 
  // Return the geometric dimension of the coordinates this dof map provides
2480
 
  virtual unsigned int geometric_dimension() const
2481
 
  {
2482
 
    return 2;
2483
 
  }
2484
 
 
2485
 
  /// Return the number of dofs on each cell facet
2486
 
  virtual unsigned int num_facet_dofs() const
2487
 
  {
2488
 
    return 5;
2489
 
  }
2490
 
 
2491
 
  /// Return the number of dofs associated with each cell entity of dimension d
2492
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
2493
 
  {
2494
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2495
 
  }
2496
 
 
2497
 
  /// Tabulate the local-to-global mapping of dofs on a cell
2498
 
  virtual void tabulate_dofs(unsigned int* dofs,
2499
 
                             const ufc::mesh& m,
2500
 
                             const ufc::cell& c) const
2501
 
  {
2502
 
    dofs[0] = c.entity_indices[0][0];
2503
 
    dofs[1] = c.entity_indices[0][1];
2504
 
    dofs[2] = c.entity_indices[0][2];
2505
 
    unsigned int offset = m.num_entities[0];
2506
 
    dofs[3] = offset + 3*c.entity_indices[1][0];
2507
 
    dofs[4] = offset + 3*c.entity_indices[1][0] + 1;
2508
 
    dofs[5] = offset + 3*c.entity_indices[1][0] + 2;
2509
 
    dofs[6] = offset + 3*c.entity_indices[1][1];
2510
 
    dofs[7] = offset + 3*c.entity_indices[1][1] + 1;
2511
 
    dofs[8] = offset + 3*c.entity_indices[1][1] + 2;
2512
 
    dofs[9] = offset + 3*c.entity_indices[1][2];
2513
 
    dofs[10] = offset + 3*c.entity_indices[1][2] + 1;
2514
 
    dofs[11] = offset + 3*c.entity_indices[1][2] + 2;
2515
 
  }
2516
 
 
2517
 
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
2518
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
2519
 
                                   unsigned int facet) const
2520
 
  {
2521
 
    switch ( facet )
2522
 
    {
2523
 
    case 0:
2524
 
      dofs[0] = 1;
2525
 
      dofs[1] = 2;
2526
 
      dofs[2] = 3;
2527
 
      dofs[3] = 4;
2528
 
      dofs[4] = 5;
2529
 
      break;
2530
 
    case 1:
2531
 
      dofs[0] = 0;
2532
 
      dofs[1] = 2;
2533
 
      dofs[2] = 6;
2534
 
      dofs[3] = 7;
2535
 
      dofs[4] = 8;
2536
 
      break;
2537
 
    case 2:
2538
 
      dofs[0] = 0;
2539
 
      dofs[1] = 1;
2540
 
      dofs[2] = 9;
2541
 
      dofs[3] = 10;
2542
 
      dofs[4] = 11;
2543
 
      break;
2544
 
    }
2545
 
  }
2546
 
 
2547
 
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
2548
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
2549
 
                                    unsigned int d, unsigned int i) const
2550
 
  {
2551
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2552
 
  }
2553
 
 
2554
 
  /// Tabulate the coordinates of all dofs on a cell
2555
 
  virtual void tabulate_coordinates(double** coordinates,
2556
 
                                    const ufc::cell& c) const
2557
 
  {
2558
 
    const double * const * x = c.coordinates;
2559
 
    coordinates[0][0] = x[0][0];
2560
 
    coordinates[0][1] = x[0][1];
2561
 
    coordinates[1][0] = x[1][0];
2562
 
    coordinates[1][1] = x[1][1];
2563
 
    coordinates[2][0] = x[2][0];
2564
 
    coordinates[2][1] = x[2][1];
2565
 
    coordinates[3][0] = 0.75*x[1][0] + 0.25*x[2][0];
2566
 
    coordinates[3][1] = 0.75*x[1][1] + 0.25*x[2][1];
2567
 
    coordinates[4][0] = 0.5*x[1][0] + 0.5*x[2][0];
2568
 
    coordinates[4][1] = 0.5*x[1][1] + 0.5*x[2][1];
2569
 
    coordinates[5][0] = 0.25*x[1][0] + 0.75*x[2][0];
2570
 
    coordinates[5][1] = 0.25*x[1][1] + 0.75*x[2][1];
2571
 
    coordinates[6][0] = 0.75*x[0][0] + 0.25*x[2][0];
2572
 
    coordinates[6][1] = 0.75*x[0][1] + 0.25*x[2][1];
2573
 
    coordinates[7][0] = 0.5*x[0][0] + 0.5*x[2][0];
2574
 
    coordinates[7][1] = 0.5*x[0][1] + 0.5*x[2][1];
2575
 
    coordinates[8][0] = 0.25*x[0][0] + 0.75*x[2][0];
2576
 
    coordinates[8][1] = 0.25*x[0][1] + 0.75*x[2][1];
2577
 
    coordinates[9][0] = 0.75*x[0][0] + 0.25*x[1][0];
2578
 
    coordinates[9][1] = 0.75*x[0][1] + 0.25*x[1][1];
2579
 
    coordinates[10][0] = 0.5*x[0][0] + 0.5*x[1][0];
2580
 
    coordinates[10][1] = 0.5*x[0][1] + 0.5*x[1][1];
2581
 
    coordinates[11][0] = 0.25*x[0][0] + 0.75*x[1][0];
2582
 
    coordinates[11][1] = 0.25*x[0][1] + 0.75*x[1][1];
2583
 
  }
2584
 
 
2585
 
  /// Return the number of sub dof maps (for a mixed element)
2586
 
  virtual unsigned int num_sub_dof_maps() const
2587
 
  {
2588
 
    return 1;
2589
 
  }
2590
 
 
2591
 
  /// Create a new dof_map for sub dof map i (for a mixed element)
2592
 
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2593
 
  {
2594
 
    return new elementrestriction_0_dof_map_0_1();
2595
 
  }
2596
 
 
2597
 
};
2598
 
 
2599
 
/// This class defines the interface for a local-to-global mapping of
2600
 
/// degrees of freedom (dofs).
2601
 
 
2602
 
class elementrestriction_0_dof_map_0: public ufc::dof_map
2603
 
{
2604
 
private:
2605
 
 
2606
 
  unsigned int __global_dimension;
2607
 
 
2608
 
public:
2609
 
 
2610
 
  /// Constructor
2611
 
  elementrestriction_0_dof_map_0() : ufc::dof_map()
2612
 
  {
2613
 
    __global_dimension = 0;
2614
 
  }
2615
 
 
2616
 
  /// Destructor
2617
 
  virtual ~elementrestriction_0_dof_map_0()
2618
 
  {
2619
 
    // Do nothing
2620
 
  }
2621
 
 
2622
 
  /// Return a string identifying the dof map
2623
 
  virtual const char* signature() const
2624
 
  {
2625
 
    return "FFC dof map for MixedElement(*[FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 4), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 4), Cell('interval', 1, Space(1)))], **{'value_shape': (2,) })";
2626
 
  }
2627
 
 
2628
 
  /// Return true iff mesh entities of topological dimension d are needed
2629
 
  virtual bool needs_mesh_entities(unsigned int d) const
2630
 
  {
2631
 
    switch ( d )
2632
 
    {
2633
 
    case 0:
2634
 
      return true;
2635
 
      break;
2636
 
    case 1:
2637
 
      return true;
2638
 
      break;
2639
 
    case 2:
2640
 
      return false;
2641
 
      break;
2642
 
    }
2643
 
    return false;
2644
 
  }
2645
 
 
2646
 
  /// Initialize dof map for mesh (return true iff init_cell() is needed)
2647
 
  virtual bool init_mesh(const ufc::mesh& m)
2648
 
  {
2649
 
    __global_dimension = 2*m.num_entities[0] + 6*m.num_entities[1];
2650
 
    return false;
2651
 
  }
2652
 
 
2653
 
  /// Initialize dof map for given cell
2654
 
  virtual void init_cell(const ufc::mesh& m,
2655
 
                         const ufc::cell& c)
2656
 
  {
2657
 
    // Do nothing
2658
 
  }
2659
 
 
2660
 
  /// Finish initialization of dof map for cells
2661
 
  virtual void init_cell_finalize()
2662
 
  {
2663
 
    // Do nothing
2664
 
  }
2665
 
 
2666
 
  /// Return the dimension of the global finite element function space
2667
 
  virtual unsigned int global_dimension() const
2668
 
  {
2669
 
    return __global_dimension;
2670
 
  }
2671
 
 
2672
 
  /// Return the dimension of the local finite element function space for a cell
2673
 
  virtual unsigned int local_dimension(const ufc::cell& c) const
2674
 
  {
2675
 
    return 24;
2676
 
  }
2677
 
 
2678
 
  /// Return the maximum dimension of the local finite element function space
2679
 
  virtual unsigned int max_local_dimension() const
2680
 
  {
2681
 
    return 24;
2682
 
  }
2683
 
 
2684
 
  // Return the geometric dimension of the coordinates this dof map provides
2685
 
  virtual unsigned int geometric_dimension() const
2686
 
  {
2687
 
    return 2;
2688
 
  }
2689
 
 
2690
 
  /// Return the number of dofs on each cell facet
2691
 
  virtual unsigned int num_facet_dofs() const
2692
 
  {
2693
 
    return 10;
2694
 
  }
2695
 
 
2696
 
  /// Return the number of dofs associated with each cell entity of dimension d
2697
 
  virtual unsigned int num_entity_dofs(unsigned int d) const
2698
 
  {
2699
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2700
 
  }
2701
 
 
2702
 
  /// Tabulate the local-to-global mapping of dofs on a cell
2703
 
  virtual void tabulate_dofs(unsigned int* dofs,
2704
 
                             const ufc::mesh& m,
2705
 
                             const ufc::cell& c) const
2706
 
  {
2707
 
    dofs[0] = c.entity_indices[0][0];
2708
 
    dofs[1] = c.entity_indices[0][1];
2709
 
    dofs[2] = c.entity_indices[0][2];
2710
 
    unsigned int offset = m.num_entities[0];
2711
 
    dofs[3] = offset + 3*c.entity_indices[1][0];
2712
 
    dofs[4] = offset + 3*c.entity_indices[1][0] + 1;
2713
 
    dofs[5] = offset + 3*c.entity_indices[1][0] + 2;
2714
 
    dofs[6] = offset + 3*c.entity_indices[1][1];
2715
 
    dofs[7] = offset + 3*c.entity_indices[1][1] + 1;
2716
 
    dofs[8] = offset + 3*c.entity_indices[1][1] + 2;
2717
 
    dofs[9] = offset + 3*c.entity_indices[1][2];
2718
 
    dofs[10] = offset + 3*c.entity_indices[1][2] + 1;
2719
 
    dofs[11] = offset + 3*c.entity_indices[1][2] + 2;
2720
 
    offset = offset + 3*m.num_entities[1];
2721
 
    dofs[12] = offset + c.entity_indices[0][0];
2722
 
    dofs[13] = offset + c.entity_indices[0][1];
2723
 
    dofs[14] = offset + c.entity_indices[0][2];
2724
 
    offset = offset + m.num_entities[0];
2725
 
    dofs[15] = offset + 3*c.entity_indices[1][0];
2726
 
    dofs[16] = offset + 3*c.entity_indices[1][0] + 1;
2727
 
    dofs[17] = offset + 3*c.entity_indices[1][0] + 2;
2728
 
    dofs[18] = offset + 3*c.entity_indices[1][1];
2729
 
    dofs[19] = offset + 3*c.entity_indices[1][1] + 1;
2730
 
    dofs[20] = offset + 3*c.entity_indices[1][1] + 2;
2731
 
    dofs[21] = offset + 3*c.entity_indices[1][2];
2732
 
    dofs[22] = offset + 3*c.entity_indices[1][2] + 1;
2733
 
    dofs[23] = offset + 3*c.entity_indices[1][2] + 2;
2734
 
  }
2735
 
 
2736
 
  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
2737
 
  virtual void tabulate_facet_dofs(unsigned int* dofs,
2738
 
                                   unsigned int facet) const
2739
 
  {
2740
 
    switch ( facet )
2741
 
    {
2742
 
    case 0:
2743
 
      dofs[0] = 1;
2744
 
      dofs[1] = 2;
2745
 
      dofs[2] = 3;
2746
 
      dofs[3] = 4;
2747
 
      dofs[4] = 5;
2748
 
      dofs[5] = 13;
2749
 
      dofs[6] = 14;
2750
 
      dofs[7] = 15;
2751
 
      dofs[8] = 16;
2752
 
      dofs[9] = 17;
2753
 
      break;
2754
 
    case 1:
2755
 
      dofs[0] = 0;
2756
 
      dofs[1] = 2;
2757
 
      dofs[2] = 6;
2758
 
      dofs[3] = 7;
2759
 
      dofs[4] = 8;
2760
 
      dofs[5] = 12;
2761
 
      dofs[6] = 14;
2762
 
      dofs[7] = 18;
2763
 
      dofs[8] = 19;
2764
 
      dofs[9] = 20;
2765
 
      break;
2766
 
    case 2:
2767
 
      dofs[0] = 0;
2768
 
      dofs[1] = 1;
2769
 
      dofs[2] = 9;
2770
 
      dofs[3] = 10;
2771
 
      dofs[4] = 11;
2772
 
      dofs[5] = 12;
2773
 
      dofs[6] = 13;
2774
 
      dofs[7] = 21;
2775
 
      dofs[8] = 22;
2776
 
      dofs[9] = 23;
2777
 
      break;
2778
 
    }
2779
 
  }
2780
 
 
2781
 
  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
2782
 
  virtual void tabulate_entity_dofs(unsigned int* dofs,
2783
 
                                    unsigned int d, unsigned int i) const
2784
 
  {
2785
 
    throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
2786
 
  }
2787
 
 
2788
 
  /// Tabulate the coordinates of all dofs on a cell
2789
 
  virtual void tabulate_coordinates(double** coordinates,
2790
 
                                    const ufc::cell& c) const
2791
 
  {
2792
 
    const double * const * x = c.coordinates;
2793
 
    coordinates[0][0] = x[0][0];
2794
 
    coordinates[0][1] = x[0][1];
2795
 
    coordinates[1][0] = x[1][0];
2796
 
    coordinates[1][1] = x[1][1];
2797
 
    coordinates[2][0] = x[2][0];
2798
 
    coordinates[2][1] = x[2][1];
2799
 
    coordinates[3][0] = 0.75*x[1][0] + 0.25*x[2][0];
2800
 
    coordinates[3][1] = 0.75*x[1][1] + 0.25*x[2][1];
2801
 
    coordinates[4][0] = 0.5*x[1][0] + 0.5*x[2][0];
2802
 
    coordinates[4][1] = 0.5*x[1][1] + 0.5*x[2][1];
2803
 
    coordinates[5][0] = 0.25*x[1][0] + 0.75*x[2][0];
2804
 
    coordinates[5][1] = 0.25*x[1][1] + 0.75*x[2][1];
2805
 
    coordinates[6][0] = 0.75*x[0][0] + 0.25*x[2][0];
2806
 
    coordinates[6][1] = 0.75*x[0][1] + 0.25*x[2][1];
2807
 
    coordinates[7][0] = 0.5*x[0][0] + 0.5*x[2][0];
2808
 
    coordinates[7][1] = 0.5*x[0][1] + 0.5*x[2][1];
2809
 
    coordinates[8][0] = 0.25*x[0][0] + 0.75*x[2][0];
2810
 
    coordinates[8][1] = 0.25*x[0][1] + 0.75*x[2][1];
2811
 
    coordinates[9][0] = 0.75*x[0][0] + 0.25*x[1][0];
2812
 
    coordinates[9][1] = 0.75*x[0][1] + 0.25*x[1][1];
2813
 
    coordinates[10][0] = 0.5*x[0][0] + 0.5*x[1][0];
2814
 
    coordinates[10][1] = 0.5*x[0][1] + 0.5*x[1][1];
2815
 
    coordinates[11][0] = 0.25*x[0][0] + 0.75*x[1][0];
2816
 
    coordinates[11][1] = 0.25*x[0][1] + 0.75*x[1][1];
2817
 
    coordinates[12][0] = x[0][0];
2818
 
    coordinates[12][1] = x[0][1];
2819
 
    coordinates[13][0] = x[1][0];
2820
 
    coordinates[13][1] = x[1][1];
2821
 
    coordinates[14][0] = x[2][0];
2822
 
    coordinates[14][1] = x[2][1];
2823
 
    coordinates[15][0] = 0.75*x[1][0] + 0.25*x[2][0];
2824
 
    coordinates[15][1] = 0.75*x[1][1] + 0.25*x[2][1];
2825
 
    coordinates[16][0] = 0.5*x[1][0] + 0.5*x[2][0];
2826
 
    coordinates[16][1] = 0.5*x[1][1] + 0.5*x[2][1];
2827
 
    coordinates[17][0] = 0.25*x[1][0] + 0.75*x[2][0];
2828
 
    coordinates[17][1] = 0.25*x[1][1] + 0.75*x[2][1];
2829
 
    coordinates[18][0] = 0.75*x[0][0] + 0.25*x[2][0];
2830
 
    coordinates[18][1] = 0.75*x[0][1] + 0.25*x[2][1];
2831
 
    coordinates[19][0] = 0.5*x[0][0] + 0.5*x[2][0];
2832
 
    coordinates[19][1] = 0.5*x[0][1] + 0.5*x[2][1];
2833
 
    coordinates[20][0] = 0.25*x[0][0] + 0.75*x[2][0];
2834
 
    coordinates[20][1] = 0.25*x[0][1] + 0.75*x[2][1];
2835
 
    coordinates[21][0] = 0.75*x[0][0] + 0.25*x[1][0];
2836
 
    coordinates[21][1] = 0.75*x[0][1] + 0.25*x[1][1];
2837
 
    coordinates[22][0] = 0.5*x[0][0] + 0.5*x[1][0];
2838
 
    coordinates[22][1] = 0.5*x[0][1] + 0.5*x[1][1];
2839
 
    coordinates[23][0] = 0.25*x[0][0] + 0.75*x[1][0];
2840
 
    coordinates[23][1] = 0.25*x[0][1] + 0.75*x[1][1];
2841
 
  }
2842
 
 
2843
 
  /// Return the number of sub dof maps (for a mixed element)
2844
 
  virtual unsigned int num_sub_dof_maps() const
2845
 
  {
2846
 
    return 2;
2847
 
  }
2848
 
 
2849
 
  /// Create a new dof_map for sub dof map i (for a mixed element)
2850
 
  virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2851
 
  {
2852
 
    switch ( i )
2853
 
    {
2854
 
    case 0:
2855
 
      return new elementrestriction_0_dof_map_0_0();
2856
 
      break;
2857
 
    case 1:
2858
 
      return new elementrestriction_0_dof_map_0_1();
2859
 
      break;
2860
 
    }
2861
 
    return 0;
2862
 
  }
2863
 
 
2864
 
};
2865
 
 
2866
 
/// This class defines the interface for the assembly of the global
2867
 
/// tensor corresponding to a form with r + n arguments, that is, a
2868
 
/// mapping
2869
 
///
2870
 
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
2871
 
///
2872
 
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
2873
 
/// global tensor A is defined by
2874
 
///
2875
 
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
2876
 
///
2877
 
/// where each argument Vj represents the application to the
2878
 
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
2879
 
/// fixed functions (coefficients).
2880
 
 
2881
 
class elementrestriction_form_0: public ufc::form
2882
 
{
2883
 
public:
2884
 
 
2885
 
  /// Constructor
2886
 
  elementrestriction_form_0() : ufc::form()
2887
 
  {
2888
 
    // Do nothing
2889
 
  }
2890
 
 
2891
 
  /// Destructor
2892
 
  virtual ~elementrestriction_form_0()
2893
 
  {
2894
 
    // Do nothing
2895
 
  }
2896
 
 
2897
 
  /// Return a string identifying the form
2898
 
  virtual const char* signature() const
2899
 
  {
2900
 
    return "None";
2901
 
  }
2902
 
 
2903
 
  /// Return the rank of the global tensor (r)
2904
 
  virtual unsigned int rank() const
2905
 
  {
2906
 
    return -1;
2907
 
  }
2908
 
 
2909
 
  /// Return the number of coefficients (n)
2910
 
  virtual unsigned int num_coefficients() const
2911
 
  {
2912
 
    return 0;
2913
 
  }
2914
 
 
2915
 
  /// Return the number of cell integrals
2916
 
  virtual unsigned int num_cell_integrals() const
2917
 
  {
2918
 
    return 0;
2919
 
  }
2920
 
 
2921
 
  /// Return the number of exterior facet integrals
2922
 
  virtual unsigned int num_exterior_facet_integrals() const
2923
 
  {
2924
 
    return 0;
2925
 
  }
2926
 
 
2927
 
  /// Return the number of interior facet integrals
2928
 
  virtual unsigned int num_interior_facet_integrals() const
2929
 
  {
2930
 
    return 0;
2931
 
  }
2932
 
 
2933
 
  /// Create a new finite element for argument function i
2934
 
  virtual ufc::finite_element* create_finite_element(unsigned int i) const
2935
 
  {
2936
 
    return 0;
2937
 
  }
2938
 
 
2939
 
  /// Create a new dof map for argument function i
2940
 
  virtual ufc::dof_map* create_dof_map(unsigned int i) const
2941
 
  {
2942
 
    return 0;
2943
 
  }
2944
 
 
2945
 
  /// Create a new cell integral on sub domain i
2946
 
  virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
2947
 
  {
2948
 
    return 0;
2949
 
  }
2950
 
 
2951
 
  /// Create a new exterior facet integral on sub domain i
2952
 
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
2953
 
  {
2954
 
    return 0;
2955
 
  }
2956
 
 
2957
 
  /// Create a new interior facet integral on sub domain i
2958
 
  virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
2959
 
  {
2960
 
    return 0;
2961
 
  }
2962
 
 
2963
 
};
2964
 
 
2965
 
#endif